通俗易懂的Monte Carlo积分方法（一）

通俗易懂的Monte Carlo积分方法（一）

Monte Carlo积分的投点法计算：

• Monte Carlo算法（投点法）的数学基础：
  伯努利大数定律：
  $$\begin{gathered} 设f_A为n重伯努利试验中事件A发生的次数，而 \\ P为事件A发生的概率,那么\forall ϵ>0 \\ 有下面的式子成立： \end{gathered}$$
  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={\int }_a^{b}f(x)dx$$
  其中，而其中$0<f(x)<M$,由积分的面积定义我们可以得到：
  $${\int }_a^{b}f(x)dx=|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]\times [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∫ab​f(x)dx​
  对于$\mathbb{N}$个独立的均匀随机数${\xi }_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 的次数。

由大数定律可知,$\forall ϵ>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∫ab​f(x)dx​∣<ϵ}=1

即：
$${\int }_a^{b}f(x)dx\approx \frac{N_{s}M(b-a)}{N}(N足够大)$$

Python实现的代码：

例：求${\int }_0^2{x}^{2}dx$的值：

	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$的估计值。

在这里，我们取$\hat{I}=\frac{N_{s}M(b-a)}{N}$,由于 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=1n​xi​,xi​∈{0,1},xi​∼B(1,P)（二项分布）。而：
$$E\hat{I}={\int }_a^{b}f(x)dx$$
而对应的方差为：
$$Var\hat{I}=\frac{(b-a{)}^2{M}^{2}P(1-P)}{N}$$
对于允许误差$ϵ$,以及允许$\delta$的失败概率，确定$N_{ϵ,\delta }$,使得：
P { ∣ I ^ − ∫ a b f ( x ) d x ∣ > ϵ } < δ P\{{|\hat{I}-\int_a^bf(x)dx|>\epsilon}\}<\delta P{∣I^−∫ab​f(x)dx∣>ϵ}<δ
设${\varphi }_{\delta }$满足：
$$\delta =\frac{1}{\sqrt{2\pi }}{\int }_{{\varphi }_{\delta }}^{\infty }{e}^{-\frac{x^2}{2}}dx$$
由中心极限定理：
$$\hat{I}\sim N(E\hat{I},Var\hat{I})$$
可得：
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^−∫ab​f(x)dx∣​≥ϕ2δ​​}≤δ
如果我们选择$N$使得${\varphi }_{\frac{\delta }{2}}\sqrt{Var\hat{I}}\le ϵ$即可限制误差在我们预设的范围之内：

最后求得：
$$N>\frac{(b-a{)}^2{M}^2{\varphi }_{\frac{\delta }{2}}^2}{4{ϵ}^2}$$