蒙特卡洛方法

1简述

蒙特·卡罗方法（Monte Carlo method），也称统计模拟方法，是二十世纪四十年代中期由于科学技术的发展和电子计算机的发明，而被提出的一种以概率统计理论为指导的一类非常重要的数值计算方法。是指使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。与它对应的是确定性算法。蒙特·卡罗方法在金融工程学，宏观经济学，计算物理学（如粒子输运计算、量子热力学计算、空气动力学计算）等领域应用广泛。

2基本思想

当所求解问题是某种随机事件出现的概率，或者是某个随机变量的期望值时，通过某种“实验”的方法，以这种事件出现的频率估计这一随机事件的概率，或者得到这个随机变量的某些数字特征，并将其作为问题的解。

3 蒙特卡洛方法解题步骤

（1）构造或描述概率过程

对于本身就具有随机性质的问题，如粒子输运问题，主要是正确描述和模拟这个概率过程，对于本来不是随机性质的确定性问题，比如计算定积分，就必须事先构造一个人为的概率过程，它的某些参量正好是所要求问题的解。即要将不具有随机性质的问题转化为随机性质的问题。

（2）实现从已知概率分布抽样

构造了概率模型以后，由于各种概率模型都可以看作是由各种各样的概率分布构成的，因此产生已知概率分布的随机变量（或随机向量），就成为实现蒙特卡罗方法模拟实验的基本手段，这也是蒙特卡罗方法被称为随机抽样的原因。最简单、最基本、最重要的一个概率分布是（0，1）上的均匀分布（或称矩形分布）。随机数就是具有这种均匀分布的随机变量。随机数序列就是具有这种分布的总体的一个简单子样，也就是一个具有这种分布的相互独立的随机变数序列。产生随机数的问题，就是从这个分布的抽样问题。在计算机上，可以用物理方法产生随机数，但价格昂贵，不能重复，使用不便。

另一种方法是用数学递推公式产生。这样产生的序列，与真正的随机数序列不同，所以称为伪随机数，或伪随机数序列。不过，经过多种统计检验表明，它与真正的随机数，或随机数序列具有相近的性质，因此可把它作为真正的随机数来使用。由已知分布随机抽样有各种方法，与从(0,1)上均匀分布抽样不同，这些方法都是借助于随机序列来实现的，也就是说，都是以产生随机数为前提的。由此可见，随机数是我们实现蒙特卡罗模拟的基本工具。

（3）建立各种估计量

一般说来，构造了概率模型并能从中抽样后，即实现模拟实验后，我们就要确定一个随机变量，作为所要求的问题的解，我们称它为无偏估计。建立各种估计量，相当于对模拟实验的结果进行考察和登记，从中得到问题的解。

4 蒙特卡洛的理论基础

蒙特卡洛方法的理论基础是大数定律。大数定律是描述相当多次数重复试验的结果的定律，根据这个定律知道
样本数量越多，其平均就越趋近于真实值。
首先考虑如下积分

接下来分别用蒙卡洛积分和牛顿莱布尼兹公式计算，在蒙特卡洛方法中样本很多时，它们的值应该相等。

利用蒙特卡洛方法，图像大致如下

5 什么是蒙特卡洛算法（通俗描述）

假如篮子里有1000个苹果，让你每次闭着眼睛找一个最大的，可以不限制挑选次数。于是，你可以闭着眼随机拿了一个，然后再随机拿一个与第一个比，留下大的，再随机拿一个，与前次留下的比较，又可以留下大的。循环往复这样，拿的次数越多，挑出最大苹果的可能性也就越大，但除非你把1000个苹果都挑一遍，否则你无法肯定最终挑出来的就是最大的一个。

6 box-muller

C++代码

#include <iostream>
#include <time.h>
#include <iomanip>
#include <math.h>
#include <fstream>
#define PI 3.14159
using namespace std;
void UNIFORM(double *);  //UINFORM函数声明
int x = 0;   //这里定义x一个全局变量并且初始付值0，这个的功用将会在子函数UNIFORM中得以体现
int main()
{
	int i, j;
	double A, B, C, E, D, r;
	double uni[2];
	double *p;
	srand((unsigned)time(NULL));  //随机数种子采用系统时钟
	ofstream outfile("Gauss.txt", ios::out);  //定义文件对象
	cout << "输入期望和方差:" << endl;
	cin >> E >> D;
	for (j = 0; j<1000; j++)
	{
		UNIFORM(&uni[0]);  //调用UNIFORM函数产生2个均匀分布的随机数并存入数组nui[2]
		A = sqrt((-2)*log(uni[0]));
		B = 2 * PI*uni[1];
		C = A*cos(B);
		r = E + C*D;    //E,D分别是期望和方差
		//outfile << r << "   ";  //将数据C输出到被定义的文件中
		cout << "r" << r << endl;
	}
	system("pause");
	return 0;
}
void UNIFORM(double *p)
{
	int i, a;
	double f;
	for (i = 0; i<2; i++, x = x + 689)
	{
		a = rand() + x;  //加上689是因为系统产生随机数的更换频率远远不及程序调用函数的时间
		a = a % 1000;
		f = (double)a;
		f = f / 1000.0;
		*p = f;
		p++;
	}
}

c#代码

using System;
using System.Collections.Generic;
using System.IO;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace ConsoleApp2
{
    class Program
    {
        int x = 0;
        double PI = 3.14159;
        static void Main(string[] args)
        {
            Program program = new Program();
            
            program.pds();
        }
        //srand((unsigned) time(NULL));
        Random rand = new Random();
        public void pds()
        {
            FileStream fs1 = new FileStream(@"C:\Users\Dell\Desktop\五院项目\test2.txt", FileMode.Create, FileAccess.Write);
            //创建写入文件 
            StreamWriter sw = new StreamWriter(fs1);

            int j =0, i;
            double A, B, C, r;
            double A1, B1, C1, r1;
            double E = 0, D = 1;
            int n1 = 10000;
            for(i =0;i< n1; i++)
            {
           
                A = Math.Sqrt((-2) * Math.Log(Gauss1()));
                B = 2 * PI * Gauss2();
                C = A * Math.Cos(B);
                r = E + C * D;

                sw.WriteLine(r.ToString("N6"));
 
                ++j;
            }
            sw.Close();
            fs1.Close();
            if (j == n1)
            {
                Console.WriteLine(n1 + "个数生成成功！");
            }
            else
            {
                Console.WriteLine("生成失败！");
            }
        }
        public double Gauss1()
        {

            double f,a;
            a = rand.Next() + x;
            a = a % 1000;
            f = (double)a;
            f = f / 1000.0;
            x = x + 689;
            return f;
        }
        public double Gauss2()
        {
            double f, a;
            a = rand.Next() + x;
            a = a % 1000;
            f = (double)a;
            f = f / 1000.0;
            return f;
        }
    
    }
}

算法--蒙特卡罗算法/正态分布随机数生成