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

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

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

IS

估计量
构造q(x)q(x)q(x)去适应被积函数f(x)h(x)f(x)h(x)f(x)h(x)，解决上述的不匹配问题，其中q(x)q(x)q(x)满足，当f(x)h(x)≠0f(x)h(x)\ne0f(x)h(x)=0时，q(x)≠0q(x)\ne0q(x)=0
则此时积分的表达式可以写成I=∫abh(x)f(x)q(x)q(x)dx=Eq(h(x)w(x))I=\int_a^bh(x)\frac{f(x)}{q(x)}q(x)dx=E_q(h(x)w(x))I=∫abh(x)q(x)f(x)q(x)dx=Eq(h(x)w(x))其中w(x)=f(x)q(x)w(x)=\frac{f(x)}{q(x)}w(x)=q(x)f(x)
则估计量为InIS=1n∑i=1nh(xi)w(xi)I_n^{IS}=\frac 1n\sum_{i=1}^nh(x_i)w(x_i)InIS=n1i=1∑nh(xi)w(xi)其中xix_ixi独立地从分布q(x)q(x)q(x)抽样得到

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

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

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

SNIS（self-normalized IS）

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

估计量
可以把积分写成如下形式I=∫abh(x)f(x)q(x)q∗(x)dx∫abf(x)q(x)q∗(x)dxI=\frac{\int_a^bh(x)\frac{f(x)}{q(x)}q^*(x)dx}{\int_a^b\frac{f(x)}{q(x)}q^*(x)dx}I=∫abq(x)f(x)q∗(x)dx∫abh(x)q(x)f(x)q∗(x)dx其中q∗(x)q^*(x)q∗(x)为规范化后的函数，则估计量为InSNIS=∑i=1nh(xi)w(xi)∑i=1nw(xi)I_n^{SNIS}=\frac{\sum_{i=1}^nh(x_i)w(x_i)}{\sum_{i=1}^nw(x_i)}InSNIS=∑i=1nw(xi)∑i=1nh(xi)w(xi)其中xix_ixi独立地从分布q(x)q(x)q(x)抽样得到

期望和方差
（1）有偏，但依概率收敛到III（自证）
（2）方差的估计如下Var(InSNIS)≈σq,sn2n=Eq(w(x)2(h(x)−I)2)n=∫ab(f(x)h(x)−If(x))2q(x)dxnVar(I_n^{SNIS})\approx \frac{\sigma^2_{q,sn}}{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}Var(InSNIS)≈nσq,sn2=nEq(w(x)2(h(x)−I)2)=n∫abq(x)(f(x)h(x)−If(x))2dx
从该表达式可以看出 1）此统计量方差的表达式和IS的很相似，但不同
2）没有q(x)q(x)q(x)能让方差为0
3）有人证过q(x)q(x)q(x)的最优形式（让方差最小）为q(x)∝∣h(x)−I∣f(x)q(x)\propto |h(x)-I|f(x)q(x)∝∣h(x)−I∣f(x)，可以用此公式算出方差σq,sn2\sigma^2_{q,sn}σq,sn2理论上能达到的最小值为(Ef(∣h(x)−I∣))2(E_f(|h(x)-I|))^2(Ef(∣h(x)−I∣))2，这意味着SNIS能将减小原始MC方法的方差的比例为σq,sn2σ2≥(Ef(∣h(x)−I∣))2Ef(h(x)−I)2\frac{\sigma^2_{q,sn}}{\sigma^2}\geq \frac{(E_f(|h(x)-I|))^2}{E_f(h(x)-I)^2}σ2σq,sn2≥Ef(h(x)−I)2(Ef(∣h(x)−I∣))2

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

IS	拒绝抽样
优点：没有Mq(x)需要盖住f(x)的条件，能减小方差
缺点：需要记录w(x)	缺点：抽样效率低

q(x)常用的寻找办法

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

指数分布族
当f(x)f(x)f(x)属于指数分布族f(x)=g(x)eη(θ0)TT(x)−A(θ0)f(x)=g(x)e^{\eta(\theta_0)^TT(x)-A(\theta_0)}f(x)=g(x)eη(θ0)TT(x)−A(θ0)时，通常选择q(x)=g(x)eη(θ)TT(x)−A(θ)q(x)=g(x)e^{\eta(\theta)^TT(x)-A(\theta)}q(x)=g(x)eη(θ)TT(x)−A(θ)的形式，并且找出一个θ\thetaθ使得方差较小，此时IS估计量为InIS=1neA(θ)−A(θ0)∑i=1nh(xi)e(η(θ0)−η(θ))TT(xi)I_n^{IS}=\frac1ne^{A(\theta)-A(\theta_0)}\sum_{i=1}^nh(x_i)e^{(\eta(\theta_0)-\eta(\theta))^TT(x_i)}InIS=n1eA(θ)−A(θ0)i=1∑nh(xi)e(η(θ0)−η(θ))TT(xi)其中xix_ixi独立地从分布q(x)q(x)q(x)抽样得到

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

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

自适应IS

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

假设q(x)属于分布族q(x,θ\bm\thetaθ)，回忆IS估计量的方差的表达式，可得θ=arg⁡min⁡θ∫(h(x)f(x))2q(x,θ)dx\theta=\mathop{\arg\min}\limits_{\theta}\int\frac{(h(x)f(x))^2}{q(x,\theta)}dxθ=θargmin∫q(x,θ)(h(x)f(x))2dx机智的你可以通过上式写出θ\thetaθ的递推式θk+1=arg⁡min⁡θ1nk∑i=1nk(h(xi)f(xi))2q(xi,θk)2x∼q(x,θk)\bm{\theta_{k+1}=\mathop{\arg\min}\limits_{\theta}\frac1{n_k}\sum_{i=1}^{n_k}\frac{(h(x_i)f(x_i))^2}{q(x_i,\theta_k)^2}\qquad x\sim q(x,\theta_k)}θk+1=θargminnk1i=1∑nkq(xi,θk)2(h(xi)f(xi))2x∼q(x,θk)可以发现想要优化这个式子有点困难，所以我们的想法是不用方差来评价，换一种评价方法，首先回忆一下IS方法得到的中什么样的q(x)q(x)q(x)是好的，没错，就是q(x)∝t(x)q(x)\propto t(x)q(x)∝t(x)的时候，如果将t(x)t(x)t(x)标准化成密度函数t∗(x)t^*(x)t∗(x)的话，也就是说q(x)=t∗(x)q(x)=t^*(x)q(x)=t∗(x)的时候q(x)q(x)q(x)最好，还记得我们在第一节课基础知识中提到的KL距离吗，正是度量分布之间的距离的一种方法，恰好可以用于此处的优化问题。
于是我们的优化问题变成了θ=arg⁡min⁡θDKL(t∗∥qθ)=arg⁡min⁡θEt∗ln(t∗(x)qθ(x))=arg⁡max⁡θEt∗ln(qθ(x))\begin{aligned} \theta=\mathop{\arg\min}\limits_{\theta}D_{KL}(t^*\|q_\theta)=&\mathop{\arg\min}\limits_{\theta}E_{t^*}ln(\frac{t^*(x)}{q_\theta(x)})\\ = & \bm{\mathop{\arg\max}\limits_{\theta} E_{t^*}ln(q_\theta(x))} \end{aligned}θ=θargminDKL(t∗∥qθ)==θargminEt∗ln(qθ(x)t∗(x))θargmaxEt∗ln(qθ(x))机智的你再次可以通过上式写出θ\thetaθ的递推式θk+1=arg⁡max⁡θEθkln(q(x,θ))h(x)f(x)q(x,θk)\bm{\theta_{k+1}= \mathop{\arg\max}\limits_{\theta} E_{\theta_k}\frac{ln(q(x,\theta))h(x)f(x)}{q(x,\theta_k)}}θk+1=θargmaxEθkq(x,θk)ln(q(x,θ))h(x)f(x)假设q(x)q(x)q(x)属于指数分布族q(x)=g(x)eθTx−A(θ)q(x)=g(x)e^{\theta^Tx-A(\theta)}q(x)=g(x)eθTx−A(θ)写成离散形式即为θk+1=arg⁡max⁡θ1nk∑i=1nkh(xi)f(xi)q(xi,θk)ln(q(xi,θ))=arg⁡max⁡θ1nk∑i=1nkHiln(q(xi,θ))=arg⁡max⁡θ1nk∑i=1nkHi(θTxi−A(θ))\begin{aligned} \theta_{k+1}=& \mathop{\arg\max}\limits_{\theta}\frac1{n_k}\sum_{i=1}^{n_k}\frac{h(x_i)f(x_i)}{q(x_i,\theta_k)}ln(q(x_i,\theta))\\ =& \mathop{\arg\max}\limits_{\theta}\frac1{n_k}\sum_{i=1}^{n_k}H_i ln(q(x_i,\theta))\\ =&\mathop{\arg\max}\limits_{\theta}\frac1{n_k}\sum_{i=1}^{n_k}H_i(\theta^Tx_i-A(\theta)) \end{aligned}θk+1===θargmaxnk1i=1∑nkq(xi,θk)h(xi)f(xi)ln(q(xi,θ))θargmaxnk1i=1∑nkHiln(q(xi,θ))θargmaxnk1i=1∑nkHi(θTxi−A(θ))
则可以通过求解下式得到θk+1\theta_{k+1}θk+1 ∂∂θA(θ)=∑i=1nkHixi∑i=1nkHi\bm{\frac{\partial}{\partial \theta}A(\theta)=\frac{\sum_{i=1}^{n_k}H_ix_i}{\sum_{i=1}^{n_k}H_i}}∂θ∂A(θ)=∑i=1nkHi∑i=1nkHixi
特别的，当q(x)=N(θ,I)q(x)=N(\theta,I)q(x)=N(θ,I)时，θk+1=∑i=1nkHixi∑i=1nkHi\theta_{k+1}=\frac{\sum_{i=1}^{n_k}H_ix_i}{\sum_{i=1}^{n_k}H_i}θk+1=∑i=1nkHi∑i=1nkHixi
当q(x)=N(θ,Σ)q(x)=N(\theta,\Sigma)q(x)=N(θ,Σ)时，θk+1=Σ−1∑i=1nkHixi∑i=1nkHi\theta_{k+1}=\Sigma^{-1}\frac{\sum_{i=1}^{n_k}H_ix_i}{\sum_{i=1}^{n_k}H_i}θk+1=Σ−1∑i=1nkHi∑i=1nkHixi

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

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

Jn\bm{J_n}Jn的构造
通常来说，为了让方差尽可能小，我们希望Cov(InIS,Jn)Cov(I_n^{IS},J_n)Cov(InIS,Jn)尽可能大，所以直观上我们希望JnJ_nJn的表达式与InISI_n^{IS}InIS相似，这样其相关性会变高，再兼顾上易于计算JnJ_nJn的期望，我们构造出的JnJ_nJn的表达式如下Jn=wˉ=1n∑i=1nw(xi)x∼q(x,θk)J_n=\bar{w}=\frac1n\sum_{i=1}^nw(x_i)\qquad x\sim q(x,\theta_k)Jn=wˉ=n1i=1∑nw(xi)x∼q(x,θk)此时E(Jn)=1E(J_n)=1E(Jn)=1InISCV=InIS−λ(wˉ−1)I_n^{ISCV}=I_n^{IS}-\lambda(\bar{w}-1)InISCV=InIS−λ(wˉ−1)

Rao-Blackwellization条件期望法（RB）

原理
首先考虑一般情况，我们要估计I=E(h(X,Y))I=E(h(X,Y))I=E(h(X,Y)),其中(X,Y)(X,Y)(X,Y)的联合分布为f(x,y)f(x,y)f(x,y)，假设我们可以显式地算出E(h(X,Y)∣Y)\bm{E(h(X,Y)|Y)}E(h(X,Y)∣Y)，则由公式E(E(h(X,Y)∣Y))=E(h(X,Y))\bm{E(E(h(X,Y)|Y))=E(h(X,Y))}E(E(h(X,Y)∣Y))=E(h(X,Y))，我们可以构造出RB估计量InRB=1n∑i=1nE(h(xi,yi)∣yi)I_n^{RB}=\frac1n\sum_{i=1}^nE(h(x_i,y_i)|y_i)InRB=n1i=1∑nE(h(xi,yi)∣yi)则由公式Var(h(X,Y))=Var(E(h(X,Y)∣Y))+E(Var(h(X,Y)∣Y))\bm{Var(h(X,Y))=Var(E(h(X,Y)|Y))+E(Var(h(X,Y)|Y))}Var(h(X,Y))=Var(E(h(X,Y)∣Y))+E(Var(h(X,Y)∣Y))可得RB估计量的方差要小于原始蒙特卡罗估计的方差

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

原始的拒绝方法估计量为InMC=1n∑i=1Mh(yi)1{ui<w(yi)}I_n^{MC}=\frac1n\sum_{i=1}^Mh(y_i)\bm{1}_{\{u_i<w(y_i)\}}InMC=n1i=1∑Mh(yi)1{ui<w(yi)}RB估计量为InRB=1n∑i=1Mh(yi)v(yi)I_n^{RB}=\frac1n\sum_{i=1}^Mh(y_i)v(y_i)InRB=n1i=1∑Mh(yi)v(yi)其中v(yi)=EUi(1{Ui<w(yi)}∣yi)v(y_i)=E_{U_i}(\bm{1}_{\{U_i<w(y_i)\}}|y_i)v(yi)=EUi(1{Ui<w(yi)}∣yi)

统计计算第五节课，Mante Calor方法（二）——减小估计量的方差相关推荐

五节课从零起步(无需数学和Python基础)编码实现AI人工智能框架电子书V1
五节课从零起步 (无需数学和Python 基础) 编码实现AI 人工智能框架王家林 2018/4/15 ...
马建威android视频,5.25春季班高级班第三期第五节课课堂总结
来源:雪球App,作者: 阿里爸爸,(https://xueqiu.com/4583899269/180867769) 5.25春季班高级班第三期第五节课课堂总结今天继续由党斐大师给大家带来中炮先锋 ...
【java】兴唐第十五节课
知识点: 1.定义包名不能以java.为开头. 2.获取系统时间的方法: 注意: (1) 引入Date类时,引入的是java.util.Date 而不是java.sql.Date(后者是前者的子类). ...
王学岗——H265编码原理详解与码流分析(对应第五节课)
H264在分辨率高的情况下缺点 1,94年h264开始研制. 2,现在屏幕变的越来越大.假设h264以16X16进行编码,但我们的屏幕宽达到了32000像素.编码一帧画面需要400000个宏块.94年 ...
第五节课·基于图像相似度比较分镜头
目录基本思路 5.1 将视频打散为图片 5.2比较图像相似度 5.2.1 基于相等比较图像相似度 5.2.2 基于numpy计算图片是否相等 5.2.3 基于哈希(以均值哈希算法为例) 5.3 视频 ...
【java】兴唐第二十五节课小程序学生卡转账小系统（自己写的异常）
1.StuCard.java public class StuCard {public static void TransMoney(int source, int money, int target ...
【java】兴唐第二十五节课(异常和log4j的使用)
异常 1.try catch finally语法(附带多重catch) 代码实现: public static void main(String[] args) {try {int i = 1/0;} ...
学习笔记第二十五节课
正则介绍_grep 正则就是一串有规律的字符串,包含特殊符号. 对以后的写shell脚本很大的帮助.可以实现很多复杂的需求. 第一个工具 grep grep 用来过滤关键词. -c 行数,过滤出来的这 ...
第五节 CImage和CBmp（二）
由于这个库主要用于VC MFC下开发应用,因此目前可以下载到版本只支持MFC开发.如果您需要标准Win32的库,可以在回复时留下您的邮箱,或跟我联系. Email:wuchunlei@163.com ...

统计计算第五节课，Mante Calor方法（二）——减小估计量的方差

IS