首页 > 其他分享> > 蒙特卡洛原理

蒙特卡洛原理

2020-12-26 13:01:19 作者：互联网

蒙特卡洛方法的基本原理是，事件的概率可以用大量试验中发生的频率来估计，当样本容量足够大可以认为该事件的发生频率即为其概率．因此，可以先对影响其可靠度的随机变量进行大量的随机抽样，然后把这些抽样值一组一组地代入功能函数式，确定结构是否失效，最后从中求得结构的失效概率，蒙特卡洛法正是基于此思路进行分析的。蒙特卡洛方法在金融工程学、宏观经济学、计算、空气动力学计算）等领域应用广泛。
蒙特卡洛的基本做法是通过大量重复试验，通过统计频率，来估计概率，从而得到问题的求解。举个例子，如下图“蒙特卡洛法估计不规则图形面积”所示，一个矩形内有一个不规则图案，要求解不规则图形的面积。显然，矩形的面积可以简单计算为 $A = ab\times ac$ ，点位于不规则形状内的概率为 $A = \frac{A_{shape}}{A}$ 。现在重复往矩形范围内随机的投射点，样本点有一定概率会落在不规则图形内，若复杂n次试验，落到不规则图形内的次数为k，频率为k/n，若样本数量较大则有：

$A = \frac{A_{shape}}{A} \approx \frac{n}{k}$

1、硬币投掷实验

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                          %硬币投掷实验
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
p=0.5; %设置二项式分布binornd函数的概率
N=2000; %随机模拟的次数
sum=0;
for k=1:N
    sum=sum+binornd(1,p); %投掷实验的正面统计（假设硬币朝上为正面binornd返回1）
    P(k)=sum/k; %统计每次实验的硬币正面的概率
end
 
figure
hold on;box on;
plot(1:N,P);

N=10000的模拟结果：

2、几何概率模拟实验

对于一些比较复杂的几何图像，不好采用微积分方法计算图像的面积或参数，可以采用蒙特卡洛投点模拟方法计算出图像的面积或参数。假设两个人在1:00—2:00这段时间约会，他们约定先到者等20分钟后未等到另一个人则可以离开，试问他们两个人见面的概率有多大。这道题我们可以假设条件然后画出条件对应的图像，在根据图像的面积比计算两个人见面的概率。

假设： $x$ 和 $y$ 分别为两个人A、B到达的时间

A、B会面的条件 $={(x,y)|x-y|\leqslant 20,x<60,y<60}$ ,如下图所示中间的条纹区域就是两个人会面的时间交集区域，中间条纹这块面积为 $S=60\times 60-40\times 40=2000$

所以中间条纹区域占整个面积的比例 $P=S\div 60\times 60=2000\div 3600=5\div 9\approx 0.55555$

那么我们现在可以通过蒙特卡洛的方法来计算两个人会面的概率，我们通过生成随机数，然后统计落在条纹区域的随机数目再与总的随机数目对比，即得到两个人会面的概率，Matlab仿真程序如下：

%%%%%%%%%%%%%%%%%%%%%%%%%
  %投点发模拟两人会面的概率
%%%%%%%%%%%%%%%%%%%%%%%%%
 m=10000; %模拟的次数
success=0; %成功会面的次数
 j = 1;
 k = 1;
for i=1:m
     
    x=unifrnd(0,60); %生成x随机数
    y=unifrnd(0,60); %生成y随机数
    
    if abs(x-y)<=20 %两个人成功会面
        A1(j) = x; %A统计成功会面的坐标点
        A2(j) = y;
        success=success+1; %成功会面的次数
        j=j+1;
    end
    if abs(x-y)>20 %未成功会面
        B1(k) = x; %B统计未成功会面的坐标点
        B2(k) = y;
        k=k+1;
    end
end
 
P=success/m %成会面的概率
plot(A1,A2,'. r',B1,B2,'. b','MarkerSize',5); %画图表示

设置m=10000次的模拟，得出如下结论：

红色表示落在成功会面的区间，蓝色表示落在未成功会面的区间，红色点数与所有随机点数之比就得到了概率P，如下图所示：

3、蒲丰针实验-间接计算圆周率

18世纪蒲丰提出以下问题：设我们有一个等间距的平行线如图下图所示，随意抛一支长度比平行线之间距离小的针，求针和其中一条平行线相交的概率，并以此概率提出了一种计算圆周率的方法—随机投针法。

这一方法的步骤是：

1）取一张白纸，在上面画上许多条间距为a的平行线。

2）取一根长度为l（l≤a）的针，随机地向画有平行直线的纸上掷n次，观察针与直线相交的次数，记为m。

3）计算针与直线相交的概率。

假设随机投掷的一根针与平行线的夹角为 $\alpha$ ，夹角的范围 $[0,\pi ]$

针与任意一根平行线相交需满足的条件为 $\frac{x}{\sin\alpha }\leqslant \frac{l}{2}$ ，其中 $\frac{x}{\sin\alpha }$ 为针中心点位置到平行线的距离，l为针的长度，如下图所示：

把针投掷到许多条间距为a的平行线中，针中心点与平行线的距离为 $[0,\frac{a}{2}]$ ，a为平行线的间距，针是随机投掷到平面上，所以针与平行线距离的概率密度函数为 $f_{x}=\frac{2}{a}$

把针投掷到许多条间距为a的平行线中，针与平行线的夹角为 $[0,\pi ]$ ，所以针与平行线夹角 $\alpha$ 的概率密度函数为 $f_{\alpha }=\frac{1}{\pi }$

联合概率密度函数为 $F(x)=f_{x}\times f_{\alpha }=\frac{2}{a\pi }$ ，根据联合概率密度可求出针与平行线之间相交的概率为：

$p=\int_{0}^{\pi }\int_{0}^{\frac{l}{2}\sin \alpha }f_{x}\times f_{\alpha }dxd\omega=\int_{0}^{\pi }\frac{l\sin \alpha }{a\pi }d\omega=\frac{2l}{a\pi }$

从上面的公式可以看出，只要知概率P就可以求出圆周率 $\pi$ ，下面我们采用蒙特卡洛随机投点法求出概率P，从而求出圆周率。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 用蒙特卡洛方法计算圆周率π
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function  circumference_ratio
N=30; %重复试验的次数
needles=10000; %每次实验投掷针的数目
length=0.6;   %针的长度，平行线间距为1
PAI=zeros(1,N); %存储每次实验结果数值
 
for k=1:N
    
    PAI(k)=buffon(length,needles);   %求圆周率
    
end
 
PAI_ave=mean(PAI)   %对多次的实验取平均值
 
figure 
hold on;
box on;
plot(PAI);
 
line([0,N],[PAI_ave,PAI_ave],'LineWidth',3,'Color','g');
xlabel('k');
ylabel('π的估计值');
 
function pai=buffon(length,N)
frq=0; 
for k=1:N
  
    distance = unifrnd(0,0.5); %生成针中心点与平行线之间的间距
  
    angle = unifrnd(0,pi); %生成针与平行线的夹角
 
    if (distance <= (length*sin(angle)/2) ) %判断针是否与平行线相交
        frq=frq+1;   %统计相交的数量
    end    
end 
 
p=frq/N; %求取相交的概率 
pai=2*length/p; %计算圆周率

估计圆周率值：

多次实验结果及其均值，图中蓝色部分曲线是每次实验的数值，绿色部分是实验的平均值如下图所示:

4、蒙特卡洛定积分计算

假设要计算得定积分为： $I=\int_{b}^{a}g(x)dx$ ，该公式中a、b为有限值被积分函数为 $g(x)$ 为连续变量 $\xi$ 的概率密度函数，因此 $g(x)$ 必须满足一下条件：

(1) $g(x)$ 非负数

(2) $\int_{-\infty }^{+\infty }g(x)dx = 1$

从上述的公式中可以看出 $I$ 就是一个积分，并且满足 $I=P$ { $(a\leqslant \varepsilon \leqslant b)\bigcap \xi$ }

利用蒙特卡洛原理给定一个满足上述(1)(2)条件的概率密度函数 $g(x)$ ，如何求取出概率积分 $I$

(1)首先产生服从 $g(x)$ 分布的随机数m

(2)然后统计落在{ $(a\leqslant \varepsilon \leqslant b)\bigcap \xi$ }区域的随机数n

(3)求得落在{ $(a\leqslant \varepsilon \leqslant b)\bigcap \xi$ }区域的概率为 $p=\frac{n}{m}$

被积函数一般不满足公式(2)条件，所以需要对被积分函数进行变换，假设被积函数为 $y=\int_{a}^{b}f(x)dx$ 。

设定一定上限值M，设 $\Omega =[(x,f(x)),f(x)\leqslant M,a\leqslant x\leqslant b]$ ，如下图所示把一般的被积函数变换成概率密度函数，方形的面积为 $M(b-a)$ ，而函数 $f(x)$ 所包围阴影部分区域面积的值即为函数的积分值。

从上图可以得出 $\bg_white I=\int\int _{Y\leqslant M,a\leqslant x\leqslant b}\frac{1}{M(b-a)}dxdy = 1$ ，那么可以把被积分区域的联合概率密度函数写成 $pdf(x,y)=\frac{I_{f(x)\leqslant M,a\leqslant x\leqslant b}}{M(b-a)}$

根据上述可以求出阴影部分的概率 $\bg_white p(Y\leqslant f(x),a\leqslant x\leqslant b)$ ，如下所示：

$\bg_white I=p(Y\leqslant f(x),a\leqslant x\leqslant b)=\int_{Y\leqslant y(x)}\int\frac{I_{f(x)\leqslant M,a\leqslant x\leqslant b}}{M(b-a)}dxdy$

$=\int_{a}^{b}\int_{0}^{f(x)}\frac{1}{M(b-a)}dydx=\int_{a}^{b}\frac{f(x)}{M(b-a)}dx$ 令： $\Theta =\int_{a}^{b}f(x)dx$

所以： $I=p(Y\leqslant f(x))=\frac{\Theta }{M(b-a)}\approx \frac{n}{m}$ $\Rightarrow$ $\Theta \approx \frac{n}{m}M(b-a)$ ，从这个公式可以得出我们只要利用随机投点计算出n与m的比值就可以间接算出函数 $f(x)$ 的积分值，其中m为蒙特卡洛随机投点的数目，n为落在 $\Omega =[(x,f(x)),f(x)\leqslant M,a\leqslant x\leqslant b]$ 区域内的数目。

例子：求函数 $f(x) = -x^{2} + 5\times x + 8$ ，在[1，5]区间内的积分，通过积分运算很容易算出该函数的积分值为50.66666， $f(x)$ 函数曲线如下所示：

(1)又上图我们可知函数在[1，5]区间的最大值小于15，所以我们设定 $M=15$ ， $a=1,b=5$

(2) $\Theta =\int _{a}^{b} -x^{2} + 5\times x + 8 dx$ ，令 $\Theta = \frac{n}{m}15\times (5-1)$

(3)在 $\Omega =[(x,y),0\leqslant y\leqslant 15,1\leqslant x\leqslant 5]$ 区域内随机投掷m个点，并统计落在 $[(x,y),0\leqslant y\leqslant f(x),a\leqslant x\leqslant b]$ 区域内的随机点数n。

利用matlab编程如下：

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  利用蒙特卡洛计算定积分
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a=1;  %设定积分下限
b=5;  %设定积分上限
M=15; %设置M值
N=100000; %产生随机点的数目
freq=0; %统计落在函数区域内的随机点数

for i=1:N
  
    x=unifrnd(a,b); %产生[a,b]区间的随机数
    v=unifrnd(0,M); %产生[0,M]区间的随机数
 
    if (-power(x,2) + 5*x + 8) >= v %判断随机点数是否落在函数内部
        freq=freq+1;  %统计数目
    end
end

p=freq/N; %计算落在函数区域内部的概率
result=p*(b-a)*M %求出函数的积分值

计算出来的结果：

5、平均值法求取函数的积分值

假设一组同分布分布并且独立的随机变量{ $X _{i}$ }，随机变量{ $X _{i}$ }概率密度函数为 $pdf(x)$ ，令 $g(x)^{\ast } = \frac{g(x)}{pdf(x)}$ 那么{ $g^{\ast }(\xi _{i})$ }也是一组同分布并且相互独立的随机变量，现对函数 $g(x)$ 的

积分可以表示为 $\int _{a}^{b}g(x)dx = \frac{1}{n}\sum_{k=1}^{n}\frac{g(x_{k})}{pdf(x_{k})}$ ，其证明如下所示：

$E[G_{n}(X))] = E[{\frac{1}{n}}\sum_{k=1}^{n}\frac{g(X_{k})}{pdf(X_{k})}]$ $=\frac{1}{n}\sum_{k=1}^{n}\int \frac{g(x)}{pdf(x)}pdf(x)dx$

$=\frac{1}{n}\sum_{k=1}^{n}\int g(x)dx = \int g(x)dx$

选择抽样概率密度函数 $pdf(x)$ 需要满足以下条件

(1) 在区间[a，b]当 $g(x)\neq 0$ 时， $pdf(x)\neq 0$

(2)满足 $\int _{a}^{b}pdf(x)dx = 1$

例子：以上一个例题为例，求函数 $f(x) = -x^{2} + 5\times x + 8$ ，在[1，5]区间内的积分， $a=1,b=5$

在区间[a,b]之间独立同分布均匀采样， $pdf(x)=\frac{1}{b-a}=\frac{1}{4}$

所以 $g(x)^{\ast } = \frac{g(x)}{pdf(x)} = \frac{-x^{2} + 5\times x + 8}{\frac{1}{4}}$

积分值 $I= \frac{1}{n}\sum_{k=1}^{n}\frac{g(x_{k})}{pdf(x_{k})} = \frac{1}{n}\sum_{k=1}^{n}\frac{-x_{k}^{2} + 5\times x_{k} + 8}{\frac{1}{4}})$ ，利用MATLAB编写如下所示：

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 利用蒙特卡洛平均法计算定积分
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a=1; %积分下限 
b=5; %积分上限 
N=10000; %采样次数
sum=0; %求和
xrandnum = unifrnd(a,b,1,N); %利用unifrnd在区间[a,b]生成1xN矩阵的N个随机数
for ii=1:N
    %对函数-x^2 + 5x +8进行采样
    sum = sum - power(xrandnum(1,ii),2) + 5*xrandnum(1,ii) + 8;
end    
result = (b-a)*(sum/N) %求和计算计算积分值

结果如下所示，该结果与前面结算的结果一致，与真实值非常接近

标签：平行线,函数,积分,概率,随机,原理,蒙特卡洛
来源： https://blog.csdn.net/qq_34911636/article/details/111502109