这是我上的统计计算课讲的主要内容，写在这可以互相交流，有些地方我不是很理解的会标出来（用加粗斜体*标出），求大佬在留言处表达自己的看法，另外如果有啥问题也可以在留言处留言，如果我看到了会回复

这次的内容较为分散，知识点较多，这次的文章结构采用顺序写法，知识点之间的结构是以一个知识点为中心，向不同方向扩展，但达成的目的相同

背景：回忆MC方法要解决的问题是计算积分$I={\int }_a^{b}f(x)h(x)dx$，用到的估计量为$$\widehat{I_n}=\frac{1}{n}\sum _{i=1}^{n}h(x_i)$$其中$x_i$独立地从分布$f(x)$抽样得到，一个潜在的问题是：如果f(x)和h(x)不匹配，即如果积分是由f(x)值很小的一片区域A决定，那么可能会出问题，因为我们可能无法从区域A中抽样出足够多的样本点，会使方差变大。 解决的办法也很简单：就是从A中抽出足够多的样本点，这就产生了importance sampling（IS） 方法。

IS

估计量
构造$q(x)$去适应被积函数$f(x)h(x)$，解决上述的不匹配问题，其中$q(x)$满足，当$f(x)h(x)\ne 0$时，$q(x)\ne 0$
则此时积分的表达式可以写成$$I={\int }_a^{b}h(x)\frac{f(x)}{q(x)}q(x)dx=E_q(h(x)w(x))$$其中$w(x)=\frac{f(x)}{q(x)}$
则估计量为$$I_n^{IS}=\frac{1}{n}\sum _{i=1}^{n}h(x_i)w(x_i)$$其中$x_i$独立地从分布$q(x)$抽样得到

估计量的均值和方差
记$t(x)=h(x)w(x)$,则容易算出$$E(I_n^{IS})=I$$$$\begin{array}{rl}Var(I_n^{IS})=& \frac{1}{n}\int \frac{(f(x)h(x){)}^2}{q(x)}dx-I^2\\ =& \frac{1}{n}\int \frac{(f(x)h(x)-Iq(x){)}^2}{q(x)}dx\end{array}$$从表达式可以看出以下事情：（1）该估计是无偏的
（2）当$q(x)\propto f(x)h(x)$时，方差最小

回忆（一）中的收敛速率$$P(|\widehat{I_n}-I|<\frac{\sigma }{\sqrt{n}\delta })>1-\delta$$上面的IS估计量就是让$\sigma$减小，以加快收敛

特别的，当f(x)，h(x)都是正态分布时，**的q(x)也是正态分布，并且均值在f和h的均值之间（简单计算即得）

SNIS（self-normalized IS）

背景：当f或q不为规范时（即积分不为1）的情况

估计量
可以把积分写成如下形式$$I=\frac{{\int }_a^{b}h(x)\frac{f(x)}{q(x)}{q}^{\ast }(x)dx}{{\int }_a^b\frac{f(x)}{q(x)}{q}^{\ast }(x)dx}$$其中$q^{\ast }(x)$为规范化后的函数，则估计量为$$I_n^{SNIS}=\frac{{\sum }_{i=1}^{n}h(x_i)w(x_i)}{{\sum }_{i=1}^{n}w(x_i)}$$其中$x_i$独立地从分布$q(x)$抽样得到

期望和方差
（1）有偏，但依概率收敛到$I$（自证）
（2）方差的估计如下$$\begin{gathered} Var(I_n^{SNIS})\approx \frac{{\sigma }_{q,sn}^2}{n}=\frac{E_q(w(x{)}^2(h(x)-I{)}^2)}{n} \\ =\frac{{\int }_a^b\frac{(f(x)h(x)-If(x){)}^2}{q(x)}dx}{n} \end{gathered}$$
从该表达式可以看出 1）此统计量方差的表达式和IS的很相似，但不同
2）没有$q(x)$能让方差为0
3）有人证过$q(x)$的最优形式（让方差最小）为$q(x)\propto |h(x)-I|f(x)$，可以用此公式算出方差${\sigma }_{q,sn}^2$理论上能达到的最小值为$(E_f(|h(x)-I|){)}^2$，这意味着SNIS能将减小原始MC方法的方差的比例为$$\frac{{\sigma }_{q,sn}^2}{{\sigma }^2}\ge \frac{(E_f(|h(x)-I|){)}^2}{E_f(h(x)-I{)}^2}$$

IS和拒绝抽样方法的比较（待完善）

q(x)常用的寻找办法

上述只是理论上找出能让方差最小的$q(x)$的形式，但如果找出的$q(x)$的形式如果较为复杂，则不便于从$q(x)$中抽样$x$，所以下面介绍一些实际中常用$q(x)$的寻找方法

讯享网

指数分布族
当$f(x)$属于指数分布族$f(x)=g(x)e^{\eta ({\theta }_0{)}^{T}T(x)-A({\theta }_0)}$时，通常选择$q(x)=g(x)e^{\eta (\theta {)}^{T}T(x)-A(\theta )}$的形式，并且找出一个$\theta$使得方差较小，此时IS估计量为$$I_n^{IS}=\frac{1}{n}{e}^{A(\theta )-A({\theta }_0)}\sum _{i=1}^{n}h(x_i)e^{(\eta ({\theta }_0)-\eta (\theta ){)}^{T}T(x_i)}$$其中$x_i$独立地从分布$q(x)$抽样得到

海森矩阵与正态分布方法
假设我们已经找到了$t(x)$(之前定义过)的极大值点$x^{\ast }$，则由泰勒展开可以得到$$\begin{gathered} ln(t(x))\approx ln(t(x^{\ast }))-\frac{1}{2}(x-x^{\ast }{)}^T(-H^{\ast })(x-x^{\ast }) \\ 则有t(x)\approx t(x^{\ast })e^{-\frac{1}{2}(x-x^{\ast }{)}^T(-H^{\ast })(x-x^{\ast })} \end{gathered}$$所以可以取$q(x)=N(x^{\ast },-H^{\ast })$，其中$H^{\ast }$为$ln(t(x))$在$x^{\ast }$处的海森矩阵

混合正态方法
由于被积函数可能出现一些多峰的情况，所以为了让样本的密度函数与被积函数的形状更匹配，可以用混合正态作为$q(x)$，此时从$q(x)$中抽样有两种抽样方法，第一种是直接抽，第二种是先抽角标，再从对应角标的正态中抽样，第二种方法的优点是抽样速度更快，缺点是方差会比第一种更大

自适应IS

通过上面的讨论，我们发现找一个让估计量方差较小并且便于抽样的$q(x)$是一件困难的事情，尤其是对于刚才讨论的“指数分布族”和“混合正态方法”，我们已经确定了$q(x)$的形式，但还未确定参数，如何找一个较好的参数似乎是首要问题，本节将介绍寻找参数的方法（本质上仍然是优化问题，所以很自然地我们想到应该迭代地去找）

假设q(x)属于分布族q(x,$\mathbf{\theta }$)，回忆IS估计量的方差的表达式，可得$$\theta =\underset{\theta }{\mathrm{argmin}}\int \frac{(h(x)f(x){)}^2}{q(x,\theta )}dx$$机智的你可以通过上式写出$\theta$的递推式$${\mathbf{\theta }}_{\mathbf{k}+1}=\underset{\mathbf{\theta }}{\mathrm{argmin}}\frac{1}{{\mathbf{n}}_{\mathbf{k}}}\sum _{\mathbf{i}=1}^{{\mathbf{n}}_{\mathbf{k}}}\frac{(\mathbf{h}({\mathbf{x}}_{\mathbf{i}})\mathbf{f}({\mathbf{x}}_{\mathbf{i}}){)}^2}{\mathbf{q}({\mathbf{x}}_{\mathbf{i}},{\mathbf{\theta }}_{\mathbf{k}}{)}^2}\mathbf{x}\sim \mathbf{q}(\mathbf{x},{\mathbf{\theta }}_{\mathbf{k}})$$可以发现想要优化这个式子有点困难，所以我们的想法是不用方差来评价，换一种评价方法，首先回忆一下IS方法得到的中什么样的$q(x)$是好的，没错，就是$q(x)\propto t(x)$的时候，如果将$t(x)$标准化成密度函数$t^{\ast }(x)$的话，也就是说$q(x)=t^{\ast }(x)$的时候$q(x)$最好，还记得我们在第一节课基础知识中提到的KL距离吗，正是度量分布之间的距离的一种方法，恰好可以用于此处的优化问题。
于是我们的优化问题变成了$$\begin{array}{rl}\theta =\underset{\theta }{\mathrm{argmin}}{D}_{KL}(t^{\ast }\Vert q_{\theta })=& \underset{\theta }{\mathrm{argmin}}{E}_{t^{\ast }}ln(\frac{t^{\ast }(x)}{q_{\theta }(x)})\\ =& \underset{\mathbf{\theta }}{\mathrm{argmax}}{\mathbf{E}}_{{\mathbf{t}}^{\ast }}\mathbf{l}\mathbf{n}({\mathbf{q}}_{\mathbf{\theta }}(\mathbf{x}))\end{array}$$机智的你再次可以通过上式写出$\theta$的递推式$${\mathbf{\theta }}_{\mathbf{k}+1}=\underset{\mathbf{\theta }}{\mathrm{argmax}}{\mathbf{E}}_{{\mathbf{\theta }}_{\mathbf{k}}}\frac{\mathbf{l}\mathbf{n}(\mathbf{q}(\mathbf{x},\mathbf{\theta }))\mathbf{h}(\mathbf{x})\mathbf{f}(\mathbf{x})}{\mathbf{q}(\mathbf{x},{\mathbf{\theta }}_{\mathbf{k}})}$$假设$q(x)$属于指数分布族$q(x)=g(x)e^{{\theta }^{T}x-A(\theta )}$写成离散形式即为$$\begin{array}{rl}{\theta }_{k+1}=& \underset{\theta }{\mathrm{argmax}}\frac{1}{n_k}\sum _{i=1}^{n_k}\frac{h(x_i)f(x_i)}{q(x_i,{\theta }_k)}ln(q(x_i,\theta ))\\ =& \underset{\theta }{\mathrm{argmax}}\frac{1}{n_k}\sum _{i=1}^{n_k}{H}_{i}ln(q(x_i,\theta ))\\ =& \underset{\theta }{\mathrm{argmax}}\frac{1}{n_k}\sum _{i=1}^{n_k}{H}_i({\theta }^T{x}_i-A(\theta ))\end{array}$$
则可以通过求解下式得到${\theta }_{k+1}$ $$\frac{\partial }{\partial \mathbf{\theta }}\mathbf{A}(\mathbf{\theta })=\frac{{\sum }_{\mathbf{i}=1}^{{\mathbf{n}}_{\mathbf{k}}}{\mathbf{H}}_{\mathbf{i}}{\mathbf{x}}_{\mathbf{i}}}{{\sum }_{\mathbf{i}=1}^{{\mathbf{n}}_{\mathbf{k}}}{\mathbf{H}}_{\mathbf{i}}}$$
特别的，当$q(x)=N(\theta ,I)$时，${\theta }_{k+1}=\frac{{\sum }_{i=1}^{n_k}{H}_i{x}_i}{{\sum }_{i=1}^{n_k}{H}_i}$
当$q(x)=N(\theta ,\Sigma )$时，${\theta }_{k+1}={\Sigma }^{-1}\frac{{\sum }_{i=1}^{n_k}{H}_i{x}_i}{{\sum }_{i=1}^{n_k}{H}_i}$

控制变量法(Control Variates)(CV)

原理
构造新的估计量$$I_n^{ISCV}=I_n^{IS}-\lambda (J_n-J)$$ 其中$E(J_n)=J$，则$I_n^{ISCV}$显然是$I$的一个无偏估计，其方差为$$Var(I_n^{ISCV})=Var(I_n^{IS})-2\lambda Cov(I_n^{IS},J_n)+{\lambda }^{2}Var(J_n)$$取$\hat{\lambda }=\frac{Cov(I_n^{IS},J_n)}{Var(J_n)}$可得估计量$I_n^{ISCV}$的方差的最小值，易算出此最小方差比原始方差$Var(I_n^{IS})$小，实际操作中$\hat{\lambda }$由对应分布的样本估计

${\mathbf{J}}_{\mathbf{n}}$的构造
通常来说，为了让方差尽可能小，我们希望$Cov(I_n^{IS},J_n)$尽可能大，所以直观上我们希望$J_n$的表达式与$I_n^{IS}$相似，这样其相关性会变高，再兼顾上易于计算$J_n$的期望，我们构造出的$J_n$的表达式如下$$J_n=\bar{w}=\frac{1}{n}\sum _{i=1}^{n}w(x_i)x\sim q(x,{\theta }_k)$$此时$E(J_n)=1$$$I_n^{ISCV}=I_n^{IS}-\lambda (\bar{w}-1)$$

Rao-Blackwellization条件期望法（RB）

原理
首先考虑一般情况，我们要估计$I=E(h(X,Y))$,其中$(X,Y)$的联合分布为$f(x,y)$，假设我们可以显式地算出$\mathbf{E}(\mathbf{h}(\mathbf{X},\mathbf{Y})|\mathbf{Y})$，则由公式$\mathbf{E}(\mathbf{E}(\mathbf{h}(\mathbf{X},\mathbf{Y})|\mathbf{Y}))=\mathbf{E}(\mathbf{h}(\mathbf{X},\mathbf{Y}))$，我们可以构造出RB估计量$$I_n^{RB}=\frac{1}{n}\sum _{i=1}^{n}E(h(x_i,y_i)|y_i)$$则由公式$\mathbf{V}\mathbf{a}\mathbf{r}(\mathbf{h}(\mathbf{X},\mathbf{Y}))=\mathbf{V}\mathbf{a}\mathbf{r}(\mathbf{E}(\mathbf{h}(\mathbf{X},\mathbf{Y})|\mathbf{Y}))+\mathbf{E}(\mathbf{V}\mathbf{a}\mathbf{r}(\mathbf{h}(\mathbf{X},\mathbf{Y})|\mathbf{Y}))$可得RB估计量的方差要小于原始蒙特卡罗估计的方差

与拒绝抽样方法结合构造出具体的算法
假设拒绝抽样在随机时间M(随机时间的定义此处不给出，从名字直观理解就行）时停止抽样，此时得到n个样本$x_1,......,x_n$，记进行拒绝操作前的原始样本为$y_1,......,y_n$

原始的拒绝方法估计量为 I n M C = 1 n ∑ i = 1 M h ( y i ) 1 { u i < w ( y i ) } I_n^{MC}=\frac1n\sum_{i=1}^Mh(y_i)\bm{1}_{\{u_i<w(y_i)\}} InMC​=n1​i=1∑M​h(yi​)1{ ui​<w(yi​)}​RB估计量为$$I_n^{RB}=\frac{1}{n}\sum _{i=1}^{M}h(y_i)v(y_i)$$其中 v ( y i ) = E U i ( 1 { U i < w ( y i ) } ∣ y i ) v(y_i)=E_{U_i}(\bm{1}_{\{U_i<w(y_i)\}}|y_i) v(yi​)=EUi​​(1{ Ui​<w(yi​)}​∣yi​)