蒙特卡罗方法 概率论本学期的大作业是关于蒙特卡罗积分，所以本文简单总结了蒙特卡罗积分。
概述
蒙特卡罗法也称统计模拟法、统计试验法，由美国数学家乌拉母(Ulam, S. M.)与美籍匈牙利数学家冯·诺依曼(von Neumann,J.)在20世纪40年代中叶提出。蒙特卡罗法是把概率现象作为研究对象的数值模拟方法，是按抽样调查法求取统计值来推定未知特性量的计算方法。该方法体现出了随机抽样的本质，利用计算机的快速计算能力，通过随机抽样（模拟）的过程，可以解决一些实际问题。我们早在初中就知道了抽样调查和普查的区别，正因为随机抽样的优点，它在实际的生产生活中被广泛运用。有时候我们无法做到对每个样本都进行研究，例如：工厂生产了一批灯泡，规定：只有当这批灯泡的合格率(合格即通电可正常使用，且寿命不低于设计寿命)大于等于0.9时，才可以出厂售卖。显然，我们对这批的每一个灯泡进行测试是不可能的，因为无论结果合格与否，这批灯泡都将报废。这时，我们必须随机的从这批灯泡中抽取一定数量的灯泡进行测试，抽取的灯泡的合格率就可以近似认为是这批灯泡的合格率，这样的方法，成本更低，结果准确。这与本文所讨论的蒙特卡罗法息息相关，即：按抽样调查法求取统计值来推定未知特性量的计算方法 。下面我通过一个具体的例子来更好的认识和理解蒙特卡罗方法。
利用蒙特卡罗方法解决简单的数学问题
估计π的值
π的定义为：圆的周长与直径的比值 。π值的计算有着很长的历史，随着计算方法的不断改进，精确度也越来越高，比如我们熟知的祖冲之，在480年，计算得π的近似值为3.1415926.首次将圆周率的计算精确到了小数点后7位。利用蒙特卡罗方法的思想，我们也可以利用计算机计算力高的优势，计算π的近似值。从圆面积的计算公式：S = π r 2 S=\pi r^2 S = π r 2 1 4 \frac14 4 1 1 4 \frac14 4 1 1 4 \frac14 4 1
假设“撒点”的总数为N，落在圆内的数量为n，正方形边长c=圆的半径r=1。则有：n N = 1 4 π r 2 c 2 (1-1)
\frac{n}{N}= \frac{\frac14 \pi r^2}{c^2} \tag{1-1}
N n = c 2 4 1 π r 2 ( 1-1 ) π = 4 n N (1-2)
\pi=\frac {4n}{N} \tag{1-2}
π = N 4 n ( 1-2 )
1 2 3 4 5 graph LR A(开始) --> 产生两个随机数x和y --> if{x,y到原点的距离小于1} if --->|yes| f1[n=n+1] --> B(结束) if --->|no| B(结束)
图 1 模拟撒点的pyhton算法流程图 对于不同的N值，我计算得到的π的近似值如下表所示。
表 1 不同N值下π的近似值
N
π的近似值
100
3.08
1000
3.156
10000
3.1516
100000
3.141
1000000
3.14712
10000000
3.141644
可见，当N不断增大时，π的结果也越来越准确。这个结果已经比较准确，所以这个方法估计π的近似值是靠谱的。
近似计算定积分
蒙特卡罗积分是蒙特卡罗方法在数学领域的一个重要的应用，下面讨论分析利用蒙特卡罗方法计算定积分的方法，并给出一般结论。在概率论中我们学过伯努利大数定理 ，它表明当试验次数足够大时，事件发生的频率依概率收敛到事件在一次试验中发生的概率，这便是上述计算π近似值方法的理论依据，也是下面讨论积分计算方法的重要理论支撑。
蒙特卡罗积分的理论推导及分析
假设我们要求定积分：∫ a b g ( x ) d x \int_a^b g(x)dx ∫ a b g ( x ) d x X 1 , X 2 , X 3 , ⋯ , X n X_1,X_2,X_3,\cdots,X_n X 1 , X 2 , X 3 , ⋯ , X n f ( x ) f(x) f ( x ) ∫ a b g ( x ) d x = ∫ a b g ( x ) f ( x ) f ( x ) d x (2-1)
\int_a^b g(x)dx=\int_a^b \frac{g(x)}{f(x)}f(x)dx \tag{2-1}
∫ a b g ( x ) d x = ∫ a b f ( x ) g ( x ) f ( x ) d x ( 2-1 ) X ∼ U ( a , b ) X\sim U(a,b) X ∼ U ( a , b ) f ( x ) = { 1 b − a x ∈ ( a , b ) 0 其他 (2-2)
f(x) =
\begin{cases}
\frac{1}{b-a} \qquad & x \in (a,b) \\
0 \qquad & 其他
\end{cases} \tag{2-2}
f ( x ) = { b − a 1 0 x ∈ ( a , b ) 其他 ( 2-2 ) ∫ a b g ( x ) f ( x ) f ( x ) d x = E [ g ( x ) f ( x ) ] (2-3)
\int_a^b \frac{g(x)}{f(x)}f(x)dx = E[\frac{g(x)}{f(x)}] \tag{2-3}
∫ a b f ( x ) g ( x ) f ( x ) d x = E [ f ( x ) g ( x ) ] ( 2-3 ) h ( x ) = g ( x ) f ( x ) h(x)=\frac{g(x)}{f(x)} h ( x ) = f ( x ) g ( x ) 辛钦大数定律 我们知道，当试验次数足够大时，样本的算数平均值依概率收敛于总体的均值。也就是说我们可以利用h ( x 1 ) , h ( x 2 ) , ⋯ , h ( x n ) h(x_1),h(x_2),\cdots,h(x_n) h ( x 1 ) , h ( x 2 ) , ⋯ , h ( x n ) E [ g ( x ) f ( x ) ] E[\frac{g(x)}{f(x)}] E [ f ( x ) g ( x ) ] h ‾ n = 1 n ∑ i = 1 n h ( x i ) = 1 n ∑ i = 1 n g ( x i ) f ( x i ) = E [ g ( x ) f ( x ) ] = ∫ a b g ( x ) d x (2-4)
\overline h_n = \frac1n \sum_{i=1}^n h(x_i) = \frac1n \sum_{i=1}^n \frac{g(x_i)}{f(x_i)} =E[\frac{g(x)}{f(x)}]=\int_a^b g(x)dx \tag{2-4}
h n = n 1 i = 1 ∑ n h ( x i ) = n 1 i = 1 ∑ n f ( x i ) g ( x i ) = E [ f ( x ) g ( x ) ] = ∫ a b g ( x ) d x ( 2-4 ) f ( x ) = 1 b − a f(x)=\frac{1}{b-a} f ( x ) = b − a 1 ∫ a b g ( x ) d x = b − a n ∑ i = 1 n g ( x i ) (2-5)
\int_a^b g(x)dx =\frac{b-a}{n} \sum_{i=1}^n g(x_i) \tag{2-5}
∫ a b g ( x ) d x = n b − a i = 1 ∑ n g ( x i ) ( 2-5 ) 上式便是利用蒙特卡罗法求解积分的近似解时的表达式，它的推导依赖于我们概率论中所学的大数定律。
对于式(2-4)，我们可以这样理解：( b − a ) g ( x i ) (b-a)g(x_i) ( b − a ) g ( x i )
随机数产生的原理
本文中的试验所涉及的随机数我都是利用python 中的random 库来产生的，而python 中random 库产生随机数的原理是利用了马特赛特旋转演算法。在计算机编程里生成的随机数是用确定性的算法计算出的一种在[0,1]上产生均匀分布 的随机数序列的算法，所以也叫“伪随机数”，但是在一般的统计意义上，我们仍可以视其为真随机数，所以，本文利用python 中random 库产生的随机数完全符合统计要求，做到了真随机，从而确保了蒙特卡罗方法的正确性。简言之，**利用python产生的随机数服从在指定区间上的均匀分布，且每个随机数之间相互独立。**这里就解释了推导式(2-2)时，让X服从均匀分布的原因：均匀分布计算简单且合理。
误差估计与分析
在数学上，我们一般用方差或者均方差来衡量误差程度。对于使用蒙特卡罗方法估算的积分，根据式(2-4)可知：D [ 1 n ∑ i = 1 n g ( x i ) f ( x i ) ] = 1 n 2 ∑ i = 1 n D [ g ( x i ) f ( x i ) ] = 1 n 2 ∑ i = 1 n ∫ a b [ g ( x i ) f ( x i ) ] 2 f ( x ) d x − E [ ( g ( x i ) f ( x i ) ) 2 ] = 1 n [ ∫ a b [ g ( x ) 2 f ( x ) ] d x − E [ ( g ( x ) f ( x ) ) 2 ] ]
\begin{equation}
\begin{split}
\qquad \qquad \qquad &D[\frac1n \sum_{i=1}^n \frac{g(x_i)}{f(x_i)}]\\
=&\frac 1{n^2} \sum_{i=1}^n D[\frac{g(x_i)}{f(x_i)}] \\
=&\frac 1{n^2} \sum_{i=1}^n \int_a^b[\frac{g(x_i)}{f(x_i)}]^2f(x)dx-E[(\frac{g(x_i)}{f(x_i)})^2]\\
=&\frac 1n[\int_a^b[\frac{g(x)^2}{f(x)}]dx-E[(\frac{g(x)}{f(x)})^2]]
\end{split}
\tag{2-6}
\end{equation}
= = = D [ n 1 i = 1 ∑ n f ( x i ) g ( x i ) ] n 2 1 i = 1 ∑ n D [ f ( x i ) g ( x i ) ] n 2 1 i = 1 ∑ n ∫ a b [ f ( x i ) g ( x i ) ] 2 f ( x ) d x − E [( f ( x i ) g ( x i ) ) 2 ] n 1 [ ∫ a b [ f ( x ) g ( x ) 2 ] d x − E [( f ( x ) g ( x ) ) 2 ]] ( 2-6 )
上述式子可以看出，估计值的标准差D ∝ 1 n \sqrt D \propto \frac{1}{\sqrt n} D ∝ n 1 1 n \frac{1}{\sqrt n} n 1
计算三重积分:∭ Ω ( x + e y + s i n z ) d x d y d z \iiint\limits_\Omega {(x+e^y+sinz)}\,{\rm d}x{\rm d}y{\rm d}z Ω ∭ ( x + e y + s in z ) d x d y d z
其中Ω \Omega Ω z = 8 − x 2 − y 2 z=8-x^2-y^2 z = 8 − x 2 − y 2 x 2 + y 2 = 4 x^2+y^2=4 x 2 + y 2 = 4 z = 0 z=0 z = 0
先求积分区域的体积，积分区域由一个圆柱体和一个以抛物面为边界的立体图形组成。通过圆柱体的体积计算公式和1的三重积分就是体积这一基本思想，不难算出积分区域的体积为24 π 24\pi 24 π
由重积分的对称性可知被积函数关于yoz平面对称，所以积分可以简化为：∭ Ω ( e y + s i n z ) d x d y d z \iiint\limits_\Omega {(e^y+sinz)}\,{\rm d}x{\rm d}y{\rm d}z Ω ∭ ( e y + s in z ) d x d y d z ∭ Ω ( e y + s i n z ) d x d y d z = 24 π 1 n ∑ i = 1 n ( e y i + s i n z i ) (3-1)
\iiint\limits_\Omega {(e^y+sinz)}\,{\rm d}x{\rm d}y{\rm d}z=24\pi \frac1n \sum_{i=1}^n (e^{y_i}+sinz_i) \tag{3-1}
Ω ∭ ( e y + s in z ) d x d y d z = 24 π n 1 i = 1 ∑ n ( e y i + s in z i ) ( 3-1 )
1 2 3 4 5 graph LR A(开始) --> input[/产生3个随机数x,y,z,/] --> if{x,y,z落在积分区域内} if --->|yes| f1[得到密度且n=n+1] --> B(结束) if --->|no| B
图 2 估算三重积分的pyhton算法流程图 对于不同的n值，计算得到的结果如下表所示。
表 2 不同n值下三重积分的近似值
n
三重积分的近似值
100
107.041
1000
122.030
10000
122.470
100000
122.063
1000000
121.573
前面说到，这个积分没有初等形式的解，我利用网络上的工具近似算得这个三重积分的值约为121.665. 可见，当n=1000000时，估算的值已经非常接近这个值了。
总结
蒙特卡罗方法的思想为我们解决一些正面难以入手的问题提供了一个很好的方法，在实际中，我们要求的面积或者多重积分往往形式复杂，难以解出其初等形式的解，利用蒙特卡罗方法和计算机模拟我们就可以很轻松的求得这类问题的解，给解决问题带来了很大的便利。从估算π的近似值到计算三重积分，深刻体会到了蒙特卡罗方法的核心思想和利用其求解问题的一般步骤。蒙特卡罗法的应用领域非常广泛，以上只是在求解积分领域的简单讨论。