统计计算 Monte Carlo积分和方差减少 学习笔记

Monte Carlo Integration and Variance Reduction

计算期望的数学方法：
$$E[g(X)]={\int }_{-\infty }^{\infty }g(x)f(x)dx$$

Simple Monte Carlo estimator 简单MC估计量

对于单位区间上的积分：
$$\hat{\theta }={\int }_0^{1}g(x)dx$$
可以有以下估计量。

$$\hat{\theta }=\overline{g_m(X)}=\frac{1}{m}g(X_i)$$
由强大数定律，\hat \theta会以1的概率收敛到
$$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]上的情况，可以用换元法。

$$\begin{gathered} y=\frac{t-a}{b-a} \\ t=(b-a)y+a \end{gathered}$$

$${\int }_a^{b}g(t)dt=(b-a){\int }_0^{1}g((b-a)y+a)dy$$

也可以用其它方式生成[a,b]区间上的随机数进行积分。比如
$${\int }_a^{b}g(t)dt=(b-a){\int }_a^{b}g(t)\frac{1}{b-a}dt$$
右式的积分便可以看成是对[a,b]上均匀分布的随机变量求MC积分，所以
$$\hat{\theta }=(b-a)\overline{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积分

$$\Phi (x)={\int }_{-\infty }^x\frac{1}{\sqrt{2\pi }}{e}^{-t^2/2}dt$$

对于无界的情况，不能直接套公式。但是，根据对称性，可以有
$$\begin{gathered} \Phi (x)={\int }_{-\infty }^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 \end{gathered}$$
又回到了有界区间积分的情况。如果x小于0，则可以考察
$$\Phi (x)=1-\Phi (-x)$$
再根据-x>0归入下面的情况。

对于如下积分，可以用[0,x]区间上的均匀分布随机数来做MC积分。
$${\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$$
但是，这样做每次都要生成一组x，不能复用。所以，可以考虑换元法.
$$y=\frac{t}{x}\in [0,1],t=xy$$

$${\int }_{-\infty }^x\frac{1}{\sqrt{2\pi }}{e}^{-t^2/2}dt={\int }_0^1\frac{1}{\sqrt{2\pi }}{e}^{-(xy{)}^2/2}xdy$$
现在，可以写出MC积分式(暂时略去常系数)
$$\begin{gathered} \theta =E_Y\left[x{e}^{-(xY{)}^2/2}\right] \\ \hat{\theta }=\overline{g_m(y)}=\frac{1}{m}\sum _{i=1}^{m}x{e}^{-(y_{i}x{)}^2/2} \end{gathered}$$

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\le x)]=P(Z\le x)=\Phi (x)$$

由此，可以用样本均值作为上式的估计量
$$\widehat{\Phi (x)}=\frac{1}{m}\sum _{i=1}^{m}I(z_i\le 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积分的总结

对于如下形式的积分，
$$\theta ={\int }_{A}g(x)f(x)dx，其中f(x)\ge 0,{\int }_{A}f(x)dx=1$$
可以先按照分布f(x)生成一组样本x1,x2,…xm，再求R.V.的函数的样本均值。
$$\hat{\theta }=\frac{1}{m}\sum _{i=1}^{m}g(x_i)$$
由强大数定律，它会收敛到\theta。

MC积分的标准差

对于
$${\sigma }^2=Va{r}_f(g(X))$$
我们可以写出
$$Var(\hat{\theta })=\frac{{\sigma }^2}{m}$$
然而，在估计的时候，我们并不能真正获得X的分布，因此只能靠样本来估计\sigma^2。因此，\hat \theta 的方差也只能是估计。
$$\frac{{\hat{\sigma }}^2}{m}=\frac{1}{m^2}\sum _{i=1}^m[g(x_i)-\overline{g(x)}{]}^2$$
右侧是用了plug-in估计的g(X)的方差。

由此，可以写出估计量的标准差
$$\widehat{se}(\hat{\theta })=\frac{\hat{\sigma }}{\sqrt{m}}=\frac{1}{m}\{\sum _{i=1}^m[g(x_i)-\overline{g(x)}{]}^2{\}}^{1/2}$$
这可以用在中心极限定理中。
$$\frac{\hat{\theta }-E[\hat{\theta }]}{\widehat{se}(\hat{\theta })}\sim N(0,1)asm\to \infty$$
接近于正态的\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 方差和效率

$$Var(\overline{g(X)})=\frac{1}{m}Var(g(X))$$

$$E(\overline{g(X)})=\frac{\theta }{b-a}$$

所以，
$$Var(\hat{\theta })=(b-a{)}^{2}Var(\overline{g(X)})=\frac{(b-a{)}^2}{m}Var(g(X))$$
由中心极限定理，当m很大时，\overline(g(X))近似于正态分布，所以\hat \theta也近似于正态分布，方差即由上式给出。

对于hit-or miss方法，有另外的方差表达式：
$$\widehat{F(x)}=\overline{g(X)}=\frac{1}{m}\sum _{i=1}^{m}I(X_i\le x)$$

$$\begin{gathered} Var(\widehat{F(x)})=p(1-p)/m=F(x)(1-F(x))/m \\ \widehat{se}(\widehat{F(x)})=\sqrt{\{\hat{F}(x)(1-\hat{F}(x))/m\}} \end{gathered}$$

Efficiency

一个估计量比另一个估计量更有效：
$$\frac{Var(\widehat{{\theta }_1})}{Var(\widehat{{\theta }_2})}<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$$
Y1,…,Ym就是独立同分布的函数，
$$E[\bar{Y}]=E\left[\frac{1}{m}\sum _{j=1}^m{Y}_j\right]=\theta$$
方差就有
$$Var(\hat{\theta })=Var\bar{Y}=\frac{Va{r}_{f}g(X)}{m}$$
增大采样数目可以减小MC估计量的方差。但是很不划算。一般来说，m的下界由下式给出
$$\begin{gathered} m\ge \lceil {\sigma }^2/e^2\rceil \\ m\ge \lfloor {\sigma }^2/e^2\rfloor \end{gathered}$$
其中sigma^2是目标随机变量的方差，e则是对于估计量我们所能接受的标准差。

Antithetic Variables Approach 对偶变量方法

考虑两个变量和的方差：
$$Var(\frac{U_1+U_2}{2})=\frac{1}{4}(Var(U_1)+Var(U_2)+2Cov(U_1,U_2))$$
当U1，U2独立时，方差式比较简单。当U1，U2 负相关时，方差可以减小。

This fact leads us to consider negtively correlated variables as a possible method for reducing variance.

考虑均匀分布：
$$U_j\sim Uinf(0,1)$$
则有
$$1-U_j\sim Unif(0,1)$$
且以上两组随机变量是负相关的。由此，我们考虑用逆变换：
$$\begin{gathered} Y_j=g(F_X^{-1}(U_1^{(j)}),...,F_X^{-1}(U_n^{(j)})) \\ {Y}_j^{\prime }=g(F_X^{-1}(1-U_1^{(j)}),...,F_X^{-1}(1-U_n^{(j)})) \end{gathered}$$

什么时候Y和Y’负相关呢？当g为单调函数的时候，以上两个随机变量负相关。

PROPOSITION：若X1,…Xn是独立的，f，g具有相同的单调性，那么
$$E[f(X)g(X)]\ge E[f(X)]E[g(X)]$$
证明：

因为f,g具有相同的单调性，所以，对任意x,y
$$\begin{gathered} (f(x)-f(y))(g(x)-g(y))\ge 0 \\ E[(f(X)-f(Y))(g(X)-g(Y))]\ge 0 \end{gathered}$$
因此，
$$E[f(X)g(X)]+E[f(Y)g(Y)]\ge 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)是单调函数，那么
$$\begin{gathered} Y_j=g(F_X^{-1}(U_1^{(j)}),...,F_X^{-1}(U_n^{(j)})) \\ {Y}_j^{\prime }=g(F_X^{-1}(1-U_1^{(j)}),...,F_X^{-1}(1-U_n^{(j)})) \end{gathered}$$
负相关。

证明：

不失一般性，认为g是单调非减函数，那么
$$\begin{gathered} Y_j=g(F_X^{-1}(U_1^{(j)}),...,F_X^{-1}(U_n^{(j)})) \\ -Y_j^{\prime }=-g(F_X^{-1}(1-U_1^{(j)}),...,F_X^{-1}(1-U_n^{(j)})) \end{gathered}$$
都是非减函数。因此，
$$E[Y(-Y^{\prime })]\ge E[Y]E[-Y^{\prime }]⟶E[Y{Y}^{\prime }]\le E[Y]E[Y^{\prime }]$$
所以，
$$Cov(Y,Y^{\prime })=E[Y{Y}^{\prime }]-E[Y]E[Y^{\prime }]<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′​}=m2​j=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%，效果拔群。