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

蒙特卡洛积分的基本思想

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

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

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

一、内置函数+平均值法

例1 计算${\int }_0^{\infty }{e}^{-x}sinxdx$

#内置函数
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 计算${\iint }_D{e}^{(x+y{)}^2}dxdy$，其中$D=[0,1]\times [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))$类型的积分。既然是估计，就有优劣之分。一般地，用方差来度量两个估计方法地优良性。定义如下：

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

二、对偶变量法

基本思想：

例1 使用对偶变量法计算${\int }_0^1{e}^{-x}sinxdx$，并与平均值法进行比较

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) else
      v <- 1 - u
    u <- 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 使用对偶变量法计算${\iint }_D{e}^{(x+y{)}^2}dxdy$，其中$D=[0,1]\times [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)$有很大的联系

例1 计算积分${\int }_0^1\frac{e^{-x}}{1+x^2}dx$，并与均值法比较

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 使用控制变量法计算积分${\iint }_D{e}^{(x+y{)}^2}dxdy$，其中$D=[0,1]\times [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 使用不同的重要性函数计算${\int }_0^1\frac{e^{-x}}{1+x^2}dx$
考虑如下两个重要性函数：
$$f_1(x)=e^{-1}(1-e^{-1}{)}^{-1}$$

$$f_2(x)=\frac{4}{\pi (1+x^2)}$$

以下代码中，在生成随机数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)$对减少估计量的方差尤为重要。

例2 使用重要抽样计算二重积分${\iint }_D{e}^{(x+y{)}^2}dxdy$，其中$D=[0,1]\times [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%

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