R语言实践利用概率统计方法计算三重积分

利用概率统计方法计算三重积分

一、蒙特卡罗尔平均值法

考虑K重积分 $I=\underset{\Omega }{\iiint }..\int f(x_1,x_2,...,x_k)d{x}_{1}d{x}_2..d{x}_k$ ，选取Ω范围内概率密度函数$g(x_1,x_2,x_3,...,x_k)$,则 $I$ 可以进行以下转换：
$$\begin{array}{rl}I& =\underset{\Omega }{\iiint }..\int \frac{f(x_1,x_2,x_3,...,x_k)}{g(x_1,x_2,x_3,...,x_k)}g(x_1,x_2,x_3,...,x_k)d{x}_{1}d{x}_{2}d{x}_3...d{x}_k\\ & =E(\frac{f(x_1,x_2,x_3,...,x_k)}{g(x_1,x_2,x_3,...,x_k)})\end{array}$$令$F(x_1,x_2,x_3,...,x_k)=\frac{f(x_1,x_2,x_3,...,x_k)}{g(x_1,x_2,x_3,...,x_k)}$ ，则
$I=E(F(x_1,x_2,x_3,...,x_k))$

根据辛钦大数定律,设随机变量X₁,X₂,X₃…X_k，相互独立服从同一分布，且具有数学期望$E(X_k)=\mu （k=1,2,3...)$,则对于任何正数 $\epsilon$，有： lim ⁡ n → ∞ P { ∣ 1 n ∑ k = 1 n X k − μ ∣ < ε } = 1 \lim\limits_{n\to ∞}P\{|\dfrac{1}{n}\sum\limits_{k=1}^{n}X_k-μ|<\varepsilon\}=1 n→∞lim​P{∣n1​k=1∑n​Xk​−μ∣<ε}=1

所以$I=E(F(x_1,x_2..,x_k))=\frac{1}{n}\sum _{i=1}^{n}F(x_{1i},x_{2i}...,x_{ki})$
简化： 取$g(x_1,x_2,x_3,...,x_k)=\frac{1}{V_{\Omega }}，(x_1,x_2,...,x_k)ϵ\Omega$ 服从K维均匀分布，则：
$$\begin{array}{rl}I& =E(F(x_1,x_2,x_3,...,x_k))\\ & =E(\frac{f(x_1,x_2,x_3,...,x_k)}{g(x_1,x_2,x_3,...,x_k)}\\ & =\frac{V_{\Omega }}{n}\sum _{i=1}^{n}f(x_{1i},x_{2i},x_{3i}...,x_{ki})\end{array}$$

二、R语言代码实践

题目： 请使用概率统计知识求解以下积分的近似值
$${\int }_1^3{\int }_2^4{\int }_3^6(x{e}^{-y^2+z^2}+xyz)dxdydz\approx 324\pm 2.9\ast 10^{-7}$$

#---------代码块---------
f=function(x,y,z) x*exp(-(y^2+z^2))+x*y*z
arr=array(1:1000)  #创建一维数组，长度1000
x1=3;x2=6
y1=2;y2=4
z1=1;z2=3
v=(x2-x1)*(y2-y1)*(z2-z1)  #空间体积
for(i in 1:1000)   #循环运算结果1000次，存入数组
{
	r1=runif(10^5,min=x1,max=x2)  #生成均匀分布随机数 x
	r2=runif(10^5,min=y1,max=y2)  #生成均匀分布随机数 y
	r3=runif(10^5,min=z1,max=z2)  #生成均匀分布随机数 z
	h=f(r1,r2,r3)
	out=v*mean(h)   #计算得到的积分值
	arr[i]=out	    #将积分值存入数组
}
print(mean(arr))    #输出计算1000次积分值后的平均值
#---------运行结果---------
[1] 324.0016

下面给出通用代码模板

fun=function(x) x[1]*exp(-x[2]^2-x[3]^2)+x[1]*x[2]*x[3]
start=c(1,2,3);end=c(3,4,6);n=10^5  #分别存放各个变量的起始值和结束值
simulation=function(fun,start,end,n)
{
	V=prod(end-start);      #(3-1)*(4-2)*(6-4)=12  求体积
	mat=matrix(0,nrow=n,ncol=length(start))  #n行3列矩阵，初值各元素为0
	for(i in 1:length(start)) #为每一列既3个变量生成随机数
	 	mat[,i]=runif(n,min=start[i],max=end[i])
	mean(apply(mat,1,fun))*V  #将矩阵的n行值代入fun中计算后取平均值最后乘以体积
}
simulation(fun,start,end,n)