通俗易懂的Monte Carlo积分方法(一)
Monte Carlo积分的投点法计算:
- Monte Carlo算法(投点法)的数学基础:
伯努利大数定律:
设 f A 为 n 重 伯 努 利 试 验 中 事 件 A 发 生 的 次 数 , 而 P 为 事 件 A 发 生 的 概 率 , 那 么 ∀ ϵ > 0 有 下 面 的 式 子 成 立 : 设f_{A}为n重伯努利试验中事件A发生的次数, 而 \\P为事件A发生的概率,那么\forall\epsilon > 0\\有下面的式子成立: 设fA为n重伯努利试验中事件A发生的次数,而P为事件A发生的概率,那么∀ϵ>0有下面的式子成立:
l i m n → ∞ P { ∣ f A n − P ∣ < ϵ } = 1 lim_{n\rightarrow\infty}P\{|\frac{f_A}{n}-P|<\epsilon\}=1 limn→∞P{∣nfA−P∣<ϵ}=1
这个定律其实告诉我们的是当样本次数足够大的时候,频率是可以逼
近概率的。那么现在的关键问题是如何逼近?以及需要多大的样本才能准确的逼近?
-
Monte Carlo积分的计算:
1.解决如何逼近的问题:
首先,我们假设要计算的定积分是:
I = ∫ a b f ( x ) d x I=\int_a^bf(x)dx I=∫abf(x)dx
其中,而其中 0 < f ( x ) < M 0<f(x)<M 0<f(x)<M,由积分的面积定义我们可以得到:
∫ a b f ( x ) d x = ∣ S ∣ \int_a^bf(x)dx=|S| ∫abf(x)dx=∣S∣
其中: ∣ S ∣ |S| ∣S∣为 S = { ( x , y ) : a ⩽ x ⩽ b , f ( x ) > y } S=\{(x,y):a\leqslant x\leqslant b,f(x)>y \} S={(x,y):a⩽x⩽b,f(x)>y}的面积
考虑在平面区域 [ a , b ] × [ 0 , M ] [a,b]\times[0,M] [a,b]×[0,M]上的随机变量 ξ \xi ξ那么:
P { ξ ∈ S } = ∫ a b f ( x ) d x ( b − a ) M P\{\xi\in S\} = \frac{\int_a^bf(x)dx}{(b-a)M} P{ξ∈S}=(b−a)M∫abf(x)dx
对于 N \N N个独立的均匀随机数 ξ i , i = 1 , 2.. N \xi_i,i=1,2..N ξi,i=1,2..N,
记: N s 为 { ξ 1 , ξ 2 , . . . , ξ n } ∈ S N_s为\{\xi_1,\xi_2,...,\xi_n\} \in S Ns为{ξ1,ξ2,...,ξn}∈S 的次数。
由大数定律可知,
∀
ϵ
>
0
\forall\epsilon>0
∀ϵ>0:
l
i
m
N
→
∞
P
{
∣
N
s
N
−
∫
a
b
f
(
x
)
d
x
(
b
−
a
)
M
∣
<
ϵ
}
=
1
lim_{N\rightarrow\infty}P\{|\frac{N_s}{N}-\frac{\int_a^bf(x)dx}{(b-a)M}|<\epsilon\} = 1
limN→∞P{∣NNs−(b−a)M∫abf(x)dx∣<ϵ}=1
即:
∫
a
b
f
(
x
)
d
x
≈
N
s
M
(
b
−
a
)
N
(
N
足
够
大
)
\int_a^bf(x)dx \approx\frac{N_sM(b-a)}{N}(N足够大)
∫abf(x)dx≈NNsM(b−a)(N足够大)
Python实现的代码:
例:求 ∫ 0 2 x 2 d x \int_{0}^2x^2dx ∫02x2dx的值:
import random,time
start = time.perf_counter()
random.seed(10)
Num = 0
numbers = 100000
Num_s = 0
b = 2
a = 0
M = 2**2
for i in range(numbers):
x_value = (b-a)*random.random()
y_value = M*random.random()
if (y_value <= x_value**2)&(y_value>0):
Num_s += 1
Num += 1
P = Num_s/Num*M*(b-a)
print("\r目前的积分值是:{:.2f},时间是{:.5f}s".format(P,time.perf_counter()-start),end = ' ')
2.解决需要多大样本的问题:
中心极限定理估计样本个数:
中心极限定理就不赘述了,不懂得建议看教材《概率论与数理统计》。
在 N N N足够大的时候,中心极限定理可以得到 N N N的估计值。
在这里,我们取
I
^
=
N
s
M
(
b
−
a
)
N
\hat{I}=\frac{N_sM(b-a)}{N}
I^=NNsM(b−a),由于
N
s
=
∑
i
=
1
n
x
i
,
x
i
∈
{
0
,
1
}
,
x
i
∼
B
(
1
,
P
)
N_s=\sum_{i=1}^nx_i,x_i\in\{0,1\},x_i\sim B(1,P)
Ns=∑i=1nxi,xi∈{0,1},xi∼B(1,P)(二项分布)。而:
E
I
^
=
∫
a
b
f
(
x
)
d
x
E\hat{I} = \int_a^bf(x)dx
EI^=∫abf(x)dx
而对应的方差为:
V
a
r
I
^
=
(
b
−
a
)
2
M
2
P
(
1
−
P
)
N
Var\hat{I} = \frac{(b-a)^2M^2P(1-P)}{N}
VarI^=N(b−a)2M2P(1−P)
对于允许误差
ϵ
\epsilon
ϵ,以及允许
δ
\delta
δ的失败概率,确定
N
ϵ
,
δ
N_{\epsilon,\delta}
Nϵ,δ,使得:
P
{
∣
I
^
−
∫
a
b
f
(
x
)
d
x
∣
>
ϵ
}
<
δ
P\{{|\hat{I}-\int_a^bf(x)dx|>\epsilon}\}<\delta
P{∣I^−∫abf(x)dx∣>ϵ}<δ
设
ϕ
δ
\phi_{\delta}
ϕδ满足:
δ
=
1
2
π
∫
ϕ
δ
∞
e
−
x
2
2
d
x
\delta = \frac{1}{\sqrt{2\pi}}\int_{\phi_\delta}^{\infty}e^{-\frac{x^2}{2}}dx
δ=2π1∫ϕδ∞e−2x2dx
由中心极限定理:
I
^
∼
N
(
E
I
^
,
V
a
r
I
^
)
\hat{I}\sim N(E\hat{I},Var\hat{I})
I^∼N(EI^,VarI^)
可得:
P
{
∣
I
^
−
∫
a
b
f
(
x
)
d
x
∣
V
a
r
I
^
≥
ϕ
δ
2
}
≤
δ
P\{{\frac{|\hat{I}-\int_a^bf(x)dx|}{\sqrt{Var\hat{I}}}}\geq\phi_{\frac{\delta}{2}}\}\leq\delta
P{VarI^∣I^−∫abf(x)dx∣≥ϕ2δ}≤δ
如果我们选择
N
N
N使得
ϕ
δ
2
V
a
r
I
^
≤
ϵ
\phi_{\frac{\delta}{2}}\sqrt{Var\hat{I}}\leq\epsilon
ϕ2δVarI^≤ϵ即可限制误差在我们预设的范围之内:
最后求得:
N
>
(
b
−
a
)
2
M
2
ϕ
δ
2
2
4
ϵ
2
N>\frac{(b-a)^2M^2\phi_{\frac{\delta}{2}}^2}{4\epsilon^2}
N>4ϵ2(b−a)2M2ϕ2δ2