R语言(二) 多种蒙特卡洛法计算一二重积分

蒙特卡洛积分的基本思想

考虑可积函数g(x)g(x)g(x)，现要计算∫abg(x)dx\int_a^bg(x)dx∫abg(x)dx。若能够将g(x)g(x)g(x)拆分为两个函数的乘积，即g(x)=f(x)t(x)g(x)=f(x)t(x)g(x)=f(x)t(x)，则原积分形式可表示为∫abf(x)t(x)dx\int_a^bf(x)t(x)dx∫abf(x)t(x)dx。记此时的x的密度函数为f(x)f(x)f(x)，那么即可表示为E(t(x))=∫abf(x)t(x)dxE(t(x))=\int_a^bf(x)t(x)dxE(t(x))=∫abf(x)t(x)dx。不难看出，E(t(x))E(t(x))E(t(x))是很容易利用随机数进行估计的，它的值在样本量足够大的时候，可以用∑i=1nt(xi)/n\sum\limits_{i=1}^nt(x_i)/ni=1∑nt(xi)/n进行估计（其中，xix_ixi是来源于密度函数为f(x)f(x)f(x)的样本）

以上思想的关键在于密度函数f(x)f(x)f(x)的寻找（因为我们只有知道了x的分布，才能产生相对应的随机数），这就需要我们对一些基本分布的密度函数有所了解(例如：指数分布、正态分布、均匀分布等)。此外，在找f(x)f(x)f(x)的时候还需注意积分区间是否一致。以∫0∞e−xsinxdx\int_0^{\infty}e^{-x}sinxdx∫0∞e−xsinxdx为例，看到e−xe^{-x}e−x可联想到指数分布的密度函数，而且此时的积分区间也恰好是0到正无穷，故是选参数为1的指数分布作为f(x)f(x)f(x)。

以下在R语言中，以例子进行展示不同的方法计算一重积分和二重积分

一、内置函数+平均值法

例1 计算∫0∞e−xsinxdx\int_0^{\infty}e^{-x}sinxdx∫0∞e−xsinxdx

#内置函数
f<-function(x){exp(-x)*sin(x)}
integrate(f,0,'inf')#生成1000个参数为1的指数分布随机数
r<-rexp(1000,1)
mean(sin(r))

从结果来看，两者计算结果还是十分接近的。内置函数integrate计算的积分值为0.5，而平均值法计算的积分值为0.4965115。

例2 计算∬De(x+y)2dxdy\iint_D{e^{(x+y)^2}}dxdy∬De(x+y)2dxdy，其中D=[0,1]×[0,1]D=[0,1]\times[0,1]D=[0,1]×[0,1]

#内置函数法
library(cubature)
double<-function(x){exp((x[1]+x[2])^2)}
adaptIntegrate(double,c(0,0),c(1,1))$integral#平均值法
n<-1000
f<-function(x,y){exp((x+y)^2)}
x<-runif(n)
y<-runif(n)
mean(f(x,y))

二重积分的蒙特卡洛法思想可由一重积分进行推广。从结果来看，该积分值约为4.8左右。

从上面两个例子不难看出，可以利用蒙特卡洛法积分法来估计E(g(X))E(g(X))E(g(X))类型的积分。既然是估计，就有优劣之分。一般地，用方差来度量两个估计方法地优良性。定义如下：

以下的三种方法都是为了减少估计量的方差而提出的。为更加客观比较，在此重复模拟1000次。

二、对偶变量法

基本思想：

例1 使用对偶变量法计算∫01e−xsinxdx\int_0^{1}e^{-x}sinxdx∫01e−xsinxdx，并与平均值法进行比较

n<-10000
k<-1000
set.seed(123)
MC.Phi <- function(n,k,anti = TRUE) {f<-function(x){exp(-x)*sin(x)}theat_hat <- numeric(length(k))for (i in 1:k) {u <- runif(n/2)if (!anti) v <- runif(n/2) elsev <- 1 - uu <- c(u, v)theat_hat[i]<-mean(f(u))}theat_hat
}
mean(MC.Phi(n,k))  #对偶变量积分
mean(MC.Phi(n,k,anti=FALSE)) #平均值法积分
print((1-var(MC.Phi(n,k))/var(MC.Phi(n,k,anti = FALSE)))*100) #方差减少百分比

计算结果显示，使用对偶变量法计算的积分与均值法基本无差，但估计量的方差减少了61.93%

例2 使用对偶变量法计算∬De(x+y)2dxdy\iint_D{e^{(x+y)^2}}dxdy∬De(x+y)2dxdy，其中D=[0,1]×[0,1]D=[0,1]\times[0,1]D=[0,1]×[0,1]，并与均值法进行比较

n<-10000
k<-1000
set.seed(123)
f<-function(x,y){exp((x+y)^2)}
#均值法
simple_carlo<-function(n,k){theat_hat<-numeric()for(i in 1:k){  x<-runif(n)y<-runif(n)theat_hat[i]<-mean(f(x,y))}return(theat_hat)
}
theat_simple<-simple_carlo(n,k)
mean(theat_simple)#对偶变量法
MC.phi<-function(n,k){theat_hat=numeric()for(i in 1:k){x<-runif(n/4)y<-runif(n/4)T1<-f(x,y)T2<-f(x,1-y)T3<-f(1-x,y)T4<-f(1-x,1-y)theat_hat[i]<-mean((T1+T2+T3+T4)/4)}return(theat_hat)
}
theat_MC<-MC.phi(n,k)
mean(theat_MC)
print((1-var(theat_MC)/var(theat_simple))*100)

从对二重积分的计算结果来看，对偶变量法依旧比平均值法减少方差约56.30%。

三、控制变量法

基本思想:

注意：减小的方差百分比与选择的相关函数f(x)f(x)f(x)有很大的联系

例1 计算积分∫01e−x1+x2dx\int_0^{1}\frac{e^{-x}}{1+x^2}dx∫011+x2e−xdx，并与均值法比较

n<-10000
k<-1000
set.seed(123)
g<-function(x){exp(-x)/(1+x^2)}#均值法
simple_carlo<-function(n,k){theat_hat<-numeric()for(i in 1:k){x<-runif(n)theat_hat[i]<-mean(g(x))}return(theat_hat)
}
theat_simple<-simple_carlo(n,k)
mean(theat_simple)#控制变量法
control_carlo<-function(n,k){f <- function(x){exp(-.5)/(1+x^2)}theat_hat<-numeric()for(i in 1:k){x<-runif(n)B <- f(x)A <- g(x)a <- -cov(A,B) / var(B) #est of c*T1 <- g(x)theat_hat[i] <- mean(T1 + a * (f(x) - exp(-.5)*pi/4))}theat_hat
}
theat_control<-control_carlo(n,k)
mean(theat_control)
print((1-var(theat_control)/var(theat_simple))*100)

计算结果来看，控制变量法的减少方差十分明显呀，比均值法减少了近95%的方差。

例2 使用控制变量法计算积分∬De(x+y)2dxdy\iint_D{e^{(x+y)^2}}dxdy∬De(x+y)2dxdy，其中D=[0,1]×[0,1]D=[0,1]\times[0,1]D=[0,1]×[0,1]，并与均值法作比较

#均值法与前面相同,在此不再写出
n<-10000
k<-1000
set.seed(123)
f<-function(x,y){exp((x+y)^2)}
control_carlo<-function(n,k){theat_hat<-numeric()for(i in 1:k){u<-runif(n)v<-runif(n)g<-function(u,v){exp(u+v)}r=-cov(g(u,v),f(u,v))/var(g(u,v)) #估计c*T1<-f(u,v)theat_hat[i]<-mean(T1+r*(g(u,v)-(exp(1)-1)^2))}return(theat_hat)
}
theat_control<-control_carlo(n,k)
mean(theat_control)
print((1-var(theat_control)/var(theat_simple))*100)

积分值约为4.9，使用控制变量法比均值法的方差减少约79.28%

四、重要性抽样法

基本思想：

注意：减少的方差百分比依旧与选择重要性函数有关

例1 使用不同的重要性函数计算∫01e−x1+x2dx\int_0^{1}\frac{e^{-x}}{1+x^2}dx∫011+x2e−xdx
考虑如下两个重要性函数：
f1(x)=e−1(1−e−1)−1f_1(x)=e^{-1}(1-e^{-1})^{-1}f1(x)=e−1(1−e−1)−1

f2(x)=4π(1+x2)f_2(x)=\frac{4}{\pi(1+x^2)}f2(x)=π(1+x2)4

以下代码中，在生成随机数x时，用到了随机数的逆变换法。

n<-10000
k<-1000
set.seed(123)
g<-function(x){exp(-x)/(1+x^2)}
simple_carlo<-function(n,k){theat_hat<-numeric()for(i in 1:k){x<-runif(n)theat_hat[i]<-mean(g(x))}return(theat_hat)
}
theat_simple<-simple_carlo(n,k)
mean(theat_simple) #均值法#重要性函数f1
f1<-function(n,k){theat_hat<-numeric()for(i in 1:k){u<-runif(n)x<-log(1-u*(1-exp(-1)))  #逆变换法fg<-g(x)/(exp(-x)/(1-exp(-1)))theat_hat[i]<-mean(fg)}theat_hat
}
theat_f1<-f1(n,k)
mean(theat_f1)
print((1-var(theat_f1)/var(theat_simple))*100)#重要性函数f2
f2<-function(n,k){theat_hat<-numeric()for(i in 1:k){u<-runif(n)x<-tan(pi*u/4)  #逆变换法fg<-g(x)/(4/(1+x^2)*pi)theat_hat[i]<-mean(fg)}theat_hat
}
theat_f2<-f2(n,k)
mean(theat_f2)
print((1-var(theat_f2)/var(theat_simple))*100)

三种方法估计的积分值基本都在0.52左右，但是使用f1方法比均值法的方差减少约83.75%，而使用f2方法只减少了63.06%。由此可见，寻找一个合适的f(x)f(x)f(x)对减少估计量的方差尤为重要。

例2 使用重要抽样计算二重积分∬De(x+y)2dxdy\iint_D{e^{(x+y)^2}}dxdy∬De(x+y)2dxdy，其中D=[0,1]×[0,1]D=[0,1]\times[0,1]D=[0,1]×[0,1]，并与均值法进行比较

#重要抽样
n<-10000
k<-1000
set.seed(123)
f<-function(x,y){exp((x+y)^2)}
r_num<-function(n){u<-runif(n)x<-log((exp(1)-1)*u+1)return(x)
}
importance_sample<-function(n,k){g<-function(x,y){exp(x+y)/(exp(1)-1)^2}theat_hat<-numeric()for(i in 1:k){x<-r_num(n)y<-r_num(n)fg<-f(x,y)/g(x,y)theat_hat[i]<-mean(fg)}return(theat_hat)
}
theat_imp<-importance_sample(n,k)
mean(theat_imp)
print((1-var(theat_imp)/var(theat_simple))*100)

从结果来看，使用重要抽样比均值法的方差减少约71.16%

以上就是本次分享的全部内容~