《Introducing Monte Carlo Methods with R》第四章笔记

误差理论，方差监测和加速收敛

引言——蒙特卡罗积分方法的基本误差理论

马尔科夫不等式

如果X仅取非负值，则对于任意a>0：$P(X\ge a)\le \frac{E(X)}{a}$

这个定理的证明过程在任何一本统计书上都有，就不详细描述。这个定理大概是说如果X全是正数的话，取值是个很大的数的可能性受一定的约束，约束条件和均值有关。
这个定理只是一个很粗略的估计，而且其哲学意义并不明显，但是以此为基础能推出一些非常有用的公式。

切比雪夫不等式

马尔科夫不等式的直接应用（将随机变量设为$\frac{(X-\mu {)}^2}{{\sigma }^2}$）就能推出如下定理：

设X均值为$\mu$，方差为${\sigma }^2$的随机变量，对于k>0，有：$P(|X-\mu |\ge k\sigma )\le \frac{1}{k^2}$

这个公式即为切比雪夫不等式，其哲学意义非常明显：一个随机变量X偏离其均值的可能性不会太大，具体来说，偏离k倍方差的概率至少要小于$\frac{1}{k^2}$。

大数定律

大数定律分为弱大数定律和强大数定律。
弱大数定律是切比雪夫公式的应用，叙述如下：

设$X_1,X_2,...,X_n$是均值为$\mu$的独立同分布随机变量，当$n\to 0$时，则对于任意$\epsilon >0$有：$P[(\frac{X_1+X_2+...+X_n}{n}-\mu )>\epsilon ]\to 0$

这个定理用语言来描述就是随着样本数n的无限增大，样本均值偏离实际均值的概率无限降低，趋近于0。

强大数定律是弱大数定律的推广：

$\underset{n\to \infty }{\lim }\frac{X_1+X_2+...+X_n}{n}=\mu$以概率1成立

根据强大数定律，可以确定一个随机变量的样本均值在样本量无穷大情况下收敛于理论均值。

中心极限定律

设$X_1,X_2,...,X_n$是均值为$\mu$且方差为${\sigma }^2({\sigma }^2<\infty )$的独立同分布随机变量，则对任意实数有：
$\underset{n\to \infty }{\lim }P(\frac{1}{\surd n}(X_1+X_2+...+X_n-n\mu )<x)=\Phi (x)$

这个定理，简写为CLT，是蒙特卡罗误差理论的基础。当蒙特卡罗生成的样本均可以被认为是从待估计的分布中抽出时，其样本均值偏离理论均值的误差可以用样本数和标准正态分布计算。从公式中可以看出，当n越大时误差越小。

估计的不确定性

这一节简单介绍几种常用的蒙特卡罗不确定度估计方法

直接法（蒙特卡罗的蒙特卡罗）

这种方法就是最简单直接的，反复运行同一种蒙特卡罗算法多次。只要运行次数足够多，就能完美无缺地进行不确定度的估计。

实验1 直接法进行蒙特卡罗数值积分的不确定度估计

待求的积分如下：
$$h(x)=[cos(50x)+sin(20x){]}^{2}x\in U(0,1)$$
R代码如下。结果图中淡黄色区域为蒙特卡罗之蒙特卡罗得出的（2.5%，97.5%）误差线，黑色线为第一次蒙特卡罗的结果，正好位于范围中。

> h=function(x){(cos(50*x)+sin(20*x))^2} #设置函数
> x=matrix(h(runif(200*10000)),ncol=200) #进行200次蒙特卡罗，每次10000个样本点
> estint=apply(x,2,cumsum)/(1:10000) #进行求均值
> y=apply(estint,1,quantile,c(0.025,0.975))  #求出每次蒙特卡罗的2.5%和97.5%的分位数
> polygon(c(1:10000,10000:1),c(y[1,],rev(y[2,])),col="whwat") #按照分位数，画出误差限
> lines(estint[,1],lwd=2) #画出第一次蒙特卡罗的结果，作为示例

实时监测——有效样本数（effective sample size）

当利用重要性采样（importance sampling）进行蒙特卡罗时，有效样本数（effective sample size，ESS）也可以作为方差的一种表征。ESS的计算方法如下：
$$\begin{gathered} {\omega }_{norm}^i={\omega }^i/\sum _{j=1}^n{\omega }^j \\ ES{S}_n=\frac{1}{{\sum }_{i=1}^n({\omega }_{norm}^i{)}^2} \end{gathered}$$
从公式中可以看出，$ESS\in [0,n]$，当所有样本都具有一样的权重时ESS=n，当全部样本都退化（权重为0）时，ESS=0。ESS越接近n说明蒙特卡罗效果越好。

除了ESS，还有另一种表征样本质量的方式：perplexity：
$$perplexity=-\sum _{i=1}^n{\omega }_{norm}^{i}log({\omega }_{norm}^i)$$
这里定义的perplexity可以被认为是提议分布g和目标分布f的K-L散度的一种近似，表征这两个分布的差别。很显然，差别越大说明重要性采样越成功。

实验2 贝叶斯推断并计算ESS

一个均值是柯西先验的正态分布：$X\sim N(\theta ,1)$，$\theta \sim Cauchy(0,1)$，当已有一个观测值x=2.5时，$\theta$的后验均值可以按下式计算：
$${\delta }^{\pi }(x)={\int }_{-\infty }^{\infty }\frac{\theta }{1+{\theta }^2}{e}^{-\frac{-(x-\theta {)}^2}{2}}d\theta /{\int }_{-\infty }^{\infty }\frac{1}{1+{\theta }^2}{e}^{-\frac{-(x-\theta {)}^2}{2}}d\theta$$
可以采用重要性采样近似以上均值：生成${\theta }_i\sim N(0,1)$，将$\frac{1}{1+{\theta }^2}$作为权重，则self-normalized形式的均值估计如下：
$${\delta }_n^{\pi }(x)=\sum _{i=1}^n\frac{{\theta }_i}{1+{\theta }_i^2}/\sum _{i=1}^n\frac{1}{1+{\theta }_i^2}$$
进行估计的R代码如下：

> norma=matrix(rnorm(500*10000),ncol = 500)+2.5 #生成500份N(0,1)样本，每份10000个
> weit=1/(1+norma^2)
> esti=apply(norma*weit,2,cumsum)/apply(weit,2,cumsum) #求估计均值
> polt(esti[,250],ylim=c(1.6,1.9))
> band=apply(esti,1,quantile,c(0.025,0.975)) #估计上下界
> polygon(c(1:10000,10000:1),c(band[1,],rev(band[2,]))) #画出上下界

计算ESS的代码如下，可以看出当对这个问题使用10000个样本进行重要性采样均值估计时，等效的有效样本数约等于5500。

> ess=apply(weit,2,cumsum)^2/apply(weit^2,2,cumsum)
> plot(ess[,500],ylim=c(0,10000))

实时监测——基于经验公式

在书的4.5节Simultaneous Monitoring中，介绍了利用经验公式的置信区间估计方法：
$$\begin{gathered} \theta (t)\in [\mu -u^{\ast }\sqrt{n{\widehat{{\sigma }_n}}^2}/round(nt),\mu +u^{\ast }\sqrt{n{\widehat{{\sigma }_n}}^2}/round(nt)] \\ {n}^{\ast }=a+b\sqrt{b} \end{gathered}$$
其中t为当前使用的样本数是总样本数的比例，即$t=\frac{n}{N}$，round为取整函数，a和b可以分别取0.1和3.15，$\widehat{{\sigma }_n}$为使用n个样本时的样本标准差。

自助法(bootstrap)

这是一种著名的估计方法，思想是通过重复抽样估计总体分布。详细来说就是将生成的样本的$F_n(x)$当做想要生成的$F(x)$的近似，$\hat{\theta }$是$\theta$的一个估计，通过从已生成样本重复有放回抽样的方法生成这个经验CDF的其他实现$F_n^{\ast }(x)$，就可以用新生成的$F_n^{\ast }(x)$的样本采用一样的方法得到${\hat{\theta }}^{\ast }$，利用一系列${\hat{\theta }}^{\ast }$实现$\hat{\theta }$的置信区间评定。

更进一步来说，自助法估计步骤如下：

1. 对第b次自助法重复试验（$b=1,...,B$）：
   （1）从观测样本$x_1,...,x_n$中有放回地抽样生成样本$x^{\ast (b)}=(x_1^{\ast },...,x_n^{\ast })$；
   （2）对第b个自助法样本计算估计值${\hat{\theta }}^{(b)}$；
2. $F_{\hat{\theta }}$分布的自助法估计的分布为${\hat{\theta }}^{(1)},...,{\hat{\theta }}^{(b)}$的经验分布。

对一个估计量$\hat{\theta }$的标准差进行自助法估计就是将自助法重复实验${\hat{\theta }}^{(1)},...,{\hat{\theta }}^{(B)}$的样本标准差作为估计值：
$$\begin{gathered} \widehat{se}({\hat{\theta }}^{\ast })=\sqrt{\frac{1}{B-1}\sum _{b=1}^B({\hat{\theta }}^{(b)}-\overline{{\hat{\theta }}^{(\ast )}})} \\ \overline{{\hat{\theta }}^{(\ast )}}=\frac{1}{B}\sum _{b=1}^B({\hat{\theta }}^{(b)}) \end{gathered}$$
B指的是运行自助法需要的实验次数，一般取50即可，很少会超过200。

当利用自助法对$\hat{\theta }$的偏差进行估计时：
$$\begin{gathered} bias(\hat{\theta })=\overline{{\hat{\theta }}^{(\ast )}}-\hat{\theta } \\ \overline{{\hat{\theta }}^{(\ast )}}=\frac{1}{B}\sum _{b=1}^B({\hat{\theta }}^{(b)}) \end{gathered}$$
其中${\hat{\theta }}^{(b)}$是第b次自助法抽样后对样本计算的估计值。

实验 自助法进行估计值的标准差和偏差估计

bootstrap包中含有两个数据集：law和law82，都是法学生的LSAT成绩和GPA数据。law数据是law82数据的一个较小的样本（15个），采用law的数据集对LSAT成绩和GPA的相关性进行计算，并用自助法实现相关系数标准差和偏差的估计，将结果和law82数据集的LSAT成绩和GPA的相关系数进行对比。

先进行相关系数的标准差估计

library(bootstrap)
B<-200 #自助法次数，抽样200次
n<-nrow(law) 
R<-numeric(B) 
for (b in 1:B)
{
  i<-sample(1:n,size=n,replace = TRUE) 
  #利用index进行有放回抽样，每次相当于重新实现一次law数据集
  LSAT<-law$LSAT[i] 
  GPA<-law$GPA[i]
  R[b]=cor(LSAT,GPA) #计算每次随机实现的law数据集中的LSAT成绩和GPA的相关系数
}

结果如下

> print(cor(law$LSAT,law$GPA)) #直接计算law数据集中的LAST成绩和GPA相关系数
[1] 0.7763745
> print(cor(law82$LSAT,law82$GPA)) #计算law82数据集中的LAST成绩和GPA相关系数，作为参考
[1] 0.7599979
> print(se.R<-sd(R)) #自助法估计的law数据集的LAST成绩和GPA相关系数的标准差
[1] 0.1219409
> hist(R,prob=TRUE) # 画R的分布直方图

之后进行相关系数的偏差估计

library(bootstrap)
theta.hat<-cor(law$LSAT,law$GPA) 

B<-20000
n<-nrow(law)
theta.b<-numeric(B)
for (b in 1:B)
{
  i<-sample(1:n,size=n,replace = TRUE)
  LSAT<-law$LSAT[i]
  GPA<-law$GPA[i]
  theta.b[b]=cor(LSAT,GPA) #计算每次随机实现的law数据集中的LSAT成绩和GPA的相关系数
  
}

bias<-mean(theta.b-theta.hat) 
#偏差即为200次law数据随机实现的估计值的均值和直接使用原始law数据集的估计量的差值

结果如下：

> print(bias)
[1] -0.005414397

计算出的bisa值可以说明law数据集对估计量的偏差值为-0.005，说明law数据集估计值倾向于高估。从上文可以看出，law数据集的相关系数估计值为0.7763745，law82数据集的相关系数估计值为0.7599979，和自助法估计的结果基本相符。

水手刀法（Jackknife）

水手刀法是比自助法更为古老的一种重抽样方法，也可以用于估计值的不确定估计。

水手刀法类似留一法，在第i次抽样中去除第i个样本，用剩下的$（i-1）$个样本计算${\hat{\theta }}_{(i)}$，如此得到${\hat{\theta }}_{(1)},...,{\hat{\theta }}_{(b)}$。设参数$\theta =t(F)$是分布F的函数，令$F_n$为一个服从分布F的随机样本的经验分布函数，$\theta$的嵌入式估计为$\hat{\theta }=t(F_n)$。如果数据的微小改变只会造成估计值的微小改变，$\hat{\theta }$的嵌入式估计是光滑的。

当利用水手刀法估计$\hat{\theta }$的估计偏差时：
$$\begin{gathered} {\widehat{bias}}_{jack-knife}=(n-1)(\overline{{\hat{\theta }}_{(\cdot )}}-\hat{\theta }) \\ \overline{{\hat{\theta }}_{(\cdot )}}=\frac{1}{n}\sum _{i=1}^n\widehat{{\theta }_i} \end{gathered}$$
其中$（n-1）$因子的原因和采用$\frac{1}{n}{\sum }_{i=1}^n(x_i-\overline{x}{)}^2$估计方差是有偏的类似，使用$(n-1)$因子修正偏差。

当使用水手刀法估计标准差时：
$${\widehat{se}}_{jack-knife}=\sqrt{\frac{n-1}{n}\sum _{i=1}^n(\widehat{{\theta }_{(i)}}-\overline{{\hat{\theta }}_{(\cdot )}})}$$
其中$\frac{n-1}{n}$因子的原因是因为当$\hat{\theta }=\overline{x}$时，$var(\overline{x})=\sqrt{var(x)/n}$，因此$\frac{n-1}{n}$因子使得${\widehat{se}}_{jack-knife}$成为无偏估计量。

必须要注意到，水手刀法进行估计只能对光滑的统计量使用，对于非光滑的统计量比如中位数就会失效。当统计量不光滑时候可以采用弃d-水手刀法（每次重复实验丢弃d个观测数据），当$\sqrt{n}/d\to 0$且$n-d\to 0$时，弃d-水手刀法对中位数的情况是可以准确估计的，但是n和d的上升会造成计算量迅速增大。

实验 水手刀法进行估计值的偏差和标准差估计

使用bootstrap包中的patch数据集进行水手刀法实验。patch数据集是一个医学数据集，包含使用一种医用贴片后病人血液中某种激素的测量结果。其中需要计算的统计量为：
$$\theta =\frac{E(新贴片)-E(旧贴片)}{E(旧贴片)-E(安慰剂)}$$
在数据集中，统计量可以表示为$\theta =\overline{Y}/\overline{Z}$，当$|\theta |<0.2$时认为新旧贴片效果相当。

library(bootstrap)
data("patch",package = "bootstrap")
n<-nrow(patch) #数据集的样本量即为水手刀法的采样次数
y<-patch$y
z<-patch$z
theta.hat<-mean(y)/mean(z) #使用公式直接计算估计量
print(theta.hat)
theta.jack<-numeric(n)
for (i in 1:n)
{
  theta.jack[i]<-mean(y[-i])/mean(z[-i]) 
  #计算每次留一法采样的估计值，[-i]表示除去数组中第i个元素
}
bias<-(n-1)*(mean(theta.jack)-theta.hat) #使用上文的水手刀法偏差公式计算估计量的偏差
print(bias)

结果如下：

> print(theta.hat)
[1] -0.0713061
> print(bias)
[1] 0.008002488

数值结果说明直接采用公式计算的$\theta =\overline{Y}/\overline{Z}$估计量有可能低估了0.008。

下面再采用水手刀法进行估计量的标准差估计，在完成偏差估计的基础上直接利用水手刀法标准差公式进行计算：

se<-sqrt((n-1)*mean((theta.jack-mean(theta.jack))^2))

结果如下：

> print(se)
[1] 0.1055278

自助法的水手刀法（Jackknife after bootstrap）

这种方法可以对自助法估计的偏差或者标准差在进行水手刀法方差估计，评判其准确程度。具体流程是：
（1）进行B次自助法抽样，每次抽样n个样本。计算$\theta$估计量和$\widehat{se}(\hat{\theta })$；
（2）对$i=1:n$，每次选取上一步所有bootstrap样本中不含有$x_n$的样本，并重新计算${\theta }_{j(i)}$；
（3）对上一步生成的所有${\theta }_{j(i)}$计算标准差，这个值即为$\widehat{se}(\hat{\theta })$的标准差。

实验 自助法的水手刀法（patch数据集）

data(patch,package='bootstrap')
n<-nrow(patch)
y<-patch$y
z<-patch$z
B<-2000
theta.b<-numeric(B)
indices<-matrix(0,nrow=B,ncol=n)

for (b in 1:B)
{
  i<-sample(1:n,size=n,replace = TRUE)
  y<-patch$y[i]
  z<-patch$z[i]
  theta.b[b]<-mean(y)/mean(z)
  indices[b, ]<-i
}

se.jack<-numeric(n)
for (i in 1:n)
{
  keep<-(1:B)[apply(indices,MARGIN = 1,FUN = function(k){!any(k==i)})]
  se.jack[i]<-sd(theta.b[keep])
}
print(sd(theta.b))
print(sqrt((n-1)*mean((se.jack-mean(se.jack))^2)))

方差缩减的方法

对偶变量法

当采用蒙特卡罗方法估计$\theta =E[X]$时，假设产生同分布的$X_1$，$X_2$，则：
$$\begin{array}{rl}Var(X)& =Var(\frac{X_1+X_2}{2})\\ & =\frac{1}{4}(Var(X_1)+Var(X_2)+2Cov(X_1,X_2))\end{array}$$
从式中可以看出，当$X_1$，$X_2$有负相关关系时，$Var(X)$可以减小。因此如果能够产生具有负相关关系的样本$X_1$，$X_2$，则估计量$E(\frac{X_1+X_2}{2})$相比估计量$E(X)$在生成同样多的样本数时方差更小。

因此，使用对偶变量法的关键是研究一种生成同分布但是负相关的样本$X_1$，$X_2$的算法。设$X_1\sim h(u_1,u_2,...,u_n)$，$u_1,...,u_n$为在[0,1]范围内的均匀分布随机变量，若$X_2\sim h(1-u_1,1-u_2,...,1-u_n)$则样本$X_1$，$X_2$服从相同的分布且明显负相关。采用这种方法生成蒙特卡罗仿真数据实现方差降低即为对偶变量法的原理。

控制变量法

当我们已知随机变量Y的期望$E(Y)={\mu }_y$时，待仿真的变量X的期望可以写成如下形式：
$$E(X)=E(X+c(Y-{\mu }_y))$$
其中c为待定参数。此时$X+c(Y-{\mu }_y)$的方差计算如下：
$$\begin{array}{rl}Var(X+c(Y-{\mu }_y))& =Var(X+cY)\\ & =Var(X)+c^{2}Var(Y)+2cCov(X,Y)\end{array}$$
此时若以方差最小为目标函数对参数c进行优化，可以得到：
$$\begin{array}{rl}{c}^{\ast }& =argmax(Var(X+c(Y-{\mu }_y)))\\ & =argmax(Var(X)+c^{2}Var(Y)+2cCov(X,Y))\\ & =-\frac{Cov(X,Y)}{Var(Y)}\end{array}$$
取$c=c^{\ast }$时，估计量$X+c(Y-{\mu }_y)$的最小方差计算如下：
$$Var(X+c(Y-{\mu }_y))=Var(X)-\frac{[Cov(X,Y){]}^2}{Var(Y)}$$
因此，当随机变量X和Y有负相关关系时，利用统计量$E(X+c(Y-{\mu }_y))$代替$E(X)$并取$c=c^{\ast }$可以降低$E(X)$的蒙特卡罗估计方差$100Corr(X,Y)\%$。随机变量Y即为待估计变量X的控制变量，$E(X+c(Y-{\mu }_y))$即为控制估计量。控制估计量的方差估计为：
$$Var(E(X+c(Y-{\mu }_y)))=\frac{1}{n}(Var(X)-\frac{Co{v}^2(X,Y)}{Var(Y)})$$

在实际应用中，通常$[Cov(X,Y)]$和$Var(Y)$都未知，也需要利用蒙特卡罗方法进行仿真估计，然后基于估计值计算$c^{\ast }$：
$$\widehat{c^{\ast }}=\frac{{\sum }_{i=1}^n(X_i-\overline{X})(Y_i-\overline{Y})}{{\sum }_{i=1}^n(Y_i-\overline{Y}{)}^2}$$

除了利用上面的公式直接计算$c^{\ast }$，还可以在获取X和Y的蒙特卡罗样本后建立$X=a+bY+e$的线性模型，其中e为随机误差。利用线性回归求得系数$\hat{a}$，$\hat{b}$，此时$\hat{b}=-\widehat{c^{\ast }}$，且$\overline{X}+\widehat{c^{\ast }}(\overline{Y}-{\mu }_y)=\hat{a}+\hat{b}{\mu }_y$，因此取$Y={\mu }_y$时的线性模型输出即为X估计值，估计方差可以认为是$\widehat{{\sigma }^2}/n$，$\hat{\sigma }$为回归方差。

条件期望法

根据条件方差公式
$$Var(X)=E(Var(X|Y))+Var(E(X|Y))$$
由于$E(Var(X|Y))>0$，明显$Var(E(X|Y))\le Var(X)$。由此而得$E(X|Y)$作为$E(X)$的估计量比直接计算仿真样本X的均值更加优秀，因为$E(X|Y)$不仅是$E(X)$的无偏估计，而且方差更小。因此，可以使用添加变量Y并仿真估计$E(X|Y)$的方法减小$E(X)$的估计方差。这一理论结果有时候也被称为“Rao-Blackwell Theorem”。

这种方法的模拟分为两阶段进行，先仿真Y样本，然后再根据Y的值生成X在Y条件下的样本，最后就能计算$E(X|Y)$。一种典型的实现方法是进行分层抽样，先对有k个取值的离散随机变量Y进行仿真，$p_i=P(Y=y_i),i=1,...,k$全部已知。之后对每个$y_i$，仿真$y_i$条件下的X样本，最后对所有的结果可以计算如下：
$$E(X)=\sum _{i=1}^{k}E(X|Y=y_i)p_i$$