在概率计算中的应用¶
问题1¶
问题¶
设计随机试验求 $\pi$ 的近似值。
求解¶
在如图所示单位正方形中取 1000000 个随机点 $\left(x_{i}, y_{i}\right)$, $i=1,2, \cdots, 1000000$, 统计点落在 $x^{2}+y^{2} \leq 1$ 内的频数 $n$ 。则由儿何概率知, 任 取单位正方形内一点, 落在单位圆内部（第一象限部分）的概率为 $p=\frac{\pi}{4}$, 由于试验次数充分多, 频率近似于概率, 有 $\frac{n}{1000000} \approx \frac{\pi}{4}$, 所以 $\pi \approx \frac{4 n}{1000000}$ 。
代码
from numpy.random import rand
import numpy as np
N=1000000; x=rand(N); y=rand(N)
n=np.sum(x**2+y**2<1)
s=4*n/N; print(s)
问题2¶
问题¶
设计随机试验求 $\pi$ 的近似值。
下面设计另外一种类似的方法, 并用 Python 画图。
(1) 在边长为 $2 R$ 的正方形内随机取 $N$ 个点。
(2) 在正方形内画一个半径为 $R$ 的圆, 统计 $N$ 个点中在圆内点的个数 $n$ 。
(3) 圆的面积是 $\pi R^{2}$, 正方形的面积为 $(2 R)^{2}$, 因此二者的面积之比是 $\pi / 4$, 由几何概率的知识, 同样得到 $\pi \approx 4 n / N$ 。
代码
import numpy as np
import matplotlib.pyplot as plt
plt.rc('font',size=16); N=10000;
x,y=np.random.uniform(-1,1,size=(2,N))
inside=(x**2+y**2)<=1
mpi=inside.sum()*4/N #求pi的近似值
error=abs((mpi-np.pi)/np.pi)*100
outside=np.invert(inside)
plt.plot(x[inside],y[inside],'b.');plt.plot(x[outside],y[outside],'r.')
plt.plot(0,0,label='$\hat\pi$={:4.3f}\nerror={:4.3f}%'.format(mpi,error),alpha=0)
plt.axis('square'); plt.legend();
plt.savefig('images/sim0902.png')
plt.show()
问题3¶
问题¶
炮弹射击的目标为一椭圆 $\frac{x^{2}}{120^{2}}+\frac{y^{2}}{80^{2}}=1$ 所围成的区域的中心 当瞄准目标的中心发射时, 受到各种因素的影响, 炮弹着地点与目标中心有 随机偏差。设炮弹着地点围绕目标中心呈二维正态分布, 且偏差的标准差在 $x$ 和 $y$ 方向均为 100 米, 并相互独立, 用 Monte Carlo 法计算炮弹落在椭圆区域内的概率, 并与数值积分计算的概率进行比较。
求解¶
炮弹的落点为二维随机变量, 记为 $(X, Y),(X, Y)$ 的联合概率密度函 数为 $$ f(x, y)=\frac{1}{20000 \pi} e^{-\frac{x^{2}+y^{2}}{20000}} . $$
炮弹落在椭圆区域内的概率为 $$ p=\iint_{\frac{x^{2}}{120^{2}}+\frac{y^{2}}{80^{2}} \leq 1} \frac{1}{20000 \pi} e^{-\frac{x^{2}+y^{2}}{20000}} d x d y . $$
利用 Python 数值解的命令, 求得 $p=0.3754$ 。
我们也可以使用 Monte Carlo 法求概率。模拟发射了 $N$ 发炮弹，统计 炮弹落在椭圆 $\frac{x^{2}}{120^{2}}+\frac{y^{2}}{80^{2}}=1$ 内部的次数 $n$, 用炮弹落在椭圆内的频率近似所 求的概率, 模拟结果为所求的概率在 $0.3754$ 附近波动。
代码
import numpy as np
from scipy.integrate import dblquad
fxy=lambda x,y: 1/(20000*np.pi)*np.exp(-(x**2+y**2)/20000)
bdy=lambda x: 80*np.sqrt(1-x**2/120**2)
p1=dblquad(fxy,-120,120,lambda x:-bdy(x),bdy)
print("概率的数值解为：",p1)
N=1000000; mu=[0,0]; cov=10000*np.identity(2);
a=np.random.multivariate_normal(mu,cov,size=N)
n=((a[:,0]**2/120**2+a[:,1]**2/80**2)<=1).sum()
p2=n/N; print('概率的近似值为：',p2)
