布丰投针试验的仿真和误差估计

一.试验原理

1.试验步骤

​ 1.选一个长度为$l$的针。再选取一张白纸，上面划分许多平行且相隔为$a$的平行线。($l<a$)

​ 2.随机向纸投$n$次针，其中有$m$次针与平行线相交。

​ 3.计算投针的概率为$p=\frac{m}{n}$。

2.理论概率及$\pi$的计算与误差

​ 这根长为$l$的针的姿态由偏角$\theta$和针的中点到相邻平行线的最短距离$x$决定。且$x$和$\theta$的分布相互独立，其中$0<x<\frac{a}{2}$,而且其偏角$0<\theta <\pi$，二者均服从平面上的均匀分布。即：
$$\begin{gathered} x\sim U(0,\frac{a}{2}) \\ \theta \sim U(0,\pi ) \end{gathered}$$
​ 因此如果朝这个平面随机投针，那么$(x,\theta )$在平面上的二维均匀分布的概率密度就是：
$$f(x,y)=\{\begin{array}{ll}\frac{2}{\pi a}& 0<x<\frac{a}{2},0<\theta <\pi \\ 0& 其他\end{array}$$
​ 而要使得针与直线相交则要满足以下条件:
D : { ( x , θ ) : 0 ≤ x ≤ l 2 s i n θ } D:\{(x,\theta):0\leq x \leq \frac{l}{2}sin \theta\} D:{(x,θ):0≤x≤2l​sinθ}
​ 因此由二维连续随机变量概率计算公式可知:
P { D } = ∬ D f ( x , y ) d σ = ∬ 0 ≤ x ≤ l 2 s i n θ 2 π a d x d θ = 2 l π a P\{D\} = \iint_{D} f(x,y)d \sigma = \iint_{0 \leq x \leq \frac{l}{2}sin \theta}\frac{2}{\pi a}dxd \theta = \frac{2l}{\pi a} P{D}=∬D​f(x,y)dσ=∬0≤x≤2l​sinθ​πa2​dxdθ=πa2l​
​ 其中设定$0-1$变量$X$：
$$X=\{\begin{array}{ll}1& 针与平行线相交\\ 0& 针没有与平行线相交\end{array}$$
​ 我们做$n$次相互独立的试验,那么$X_i(i=1,2..n)$就服从以下参数的二项分布：
$$X_i\sim b(n,\frac{2l}{\pi a})$$
​ 我们可以知道每次伯努利试验的数学期望:$Ex=\frac{2l}{\pi a}$。同时做实验$n$次，得出由$m$次与直线相交：
$$\overline{X}=\frac{1}{n}\sum _{i=1}^n{X}_i=\frac{m}{n}$$
​ 因此，我们可以由中心极限定理的公式来估计误差:
l i m N → ∞ P { ∣ m n − 2 l π a ∣ < λ α σ n } = 2 2 π ∫ 0 λ α e − t 2 2 d t = 1 − α lim_{N \rightarrow \infty}P\{|\frac{m}{n}-\frac{2l}{\pi a}|<\frac{\lambda_{\alpha} \sigma}{\sqrt{n}}\} = \frac{2}{\sqrt{2\pi}}\int_0^{\lambda_{\alpha}}e^{-\frac{t^2}{2}}dt = 1-\alpha limN→∞​P{∣nm​−πa2l​∣<n ​λα​σ​}=2π ​2​∫0λα​​e−2t2​dt=1−α
​ 其中，当$\alpha =0.5,0.05,0.003$时，${\lambda }_{\alpha }=0.6745,1.96,3$。如果我们设$ϵ=\frac{{\lambda }_{\alpha }\sigma }{\sqrt{n}}$，那么最后$\pi$的$1-\alpha$的置信区间为：
$$[\frac{2l}{a(\frac{m}{n}+ϵ)},\frac{2l}{a(\frac{m}{n}-ϵ)}]$$

二.Python代码的实现

在这里我们取$\alpha =0.05$来用$MenteCarlo$法求$\pi$验证一下结果的正确性：

import random
import math
random.seed(5)
numbers = 1000000
m = 0
sum2 = 0
a = 4
l = 1
for i in range(numbers):
    x = a/2*random.random()
    theta = math.pi*random.random()
    if x<=(l/2*math.sin(theta)):
        m += 1
        sum2 += 1**2
Prop = m/numbers
pi_stand = 2*l/(a*Prop)
sigma = math.sqrt(1/numbers*sum2-(1/numbers*m)**2)
epsilon = 1.96*sigma /math.sqrt(numbers)
pi_max = 2*l/(a*(Prop-epsilon))
pi_min = 2*l/(a*(Prop+epsilon))
print('在{:}的试验之后，认为pi的估测值是{:.3f},有0.95的概率可以认为pi在[{:.4f} {:.4f}]'.format(numbers,pi_stand,pi_min,pi_max))

最后的结果是：

在1000000的试验之后，认为pi的估测值是3.145,有0.95的概率可以认为pi在[3.1314 3.1597]

Process finished with exit code 0

可以大概率认为我们的结果是正确的。