一.试验原理
1.试验步骤
1.选一个长度为 l l l的针。再选取一张白纸,上面划分许多平行且相隔为 a a a的平行线。( l < a l<a l<a)
2.随机向纸投 n n n次针,其中有 m m m次针与平行线相交。
3.计算投针的概率为 p = m n p = \frac{m}{n} p=nm。
2.理论概率及 π \pi π的计算与误差
这根长为
l
l
l的针的姿态由偏角
θ
\theta
θ和针的中点到相邻平行线的最短距离
x
x
x决定。且
x
x
x和
θ
\theta
θ的分布相互独立,其中
0
<
x
<
a
2
0<x<\frac{a}{2}
0<x<2a,而且其偏角
0
<
θ
<
π
0<\theta<\pi
0<θ<π,二者均服从平面上的均匀分布。即:
x
∼
U
(
0
,
a
2
)
θ
∼
U
(
0
,
π
)
x \sim U(0,\frac{a}{2})\\ \theta \sim U(0,\pi)
x∼U(0,2a)θ∼U(0,π)
因此如果朝这个平面随机投针,那么
(
x
,
θ
)
(x,\theta)
(x,θ)在平面上的二维均匀分布的概率密度就是:
f
(
x
,
y
)
=
{
2
π
a
0
<
x
<
a
2
,
0
<
θ
<
π
0
其
他
f(x,y) = \begin{cases} \frac{2}{\pi a} &\ 0<x<\frac{a}{2},0<\theta<\pi\\ 0&其他\end{cases}
f(x,y)={πa20 0<x<2a,0<θ<π其他
而要使得针与直线相交则要满足以下条件:
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≤2lsinθ}
因此由二维连续随机变量概率计算公式可知:
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}=∬Df(x,y)dσ=∬0≤x≤2lsinθπa2dxdθ=πa2l
其中设定
0
−
1
0-1
0−1变量
X
X
X:
X
=
{
1
针
与
平
行
线
相
交
0
针
没
有
与
平
行
线
相
交
X = \begin{cases} 1&针与平行线相交\\0 & 针没有与平行线相交\end{cases}
X={10针与平行线相交针没有与平行线相交
我们做
n
n
n次相互独立的试验,那么
X
i
(
i
=
1
,
2..
n
)
X_i(i=1,2..n)
Xi(i=1,2..n)就服从以下参数的二项分布:
X
i
∼
b
(
n
,
2
l
π
a
)
X_i \sim b(n,\frac{2l}{\pi a})
Xi∼b(n,πa2l)
我们可以知道每次伯努利试验的数学期望:
E
x
=
2
l
π
a
Ex = \frac{2l}{\pi a}
Ex=πa2l。同时做实验
n
n
n次,得出由
m
m
m次与直线相交:
X
‾
=
1
n
∑
i
=
1
n
X
i
=
m
n
\overline{X} = \frac{1}{n} \sum_{i=1}^n X_i = \frac{m}{n}
X=n1i=1∑nXi=nm
因此,我们可以由中心极限定理的公式来估计误差:
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−2t2dt=1−α
其中,当
α
=
0.5
,
0.05
,
0.003
\alpha = 0.5,0.05,0.003
α=0.5,0.05,0.003时,
λ
α
=
0.6745
,
1.96
,
3
\lambda_{\alpha} = 0.6745,1.96,3
λα=0.6745,1.96,3。如果我们设
ϵ
=
λ
α
σ
n
\epsilon=\frac{\lambda_{\alpha}\sigma}{\sqrt{n}}
ϵ=nλασ,那么最后
π
\pi
π的
1
−
α
1-\alpha
1−α的置信区间为:
[
2
l
a
(
m
n
+
ϵ
)
,
2
l
a
(
m
n
−
ϵ
)
]
[\frac{2l}{a(\frac{m}{n}+\epsilon)},\frac{2l}{a(\frac{m}{n}-\epsilon)}]
[a(nm+ϵ)2l,a(nm−ϵ)2l]
二.Python代码的实现
在这里我们取 α = 0.05 \alpha = 0.05 α=0.05来用 M e n t e C a r l o MenteCarlo 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
可以大概率认为我们的结果是正确的。