Monte Carlo Integration and Variance Reduction
计算期望的数学方法:
E
[
g
(
X
)
]
=
∫
−
∞
∞
g
(
x
)
f
(
x
)
d
x
E[g(X)]=\int_{-\infin}^{\infin}g(x)f(x)dx
E[g(X)]=∫−∞∞g(x)f(x)dx
Simple Monte Carlo estimator 简单MC估计量
对于单位区间上的积分:
θ
^
=
∫
0
1
g
(
x
)
d
x
\hat \theta=\int_0^1g(x)dx
θ^=∫01g(x)dx
可以有以下估计量。
θ
^
=
g
m
(
X
)
‾
=
1
m
g
(
X
i
)
\hat \theta=\overline{g_m(X)}=\frac{1}{m}g(X_i)
θ^=gm(X)=m1g(Xi)
由强大数定律,\hat \theta会以1的概率收敛到
E
[
g
(
X
)
]
E[g(X)]
E[g(X)]
而这个期望值就是\theta。
用r代码比较精确解的和MonteCarlo解。
m<-10000
x<-runif(m)
theta.hat<-mean(exp(-x))
print(theta.hat)
print(1-exp(-1))
推广1:对于积分区间在[a,b]上的情况,可以用换元法。
y = t − a b − a t = ( b − a ) y + a y=\frac{t-a}{b-a}\\t=(b-a)y+a y=b−at−at=(b−a)y+a
∫ a b g ( t ) d t = ( b − a ) ∫ 0 1 g ( ( b − a ) y + a ) d y \int_a^bg(t)dt=(b-a)\int_0^1g((b-a)y+a)dy ∫abg(t)dt=(b−a)∫01g((b−a)y+a)dy
也可以用其它方式生成[a,b]区间上的随机数进行积分。比如
∫
a
b
g
(
t
)
d
t
=
(
b
−
a
)
∫
a
b
g
(
t
)
1
b
−
a
d
t
\int_a^bg(t)dt=(b-a)\int_a^bg(t)\frac{1}{b-a}dt
∫abg(t)dt=(b−a)∫abg(t)b−a1dt
右式的积分便可以看成是对[a,b]上均匀分布的随机变量求MC积分,所以
θ
^
=
(
b
−
a
)
g
(
X
)
‾
\hat \theta=(b-a)\overline{g(X)}
θ^=(b−a)g(X)
m<-10000
x<-runif(m,min=2,max=4)
theta.hat<-mean(exp(-x))*2
print(theta.hat)
print(exp(-2)-exp(-4))
用R代码比较精确解和MC解。
推广2:无界区间上的MC积分
Φ ( x ) = ∫ − ∞ x 1 2 π e − t 2 / 2 d t \Phi(x)=\int_{-\infin}^x\frac{1}{\sqrt{2\pi}}e^{-t^2/2}dt Φ(x)=∫−∞x2π1e−t2/2dt
对于无界的情况,不能直接套公式。但是,根据对称性,可以有
Φ
(
x
)
=
∫
−
∞
0
1
2
π
e
−
t
2
/
2
d
t
+
∫
0
x
1
2
π
e
−
t
2
/
2
d
t
=
0.5
+
∫
0
x
1
2
π
e
−
t
2
/
2
d
t
\Phi(x)=\int_{-\infin}^0\frac{1}{\sqrt{2\pi}}e^{-t^2/2}dt+\int_{0}^x\frac{1}{\sqrt{2\pi}}e^{-t^2/2}dt\\ =0.5+\int_{0}^x\frac{1}{\sqrt{2\pi}}e^{-t^2/2}dt
Φ(x)=∫−∞02π1e−t2/2dt+∫0x2π1e−t2/2dt=0.5+∫0x2π1e−t2/2dt
又回到了有界区间积分的情况。如果x小于0,则可以考察
Φ
(
x
)
=
1
−
Φ
(
−
x
)
\Phi(x)=1-\Phi(-x)
Φ(x)=1−Φ(−x)
再根据-x>0归入下面的情况。
对于如下积分,可以用[0,x]区间上的均匀分布随机数来做MC积分。
∫
0
x
1
2
π
e
−
t
2
/
2
d
t
=
x
∫
0
x
1
2
π
e
−
t
2
/
2
1
x
d
t
\int_{0}^x\frac{1}{\sqrt{2\pi}}e^{-t^2/2}dt=x\int_{0}^x\frac{1}{\sqrt{2\pi}}e^{-t^2/2}\frac{1}{x}dt
∫0x2π1e−t2/2dt=x∫0x2π1e−t2/2x1dt
但是,这样做每次都要生成一组x,不能复用。所以,可以考虑换元法.
y
=
t
x
∈
[
0
,
1
]
,
t
=
x
y
y=\frac{t}{x}\in[0,1],t=xy
y=xt∈[0,1],t=xy
∫
−
∞
x
1
2
π
e
−
t
2
/
2
d
t
=
∫
0
1
1
2
π
e
−
(
x
y
)
2
/
2
x
d
y
\int_{-\infin}^x\frac{1}{\sqrt{2\pi}}e^{-t^2/2}dt=\int_0^1\frac{1}{\sqrt{2\pi}}e^{-(xy)^2/2}xdy
∫−∞x2π1e−t2/2dt=∫012π1e−(xy)2/2xdy
现在,可以写出MC积分式(暂时略去常系数)
θ
=
E
Y
[
x
e
−
(
x
Y
)
2
/
2
]
θ
^
=
g
m
(
y
)
‾
=
1
m
∑
i
=
1
m
x
e
−
(
y
i
x
)
2
/
2
\theta=E_Y\bigg[xe^{-(xY)^2/2}\bigg]\\ \hat \theta=\overline {g_m(y)}=\frac{1}{m}\sum_{i=1}^mxe^{-(y_ix)^2/2}
θ=EY[xe−(xY)2/2]θ^=gm(y)=m1i=1∑mxe−(yix)2/2
x<-seq(.1,2.5,length=10)
m<-10000
u<-runif(m)
cdf<-numeric(length(x))
for(i in 1:length(x)){
g<-x[i]*exp(-(u*x[i])^2/2)
cdf[i]<-mean(g)/sqrt(2*pi)+0.5
}
Phi<-pnorm(x)
print(round(rbind(x,cdf,Phi),3))
R代码用MC求标准正态分布的CDF。
HIt or Miss Algorithm
E [ I ( Z ≤ x ) ] = P ( Z ≤ x ) = Φ ( x ) E[I(Z\le x)]=P(Z\le x)=\Phi(x) E[I(Z≤x)]=P(Z≤x)=Φ(x)
由此,可以用样本均值作为上式的估计量
Φ
(
x
)
^
=
1
m
∑
i
=
1
m
I
(
z
i
≤
x
)
\widehat{\Phi(x)}=\frac{1}{m}\sum_{i=1}^mI(z_i\le x)
Φ(x)
=m1i=1∑mI(zi≤x)
用R代码检验Hit or Miss 算法
x<-seq(.1,2.5,length=10)
m<-10000
z<-rnorm(m)
dim(x)<-length(x)
p<-apply(x,MARGIN=1,FUN=function(x,z){mean(z<x)},z=z)
Phi<-pnorm(x)
print(round(rbind(x,p,Phi),3))
可以看出,在中心区域估计的效果较差,在尾部估计的效果较好。(加一个0就基本上都好了。)
简单MC积分的总结
对于如下形式的积分,
θ
=
∫
A
g
(
x
)
f
(
x
)
d
x
,其中
f
(
x
)
≥
0
,
∫
A
f
(
x
)
d
x
=
1
\theta=\int_Ag(x)f(x)dx,其中f(x)\ge0,\int_Af(x)dx=1
θ=∫Ag(x)f(x)dx,其中f(x)≥0,∫Af(x)dx=1
可以先按照分布f(x)生成一组样本x1,x2,…xm,再求R.V.的函数的样本均值。
θ
^
=
1
m
∑
i
=
1
m
g
(
x
i
)
\hat \theta=\frac{1}{m}\sum_{i=1}^mg(x_i)
θ^=m1i=1∑mg(xi)
由强大数定律,它会收敛到\theta。
MC积分的标准差
对于
σ
2
=
V
a
r
f
(
g
(
X
)
)
\sigma^2=Var_f(g(X))
σ2=Varf(g(X))
我们可以写出
V
a
r
(
θ
^
)
=
σ
2
m
Var(\hat \theta)=\frac{\sigma^2}{m}
Var(θ^)=mσ2
然而,在估计的时候,我们并不能真正获得X的分布,因此只能靠样本来估计\sigma^2。因此,\hat \theta 的方差也只能是估计。
σ
^
2
m
=
1
m
2
∑
i
=
1
m
[
g
(
x
i
)
−
g
(
x
)
‾
]
2
\frac{\hat \sigma^2}{m}=\frac{1}{m^2}\sum_{i=1}^m[g(x_i)-\overline{g(x)}]^2
mσ^2=m21i=1∑m[g(xi)−g(x)]2
右侧是用了plug-in估计的g(X)的方差。
由此,可以写出估计量的标准差
s
e
^
(
θ
^
)
=
σ
^
m
=
1
m
{
∑
i
=
1
m
[
g
(
x
i
)
−
g
(
x
)
‾
]
2
}
1
/
2
\widehat {se}(\hat \theta)=\frac{\hat \sigma}{\sqrt{m}}=\frac{1}{m}\bigg\{\sum_{i=1}^m[g(x_i)-\overline{g(x)}]^2\bigg\}^{1/2}
se
(θ^)=mσ^=m1{i=1∑m[g(xi)−g(x)]2}1/2
这可以用在中心极限定理中。
θ
^
−
E
[
θ
^
]
s
e
^
(
θ
^
)
∼
N
(
0
,
1
)
a
s
m
→
∞
\frac{\hat \theta-E[\hat \theta]}{\widehat{se}(\hat \theta)}\sim N(0,1)\; as\; m\rightarrow\infin
se
(θ^)θ^−E[θ^]∼N(0,1)asm→∞
接近于正态的\hat \theta可以用于MC积分的误差界分析中。
用R代码利用标准差估计MC积分的误差界
x<-2
m<-10000
z<-rnorm(m)
g<-(z<x) #the indicator function
v<-mean((g-mean(g))^2)/m
cdf<-mean(g)
c(cdf,v)
c(cdf-1.96*sqrt(v),cdf+1.96*sqrt(v))
Variance and Efficiency 方差和效率
V a r ( g ( X ) ‾ ) = 1 m V a r ( g ( X ) ) Var(\overline{g(X)})=\frac{1}{m}Var(g(X)) Var(g(X))=m1Var(g(X))
E ( g ( X ) ‾ ) = θ b − a E(\overline{g(X)})=\frac{\theta}{b-a} E(g(X))=b−aθ
所以,
V
a
r
(
θ
^
)
=
(
b
−
a
)
2
V
a
r
(
g
(
X
)
‾
)
=
(
b
−
a
)
2
m
V
a
r
(
g
(
X
)
)
Var(\hat \theta)=(b-a)^2Var(\overline{g(X)})=\frac{(b-a)^2}{m}Var(g(X))
Var(θ^)=(b−a)2Var(g(X))=m(b−a)2Var(g(X))
由中心极限定理,当m很大时,\overline(g(X))近似于正态分布,所以\hat \theta也近似于正态分布,方差即由上式给出。
对于hit-or miss方法,有另外的方差表达式:
F
(
x
)
^
=
g
(
X
)
‾
=
1
m
∑
i
=
1
m
I
(
X
i
≤
x
)
\widehat{F(x)}=\overline{g(X)}=\frac{1}{m}\sum_{i=1}^mI(X_i\le x)
F(x)
=g(X)=m1i=1∑mI(Xi≤x)
V a r ( F ( x ) ^ ) = p ( 1 − p ) / m = F ( x ) ( 1 − F ( x ) ) / m s e ^ ( F ( x ) ^ ) = { F ^ ( x ) ( 1 − F ^ ( x ) ) / m } Var(\widehat {F(x)})=p(1-p)/m=F(x)(1-F(x))/m\\ \widehat{se}(\widehat{F(x)})=\sqrt{\bigg\{\widehat F(x)(1-\widehat F(x))/m\bigg\}} Var(F(x) )=p(1−p)/m=F(x)(1−F(x))/mse (F(x) )={F (x)(1−F (x))/m}
Efficiency
一个估计量比另一个估计量更有效:
V
a
r
(
θ
1
^
)
V
a
r
(
θ
2
^
)
<
1
\frac{Var(\hat{\theta_1})}{Var(\hat{\theta_2})}<1
Var(θ2^)Var(θ1^)<1
考虑到总可以通过增大样本量的方式减少方差,效率还是需要重视的。
Variance Reduction 减小方差
我们来考虑多元随机变量的函数:
X
(
j
)
=
{
X
1
(
j
)
,
X
2
(
j
)
,
.
.
.
,
X
n
(
j
)
}
X^{(j)}=\{X_1^{(j)},X_2^{(j)},...,X_n^{(j)}\}
X(j)={X1(j),X2(j),...,Xn(j)}
考虑下面的函数
Y
j
=
g
(
X
1
(
j
)
,
X
2
(
j
)
,
.
.
.
,
X
n
(
j
)
)
,
j
=
1
,
2
,
.
.
.
m
Y_j=g(X_1^{(j)},X_2^{(j)},...,X_n^{(j)}),j=1,2,...m
Yj=g(X1(j),X2(j),...,Xn(j)),j=1,2,...m
Y1,…,Ym就是独立同分布的函数,
E
[
Y
ˉ
]
=
E
[
1
m
∑
j
=
1
m
Y
j
]
=
θ
E[\bar Y]=E\bigg[\frac{1}{m}\sum_{j=1}^mY_j\bigg]=\theta
E[Yˉ]=E[m1j=1∑mYj]=θ
方差就有
V
a
r
(
θ
^
)
=
V
a
r
Y
ˉ
=
V
a
r
f
g
(
X
)
m
Var(\hat \theta)=Var\bar Y=\frac{Var_fg(X)}{m}
Var(θ^)=VarYˉ=mVarfg(X)
增大采样数目可以减小MC估计量的方差。但是很不划算。一般来说,m的下界由下式给出
m
≥
⌈
σ
2
/
e
2
⌉
m
≥
⌊
σ
2
/
e
2
⌋
m\ge\lceil\sigma^2/e^2 \rceil\\ m\ge\lfloor\sigma^2/e^2 \rfloor
m≥⌈σ2/e2⌉m≥⌊σ2/e2⌋
其中sigma^2是目标随机变量的方差,e则是对于估计量我们所能接受的标准差。
Antithetic Variables Approach 对偶变量方法
考虑两个变量和的方差:
V
a
r
(
U
1
+
U
2
2
)
=
1
4
(
V
a
r
(
U
1
)
+
V
a
r
(
U
2
)
+
2
C
o
v
(
U
1
,
U
2
)
)
Var(\frac{U_1+U_2}{2})=\frac{1}{4}(Var(U_1)+Var(U_2)+2Cov(U_1,U_2))
Var(2U1+U2)=41(Var(U1)+Var(U2)+2Cov(U1,U2))
当U1,U2独立时,方差式比较简单。当U1,U2 负相关时,方差可以减小。
This fact leads us to consider negtively correlated variables as a possible method for reducing variance.
考虑均匀分布:
U
j
∼
U
i
n
f
(
0
,
1
)
U_j\sim Uinf(0,1)
Uj∼Uinf(0,1)
则有
1
−
U
j
∼
U
n
i
f
(
0
,
1
)
1-U_j\sim Unif(0,1)
1−Uj∼Unif(0,1)
且以上两组随机变量是负相关的。由此,我们考虑用逆变换:
Y
j
=
g
(
F
X
−
1
(
U
1
(
j
)
)
,
.
.
.
,
F
X
−
1
(
U
n
(
j
)
)
)
Y
j
′
=
g
(
F
X
−
1
(
1
−
U
1
(
j
)
)
,
.
.
.
,
F
X
−
1
(
1
−
U
n
(
j
)
)
)
Y_j=g(F_X^{-1}(U_1^{(j)}),...,F_X^{-1}(U_n^{(j)}))\\ Y_j'=g(F_X^{-1}(1-U_1^{(j)}),...,F_X^{-1}(1-U_n^{(j)}))
Yj=g(FX−1(U1(j)),...,FX−1(Un(j)))Yj′=g(FX−1(1−U1(j)),...,FX−1(1−Un(j)))
什么时候Y和Y’负相关呢?当g为单调函数的时候,以上两个随机变量负相关。
PROPOSITION:若X1,…Xn是独立的,f,g具有相同的单调性,那么
E
[
f
(
X
)
g
(
X
)
]
≥
E
[
f
(
X
)
]
E
[
g
(
X
)
]
E[f(X)g(X)]\ge E[f(X)]E[g(X)]
E[f(X)g(X)]≥E[f(X)]E[g(X)]
证明:
因为f,g具有相同的单调性,所以,对任意x,y
(
f
(
x
)
−
f
(
y
)
)
(
g
(
x
)
−
g
(
y
)
)
≥
0
E
[
(
f
(
X
)
−
f
(
Y
)
)
(
g
(
X
)
−
g
(
Y
)
)
]
≥
0
(f(x)-f(y))(g(x)-g(y))\ge 0\\ E[(f(X)-f(Y))(g(X)-g(Y))]\ge 0
(f(x)−f(y))(g(x)−g(y))≥0E[(f(X)−f(Y))(g(X)−g(Y))]≥0
因此,
E
[
f
(
X
)
g
(
X
)
]
+
E
[
f
(
Y
)
g
(
Y
)
]
≥
E
[
f
(
X
)
g
(
Y
)
]
+
E
[
f
(
Y
)
g
(
X
)
]
E[f(X)g(X)]+E[f(Y)g(Y)]\ge E[f(X)g(Y)]+E[f(Y)g(X)]
E[f(X)g(X)]+E[f(Y)g(Y)]≥E[f(X)g(Y)]+E[f(Y)g(X)]
因为X和Y是独立同分布的,所以上式可以改写为
KaTeX parse error: No such environment: eqnarray at position 8: \begin{̲e̲q̲n̲a̲r̲r̲a̲y̲}̲ 2E[f(X)g(X)]&&…
由此,上述陈述对n=1时成立。考虑数学归纳法,对n元的情况,先写出条件分布
KaTeX parse error: No such environment: eqnarray at position 8: \begin{̲e̲q̲n̲a̲r̲r̲a̲y̲}̲ E[f(X)g(X)|X_n…
再对条件分布(关于X_n的一元函数)取期望,得到
$$
\begin{eqnarray}
E[f(X)g(X)]&&=&&E[E[f(X)g(X)|X_n=n]]\
&&\ge&&E[E[f(X)|X_n]E[g(X)|X_n]]\
&&\ge&&E[E[f(X)|X_n]]E[E[g(X)|X_n]](归纳的基础,一元的情况)\
&&=&&E[f(X)]E[f(X)]\
\end{eqnarray}
$$
现在,我们就可以证明一开始想要的那个结论了(Y和Y‘负相关)
推论:If g=g(X1,X2,…Xn)是单调函数,那么
Y
j
=
g
(
F
X
−
1
(
U
1
(
j
)
)
,
.
.
.
,
F
X
−
1
(
U
n
(
j
)
)
)
Y
j
′
=
g
(
F
X
−
1
(
1
−
U
1
(
j
)
)
,
.
.
.
,
F
X
−
1
(
1
−
U
n
(
j
)
)
)
Y_j=g(F_X^{-1}(U_1^{(j)}),...,F_X^{-1}(U_n^{(j)}))\\ Y_j'=g(F_X^{-1}(1-U_1^{(j)}),...,F_X^{-1}(1-U_n^{(j)}))
Yj=g(FX−1(U1(j)),...,FX−1(Un(j)))Yj′=g(FX−1(1−U1(j)),...,FX−1(1−Un(j)))
负相关。
证明:
不失一般性,认为g是单调非减函数,那么
Y
j
=
g
(
F
X
−
1
(
U
1
(
j
)
)
,
.
.
.
,
F
X
−
1
(
U
n
(
j
)
)
)
−
Y
j
′
=
−
g
(
F
X
−
1
(
1
−
U
1
(
j
)
)
,
.
.
.
,
F
X
−
1
(
1
−
U
n
(
j
)
)
)
Y_j=g(F_X^{-1}(U_1^{(j)}),...,F_X^{-1}(U_n^{(j)}))\\ -Y_j'=-g(F_X^{-1}(1-U_1^{(j)}),...,F_X^{-1}(1-U_n^{(j)}))
Yj=g(FX−1(U1(j)),...,FX−1(Un(j)))−Yj′=−g(FX−1(1−U1(j)),...,FX−1(1−Un(j)))
都是非减函数。因此,
E
[
Y
(
−
Y
′
)
]
≥
E
[
Y
]
E
[
−
Y
′
]
⟶
E
[
Y
Y
′
]
≤
E
[
Y
]
E
[
Y
′
]
E[Y(-Y')]\ge E[Y]E[-Y']\longrightarrow E[YY']\le E[Y]E[Y']
E[Y(−Y′)]≥E[Y]E[−Y′]⟶E[YY′]≤E[Y]E[Y′]
所以,
C
o
v
(
Y
,
Y
′
)
=
E
[
Y
Y
′
]
−
E
[
Y
]
E
[
Y
′
]
<
0
Cov(Y,Y')=E[YY']-E[Y]E[Y']<0
Cov(Y,Y′)=E[YY′]−E[Y]E[Y′]<0
得证。
Antithetic variable approach 对偶变量法十分容易实现。只需要重复生成m/2组变量,剩下的一半可以自动得到。
θ
^
=
1
m
{
Y
1
+
Y
1
′
+
,
.
.
.
+
,
Y
m
/
2
+
Y
m
/
2
′
}
=
2
m
∑
j
=
1
m
/
2
(
Y
j
+
Y
j
′
2
)
\hat \theta =\frac{1}{m}\{Y_1+Y_1'+,...+,Y_{m/2}+Y_{m/2}' \}\\ =\frac{2}{m}\sum_{j=1}^{m/2}\bigg(\frac{Y_j+Y_j'}{2}\bigg)
θ^=m1{Y1+Y1′+,...+,Ym/2+Ym/2′}=m2j=1∑m/2(2Yj+Yj′)
也就是说,不仅需要生成的变量的数目变少了,方差还变小了。
MC.Phi <- function(x, R = 10000, antithetic = TRUE) {
u <- runif(R/2)
if (!antithetic) v <- runif(R/2) else
v <- 1 - u
u <- c(u, v)
cdf <- numeric(length(x))
for (i in 1:length(x)) {
g <- x[i] * exp(-(u * x[i])^2 / 2)
cdf[i] <- mean(g) / sqrt(2 * pi) + 0.5
}
cdf
}
对估计量的效果进行对比:
x <- seq(.1, 2.5, length=5)
Phi <- pnorm(x)
set.seed(123)
MC1 <- MC.Phi(x, anti = FALSE)
set.seed(123)
MC2 <- MC.Phi(x)
print(round(rbind(x, MC1, MC2, Phi), 5))
更准,就说明效率是更高了。
再对比使用antithetic方法和不使用得到的方差:
m <- 1000
MC1 <- MC2 <- numeric(m)
x <- 1.95
for (i in 1:m) {
MC1[i] <- MC.Phi(x, R = 1000, anti = FALSE)
MC2[i] <- MC.Phi(x, R = 1000)
}
方差减少了99.5%,效果拔群。