蒙特卡洛积分的基本思想
考虑可积函数 g ( x ) g(x) g(x),现要计算 ∫ a b g ( x ) d x \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),则原积分形式可表示为 ∫ a b f ( x ) t ( x ) d x \int_a^bf(x)t(x)dx ∫abf(x)t(x)dx。记此时的x的密度函数为 f ( x ) f(x) f(x),那么即可表示为 E ( t ( x ) ) = ∫ a b f ( x ) t ( x ) d x E(t(x))=\int_a^bf(x)t(x)dx E(t(x))=∫abf(x)t(x)dx。不难看出, E ( t ( x ) ) E(t(x)) E(t(x))是很容易利用随机数进行估计的,它的值在样本量足够大的时候,可以用 ∑ i = 1 n t ( x i ) / n \sum\limits_{i=1}^nt(x_i)/n i=1∑nt(xi)/n进行估计(其中, x i x_i xi是来源于密度函数为 f ( x ) f(x) f(x)的样本)
以上思想的关键在于密度函数 f ( x ) f(x) f(x)的寻找(因为我们只有知道了x的分布,才能产生相对应的随机数),这就需要我们对一些基本分布的密度函数有所了解(例如:指数分布、正态分布、均匀分布等)。此外,在找 f ( x ) f(x) f(x)的时候还需注意积分区间是否一致。以 ∫ 0 ∞ e − x s i n x d x \int_0^{\infty}e^{-x}sinxdx ∫0∞e−xsinxdx为例,看到 e − x e^{-x} e−x可联想到指数分布的密度函数,而且此时的积分区间也恰好是0到正无穷,故是选参数为1的指数分布作为 f ( x ) f(x) f(x)。
以下在R语言中,以例子进行展示不同的方法计算一重积分和二重积分
一、内置函数+平均值法
例1 计算 ∫ 0 ∞ e − x s i n x d x \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 计算 ∬ D e ( x + y ) 2 d x d y \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 使用对偶变量法计算 ∫ 0 1 e − x s i n x d x \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) 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 使用对偶变量法计算 ∬ D e ( x + y ) 2 d x d y \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 计算积分 ∫ 0 1 e − x 1 + x 2 d x \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 使用控制变量法计算积分 ∬ D e ( x + y ) 2 d x d y \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 使用不同的重要性函数计算
∫
0
1
e
−
x
1
+
x
2
d
x
\int_0^{1}\frac{e^{-x}}{1+x^2}dx
∫011+x2e−xdx
考虑如下两个重要性函数:
f
1
(
x
)
=
e
−
1
(
1
−
e
−
1
)
−
1
f_1(x)=e^{-1}(1-e^{-1})^{-1}
f1(x)=e−1(1−e−1)−1
f 2 ( x ) = 4 π ( 1 + x 2 ) 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 使用重要抽样计算二重积分 ∬ D e ( x + y ) 2 d x d y \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%
以上就是本次分享的全部内容~