蒙特卡罗方法介绍( 二)
一、蒙特卡罗求解定积分
蒙特卡洛方法求解定积分有两种方法,一种是上一节中讲的投点法,另外一种是期望法(也称平均值法)。
1.1 投点法
给出如下曲线f(x)f(x)f(x),求f(x)f(x)f(x)在a,ba,ba,b上的积分,即xxx在a,ba,ba,b上与曲线f(x)f(x)f(x)围成的面积。随机在如下矩形中进行投点,假设绿色点数为ggg,红色点数为rrr,矩阵的面积为AAA,则绿色部分的面积SSS就是:
S=A∗gr+gS = A*\frac{g}{r+g} S=A∗r+gg
注意,蒙特卡罗方法所求的数值是一个近似值,而不是精确值。
1.2 期望法(平均值法)
1.2.1 无意识统计学家定律(Law of the unconscious statistician,LOTUS)
统计学中期望描述了随机变量的平均取值的大小特征,计算期望过程:试验中每次可能结果的概率值乘以其结果的总和,也就是数值乘以其发生的概率的和。理解无意识统计学家定律就是告诉我们,事件万物相互之间是有相互联系的,一个事件的发生概率一定程度上可以用另一个事件的发生概率进行描述,比如函数F(X)=(X−3)2F(X) = (X-3)^2F(X)=(X−3)2中,XXX是一个随机变量,XXX的概率密度函数是影响到F(X)F(X)F(X)取值的概率密度函数的。
无意识统计学加定律,则是告诉已知随机变量XXX的概率分布fx(x)f_x(x)fx(x),但不知道g(X)g(X)g(X)的分布,此时用LOTUS公式能计算出函数g(X)g(X)g(X)的数学期望。
XXX是离散分布时
E[g(X)]=∑xg(X)fx(x)E[g(X)] = \sum_xg(X)f_x(x) E[g(X)]=x∑g(X)fx(x)
XXX是连续分布时
E[g(X)]=∫g(X)fx(x)dxE[g(X)] = \int g(X)f_x(x)dx E[g(X)]=∫g(X)fx(x)dx
常见的分布变换如下:
1.2.2 期望法
还是上边的例子,我们用期望法求解曲线的面积。
任意选取一组独立的、同分布的随机变量Xi{X_i}Xi,Xi{X_i}Xi在[a,b][a,b][a,b]上服从分布律fxf_xfx,也就是说fxf_xfx是随机变量XXX的概率密度函数(PDF)或者PMF。令:
g∗(x)=g(x)fx(x)g^*(x) = \frac{g(x)}{f_x(x)} g∗(x)=fx(x)g(x)
则g∗(x)g^*(x)g∗(x)也是一组独立同分布的随机变量。因为g∗(x)g^*(x)g∗(x)是xxx的函数,所以,根据LOTUS得:
E[g∗(Xi)]=∫abg∗(x)fx(x)dx=∫abg(x)dx=IE[g^*(X_i)] = \int_a^b g^*(x)f_x(x)dx = \int_a^b g(x)dx = I E[g∗(Xi)]=∫abg∗(x)fx(x)dx=∫abg(x)dx=I
由强大数定理知:
Pr(limN→∞1N∑g∗(Xi)=I)=1P_r(\lim_{N \to \infty} \frac{1}{N} \sum g^*(X_i) = I) = 1 Pr(N→∞limN1∑g∗(Xi)=I)=1
上式理解为所有概率值和为1.
假设若计算的积分有如下形式:
I=∫abg(x)dxI = \int_a^b g(x)dx I=∫abg(x)dx
其中,被积函数g(x)g(x)g(x)在区间[a,b][a,b][a,b]内可积,任意选择一个有简便方法可以抽样的函数fx(x)f_x(x)fx(x),满足以下条件:
{fx(x)≠0(a≤x≤b),g(x)≠0∫abg(x)dx=1\begin{cases} f_x(x) \neq 0& (a \le x \le b),g(x)\neq 0\\ \int_a^b g(x)dx=1& \\ \end{cases} {fx(x)=0∫abg(x)dx=1(a≤x≤b),g(x)=0
记:
g∗(x)={g(x)fx(x),fx(x)≠00,fx(x)=0g^*(x) = \begin{cases} \frac{g(x)}{f_x(x)} &,f_x(x)\neq 0\\ 0,& f_x(x) = 0\\ \end{cases} g∗(x)={fx(x)g(x)0,,fx(x)=0fx(x)=0
那么原积分可以写成:
I=∫abg∗(x)fx(x)dxI=\int_a^b g^*(x)f_x(x)dx I=∫abg∗(x)fx(x)dx
计算积分的步骤为:
- 产生服从分布fx(x)f_x(x)fx(x)的随机变量Xi,i=1,2,3...NX_i,{i=1,2,3...N}Xi,i=1,2,3...N
- 计算均值
I∗=1N∑1Ng∗(Xi)I^*=\frac{1}{N} \sum_1^N g^*(X_i) I∗=N11∑Ng∗(Xi)
将I∗I^*I∗作为III的近似值。
如果a,ba,ba,b为有限值,取fx(x)f_x(x)fx(x)为均匀分布
fx(x)={1b−a,a≤x≤b0,otherwisef_x(x) = \begin{cases} \frac{1}{b-a} &,a \le x \le b\\ 0,& otherwise\\ \end{cases} fx(x)={b−a10,,a≤x≤botherwise
则:
I=∫abg(x)dx=(b−a)∗1b−a∫abg(x)dx=(b−a)∫abg(x)1b−adxI=\int_a^b g(x)dx = (b-a) * \frac{1}{b-a}\int_a^b g(x)dx=(b-a)\int_a^b g(x)\frac{1}{b-a}dx I=∫abg(x)dx=(b−a)∗b−a1∫abg(x)dx=(b−a)∫abg(x)b−a1dx
计算积分的步骤为:
- 产生服从均匀分布fx(x)f_x(x)fx(x)的随机变量Xi,i=1,2,3...NX_i,{i=1,2,3...N}Xi,i=1,2,3...N
- 计算均值
I∗=b−aN∑1Ng(Xi)I^*=\frac{b-a}{N} \sum_1^N g(X_i) I∗=Nb−a1∑Ng(Xi)
将I∗I^*I∗作为III的近似值。
二、理解平均值法
给出如下问题,如图所示:已知曲线f(x)f(x)f(x),求
F=∫abf(x)dxF=\int_a^b f(x)dx F=∫abf(x)dx
以上问题,即求解阴影部分的面积。这里,我们使用一种近似的方法来求积分,任意选取a,ba,ba,b之间的一点xxx,函数值为f(x)f(x)f(x),则阴影部分的面积约等于:
S=f(x)∗(b−a)S=f(x)*(b-a) S=f(x)∗(b−a)
以上随机选取了一个点,如果随机选取4个点,那我们分别计算4个面积,然后取平均值,作为阴影部分的面积,如下图所示:
我们选取的点越多,计算出来的值也就越接近真实值。假设随机选取NNN个点,则根据以上,我们可以得出:
S=1N∑1Nf(xi)∗(b−a)=b−aN∑1Nf(xi)S=\frac{1}{N} \sum_1^N f(x_i)*(b-a) = \frac{b-a}{N} \sum_1^N f(x_i) S=N11∑Nf(xi)∗(b−a)=Nb−a1∑Nf(xi)
当然,这也是一种近似方法。
参考文献
- https://www.scratchapixel.com/lessons/mathematics-physics-for-computer-graphics/monte-carlo-methods-in-practice/monte-carlo-integration
- https://blog.csdn.net/qq_32618327/article/details/90452846
- https://blog.csdn.net/baimafujinji/article/details/53869358
- https://blog.csdn.net/qq_36931982/article/details/92788754