最近由于项目需求，需要c#编程实现蒙特卡罗算法。在网上找了好几天的资料，都没找到自己想要的结果，最终还是得靠自己动手哦。
关于蒙特卡罗算法的概念，意义及具体步骤介绍，可以查阅相关的资料，都有详细介绍。我这里大概分三步：
1·确定参数的一个数据分布和决定模拟次数，在服从该分布的情况下产生N个随机数。
2·将每次产生的随机数，带入表达式中算出结果，这样可以得到N个结果。
3·用统计方法把N个结果的数字特征表示出来，比如，求出平均值，方差等。
对于c#程序的实现我觉得难点在于第一步，也就是概率随机数的产生，我这里也主要讲的是随机数的产生，以服从对数正态分布的随机数为例。
在Matlab和EXCEL中都有自带的库函数，可以直接产生。c#中则要自己编写，本来想自己是数学专业毕业的，对于这个问题应该不难，但是我错了，反反复复写了几次都不对，衡量的标准是公司一兄弟用excel中的loginv函数产生的随机数来比较。在两研究生同事的帮忙推数学公式的情况下，研究了excel中的loginv函数后，终于代码落成，由于代码中使用很多数学公式，而在代码中不方便注释，请读者查看对数正态分布和正态分布相关函数以及EXCEL的帮助文档。代码在visual studio 2010中测试通过，没问题。
/// <summary>
/// 产生对数正态分布随机数
/// </summary>
/// <param name="miu">平均值（实际应用中常需要将使用Math.Log(miu)作为参数）</param>
/// <param name="sigma">方差（实际应用中常需要将使用0.5*Math.Log(sigma)作为参数）</param>
/// <returns></returns>
public static double Random_LogNorMal(double miu, double sigma)
{
var rand = new Random(Guid.NewGuid().GetHashCode());
var x = rand.NextDouble();
return Math.Exp(miu + sigma * NormalFun(x));
}
/// <summary>
/// 误差函数估算
/// </summary>
/// <param name="x"></param>
/// <returns></returns>
private static double ErfFun(double x)
{
double d = 2.0 / Math.Sqrt(Math.PI);
double t = x / Math.Sqrt(2);
double dsum = 0;
double dn = 0;
for (int n = 0; ; n++)
{
if (n == 0)
{
dn = t;
}
else
{
dn = Math.Pow(-1, n) * (Math.Pow(t, 2 * n + 1) / ((2 * n + 1) * nby(n)));
}
if (Math.Abs(dn) < 0.00001)
{
break;
}
dsum += dn;
}
return 0.5 * (d * dsum + 1);
}
/// <summary>
/// n阶乘算法
/// </summary>
/// <param name="n"></param>
/// <returns></returns>
private static double nby(int n)
{
double by = 1;
for (int i = 1; i <= n; i++)
{
by = by * i;
}
return by;
}
/// <summary>
/// 正态分布随机数产生（初次以1为迭代单位，尔后采用二分查找法搜索符合要求随机数）
/// </summary>
/// <param name="p"></param>
/// <returns></returns>
private static double NormalFun(double p)
{
double z = 0;
double sumz = 0.5;
if (p > 0.5)
{
do
{
z += 1;
sumz = ErfFun(z);
}
while (sumz < p);
if (sumz - p > 0.00001) //可根据需求自行确定精度。
{
return SearchNum(z - 1, z, p);
}
else return z;
}
else
if (p < 0.5)
{
do
{
z -= 1;
sumz = ErfFun(z);
}
while (sumz > p);
if (p - sumz > 0.00001)
{
return SearchNum(z, z + 1, p);
}
else return z;
}
return z;
}
/// <summary>
/// 递归搜索随机数
/// </summary>
/// <param name="z1"></param>
/// <param name="z2"></param>
/// <param name="p"></param>
/// <returns></returns>
private static double SearchNum(double z1, double z2, double p)//z1>z2;
{
double z = (z1 + z2) / 2;
double ef = ErfFun(z);
if (ef < p && p - ef > 0.00001)
{
return SearchNum(z, z2, p);
}
else
{
if (ef > p && ef - p > 0.00001)
{
return SearchNum(z1, z, p);
}
}
return z;
}
需要注意的是Random_LogNorMal(double miu, double sigma)中的miu和sigma并不是随机数的平均值和方差，而是随机数的对数平均值和方差，所以常用miu=Math.Log(平均值)，sigma=0.5*Math.Log(方差)作为参数传递，如有不明处可以参照对数正态分布和正态分布相关介绍。
调用一次Random_LogNorMal(double miu, double sigma)函数就产生一个服从对数正态分布的随机数，根据模拟次数，可以产生多次随机数。如有不妥的地方，欢迎朋友们指正，如有其他更好的方法，也期待共享。
