c# 蒙特卡罗算法程序----对数正态分布的随机数的产生

  最近由于项目需求，需要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)函数就产生一个服从对数正态分布的随机数，根据模拟次数，可以产生多次随机数。如有不妥的地方，欢迎朋友们指正，如有其他更好的方法，也期待共享。