蒙特卡洛方法计算三重积分

蒙特卡洛方法计算三重积分

蒙特卡洛方法简介

通过大量随机样本，了解一个系统或模拟一个过程,进而得到所要计算的值。
eg:已知正方形内部有一内切圆，计算圆周率π 解：正方形与圆的面积之比是π/ 4
现在在正方形内部买，随机产生10000个点，计算他们与中心点的距离，从而判断是否落在圆的内部。如果这些点均匀分布，那么圆内的点占到所有点的π/4。随着点的数量增加，模拟的结果越接近理论值。

蒙特卡洛算法实现步骤

蒙特卡洛算法应用在金融工程学，宏观经济学，生物医学，计算物理学等领域应用广泛。

应用蒙特卡洛方法解决实际问题的时候主要有两部分工作：1.模拟过程，产生各种随机变量。2.表征模型的数字特征，从而得到问题的近似解。

计算三重积分
使用数学方法计算：
使用Python 的scipy包的积分功能计算

使用蒙特卡洛方法计算：
使用Python 的 numpy包生成随机数并进行计算
误差分析
分别比较点数N取不同值时的误差，因为程序每次运行的结果不同，取同一N值的10次运行结果取平均数计算误差。

数学方法

from scipy import integrate
import scipy
def f(z, y, x): 
   return scipy.sin(z) + x + scipy.exp(y)

result1 = integrate.tplquad(f, -2, 2, lambda x: (-1) *
                           scipy.sqrt(4 - x * x), lambda x: scipy.sqrt(4 - x * x), 0, lambda x, y: 8 - x * x - y * y)
print(result1)
## (121.66509988033158, 7.739990467710454e-07)

蒙特卡洛方法

import numpy
N = 10000000
accum = 0
z_sum = 0
y_sum = 0
for i in range(N):
   x = numpy.random.uniform(-2, 2)
   y_max = numpy.sqrt(4 - x * x)
   y_min = (-1) * y_max
   y_sum += y_max
   y = numpy.random.uniform(y_min, y_max)
   z_max = 8 - x * x - y * y
   z_sum += z_max
   z = numpy.random.uniform(0, z_max)
   accum += numpy.sin(z) + x + numpy.exp(y)
z_mean = z_sum / N
y_mean = y_sum / N
measure = (2 - (-2)) * 2 * y_mean * (z_mean-0)
result2 =  measure * accum/float(N)
print(result2)
## 119.15345432219415
## 计算误差
rate = numpy.abs(result1 - result2) / result1
print(rate)

误差分析

• N = 1000时， 误差为2.82%
• N = 10000时，误差为2.34%
• N = 100000时，误差为1.76%
• N = 1000000时，误差为1.58%
• N = 10000000时，误差为1.55%

随机数产生原理
查阅numpy.random.uniform的官方文档。可知uniform产生随机数的分布是均匀分布。符合蒙特卡洛算法的要求。

4.蒙特卡洛算法的拓展与后续问题

由蒙特卡洛算法近似计算积分，联系到华—王方法近似计算高维积分，后续可以进一步学习。 numpy.random包中还提供了别的更复杂的概率模型买，适用于更多更复杂的场景。后续可以查看文档进行相关学习。

5.学习-参考资料

1.https://numpy.org/doc/stable/release.html
2.https://risk-engineering.org/notebook/monte-carlo-integration.html
3.https://www.investopedia.com/terms/m/montecarlosimulation.asp
4.https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.tplquad.html
5.www.ibm.com/cloud/learn/monte-carlo-simulation
6.https://theclevermachine.wordpress.com/tag/monte-carlo-integration/