Monte Carlo随机模拟

原文链接：http://blog.csdn.net/on2way/article/details/47280671

Monte Carlo方法也叫随机模拟、随机抽样或者统计实验方法。其主要用途用于模拟一些无法用数值产生的随机系统。比如当系统的各个单元的特征量已知，但系统过于复杂导致无法预测其准确数学模型，这个时候可以用随机模拟法计算系统的相关参数。

蒙特卡罗的基本思想，为了求解数学、物理、生产管理等方面的问题，首先需要建立一个概率模型或者一个抽样试验来计算所求参数的统计特征，最后通过计算机模拟求出近似解。

我们说蒙特卡罗模拟的是随机现象，那么不同的问题的随机现象是不同的。比如掷骰子，各个点数等概率出现，模拟的均匀分布。去某个地方排队，人来的过程是一个随时间陆续变化的，这种一般认为符合指数分布。再有正态分布的例子太多，一般的统计数据没特殊要求都认为是服从正态分布的。每个分布都有对应的函数来表示。

假设一个随机事件在一次实验中发生的概率为P(x),显然0≤p(x)≤1

• 一次贝努利实验：

P(ξ=x)={p,1−p,x=1x=0

• . 二项分布 
  当把贝努利实验重复多做几次就是二项分布了 
  
  P(X=ak)=Cknpk(1−p)n−k,(p<1)
• 离散均匀分布 
  以相同的概率取所有可能的值 
  
  P(X=ak)=1/n,k=1,2,3...n
• 泊松分布 
  发生概率低，同时次数无限增加的贝努利分布就是泊松分布 
  
  P(X=k)=λkk!e−λ,k=1,2,...,n

再看一些密度函数

• 均匀分布密度函数 
  在区间[a,b]上等概率出现，那么出现的概率： 
  
  P(x)=⎧⎩⎨1b−a0,x∈[a,b]other
• 正态分布密度函数 
  用的最多 
  
  P(x)=12π−−√σexp(−(x−u)22σ2)
• 指数分布密度函数 
  
  P(x)={λe−λx0x≥0x<0

下面采用matlab进行蒙特卡罗的一些小实验 
首先关于matlab产生上述具有特定分布的随机数都有特定的函数，常用的有三种：

• rand：生成均匀分布的伪随机数。分布在（0~1）之间 
  主要语法：rand(m,n)生成m行n列的均匀分布的伪随机数 
  rand(m,n,’double’)生成指定精度的均匀分布的伪随机数，参 数还可以是’single’, rand(RandStream,m,n)利用指定的RandStream(我理解为随 机种子)生成伪随机数
• randn 生成标准正态分布的伪随机数（均值为0，方差为1） 
  主要语法：和上面一样
• randi 生成均匀分布的伪随机整数 
  主要语法：randi（iMax）在开区间（0，iMax）生成均匀分布的伪随机整数randi（iMax，m，n）在开区间（0，iMax）生成m*n型随机矩 阵r = randi([iMin,iMax],m,n)在开区间（iMin，iMax）生成 m*n型随机矩阵

还有一些其他的生成函数： 
unifrnd()：和rand()类似，这个函数生成某个区间内均匀分布的随机数 
normrnd()：和randn()类似，此函数生成指定均值、标准差的正态分布的随机数 
chi2rnd()：函数生成服从卡方（Chi-square)分布的随机数。卡方分布只有一个参数：自由度v 
frnd()：函数生成服从F分布的随机数。F分布有2个参数：v1, v2 
trnd()：函数生成服从t分布的随机数。t分布有1个参数：自由度v 
betarnd()：此函数生成服从Beta分布的随机数。Beta分布有两个参数分别是A和B 
exprnd()：此函数生成服从指数分布的随机数。指数分布只有一个参数: mu 
gamrnd()：生成服从Gamma分布的随机数。Gamma分布有两个参数：A和B

例如：用蒙塔卡罗方法计算图像的面积，这里为了方便起见计算一下圆的面积，顺便用这种方法预估一下π的值。当然这是规则图像，可以用公式来计算，当遇到的图像不是规则的图形的时候，就非常有用了。 
代码如下：

clcclear%生成num个（-1,1）的均匀点num = 100000;data = 2*rand(2,num)-1;plot(data(1,:),data(2,:),'*')hold on;%num_in_circle = 0;radius = 1;%小于1%plot circletheta = 0:0.1:2*pi;plot(radius*cos(theta),radius*sin(theta),'r','LineWidth',3);%test samplesfor i=1:num    %随机点在圆内    if (data(1,i)^2 + data(2,i)^2) <= radius^2        num_in_circle = num_in_circle + 1;    endendarea = 4*num_in_circle/num;%矩形面积为4PI = area/(radius^2);%title(['num=',num2str(num),' 半径=',num2str(radius),' 模拟的PI值 = ',num2str(PI)])
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

 
 
 
 
可以看到num样本越多预测的越准，当再大些时，基本上可以和π相差不多了。