Vica's Blog
  2. Note

BMSI系列的第四篇，介绍贝叶斯框架下基本的抽样方法

手敲\(\LaTeX\)难免出现纰漏，有任何疑似错误或者不清楚的地方请直接在下方评论区留言，谢谢各位读者。

引入

我们假设读者已经对贝叶斯框架有了基本的认识，那么读者应该注意到，我们最关心的其实是一下两条

  • \(p(\theta\lvert y)\)，也就是我们感兴趣的参数的后验分布
  • \(p(\tilde{y}\lvert y)\)，也就是有了已有数据后，预测新的数据

然而这两者通常都不会像前面那样有容易求解的闭式解，就算是简单的层次化正态模型，其中的计算量也并不小。因此在实践中，我们通常需要用模拟的方式来代替这一过程。这是基于这样的一个思想：只要我能获得后验分布中的样本，我就可以给出和后验分布有关的推断的数值解。

蒙特卡洛方法

蒙特卡洛不是一个人名而是赌场名，概率论和赌博有着深厚的渊源（笑）

蒙特卡洛方法Monte Carlo，简单说就是利用计算机的随机数(伪随机数)来解决计算问题。最经典的例子就是撒豆子然后求解圆周率了。使用R几行就能写出这样一个程序

小例子

1
2
3
4
5
6
7
set.seed(42)
get.pi <- function(n) {
    x = runif(n)
    y = runif(n)
    return(sum(x^2+y^2<1) / n * 4)
}
get.pi(10000) # 3.15

我们无需知道任何关于圆的性质，除了\(x^2+y^2<1\)表示点在圆内。然后随机生成10000个点，统计落在园内的比例就可以刻画出圆周率的值了。

当然，我们这里更多的是利用其算积分或者期望，这是因为如果\(X\sim f(X)\)，我们有

\[\mathbb{E}(g(X))=\int g(x)f(x)\,\mathrm{d}x\approx \frac{1}{S}\sum_{i=1}^{S}g(x_i)\]

可以证明，样本均值是\(\mathbb{E}(g(X))\)的无偏估计，而中心极限定理central limit theorem, CLT则保证了其收敛速度为\(O(\frac{1}{\sqrt{S}})\)，如果\(x_i\)是服从\(f(x)\)的独立同分布的样本。

考虑一个简单的积分\(\int_{0}^{1}\exp\{-x\}\,\mathrm{d}x\)，学过微积分的人不难直接算出其为\(1-\exp\{-1\}\approx 0.632\)，我们使用模拟得到接近的结果

1
2
3
set.seed(42)
mean(exp(-runif(10000))) # 0.6331721
1-exp(-1) # 0.6321206

再考虑一个稍微复杂点的积分\(\int_{2}^{4}\exp\{-x\}\,\mathrm{d}x\)，其值为\(\exp\{-2\}-\exp\{-4\}\approx 0.117\)$，使用模拟的方法，首先

\[\int_a^b g(x)\cdot 1\,\mathrm{d}x=(b-a)\int_a^b g(x)\frac{1}{b-a}\,\mathrm{d}x\]

本例中我们从runif(2,4)即\(f(x)\)中抽取样本，然后乘回\(b-a=2\)即可

1
2
3
set.seed(42)
2*mean(exp(-runif(10000,2,4))) # 0.117543
exp(-2)-exp(-4) # 0.1170196

再考虑一个更复杂的例子\(\Phi(1)=\int_{-\infty}^1\frac{1}{2\pi}\exp\{-x^2/2\}\,\mathrm{d}x\)，其值为pnorm(1)约为\(0.841\)，我们使用两种抽样方式计算，请读者思考并对比其异同。(注意，我们无法使用\((-\infty,1)\)上的均匀分布，但由于\(\int_{-\infty}^0\frac{1}{2\pi}\exp\{-x^2/2\}\,\mathrm{d}x=0.5\)，故我们只需计算\(\int_{0}^1\frac{1}{2\pi}\exp\{-x^2/2\}\,\mathrm{d}x\))

1
2
3
4
set.seed(42)
sum(rnorm(10000)<1)/10000 # 0.8446
0.5+mean(exp(-runif(10000)^2/2)/sqrt(2*pi)) # 0.8414712
pnorm(1) # 0.8413447

经过这三个简单的例子，我希望读者对蒙特卡洛计算积分有了更深入的认识。

如何从给定分布抽样

均匀分布

线性同余方法linear congruential generator, (LCG) 使用一个相当简单的方式生成一系列的随机数，具体的原理和方法在这里不展开。但均匀分布几乎是我们后面所有抽样的基础——因为其抽样足够简单而有效。

格点法

方法

格点法通过模拟累积分布密度函数来抽样。假设我们有了概率密度函数\(p(x)\)，我们先做处理

  1. 在其支撑集上选择一系列等距的点\(x_1,\ldots,x_n\)
  2. 设\(p_0=0\),计算\(p_i=p(x_i)\)
  3. 归一化上面的概率，\(\tilde{p}_i=p_i/\sum_{j=0}^{n}p_j\)
  4. 模拟累积分布密度函数\(\tilde{P}_i=\sum_{j=0}^{i}\tilde{p_j}\approx \mathbb{P}(X\leq x_i)\)

抽样时

  1. 从\((0,1)\)均匀分布中抽样\(u\)
  2. 寻找\(i\)使得\(P_{i-1}< u \leq P_{i}\)，然后返回\(x_i\)即可

例子

从贝塔分布\(\mathrm{Beta}(2,3)\)采样

1
2
3
4
5
6
7
8
9
10
11
12
13
my.rbeta <- function(n, shape1, shape2) {
    X = 0:10000 / 10000
    p = dbeta(X, shape1, shape2)
    p = p / sum(p)
    P = cumsum(p)
    u = runif(n)
    idxs = findInterval(u, P)
    return(X[idxs])
}
set.seed(42)
s <- my.rbeta(10000, 2, 3)
hist(s, prob=T)
curve(dbeta(x,2,3), 0, 1, add=T)

优劣

优点

  • 方法简单，计算简单
  • 对任何概率密度都可以使用
  • 能轻松扩展到高维密度函数

缺点

  • 需要确定格点的密度
  • 对于连续型密度函数，该方法产生的样本只能有有限种取值
  • 事实上，其本质还是从一个离散的分布里抽，只是这个离散的分布对我们的目标分布近似得比较好

逆密度法

该方法基于以下原理，若\(F(x)\)是随机变量\(X\)的累计分布密度函数，我们定义其的逆为

\[F^{-1}(u)=\inf\{x\lvert F(x)\leq u\}\]

对于连续型随机变量，必有

\[\mathbb{P}(F(X)\leq u)=\mathbb{P}(X\leq F^{-1}(u))=F(F^{-1}(u))=u\]

也就是说\(F(x)\sim U(0,1)\)

抽样方法

  1. 抽取独立同分布样本\(u_1,\ldots,u_n\sim U(0,1)\)
  2. 计算\(x_i=F^{-1}(u_i)\)

注意，有时候也用\(S(x)=1-F(x)\)代替\(F(x)\)，因为如果\(u\sim U(0,1)\)，那么\((1-u)\sim U(0,1)\)，当然我们仅当\(S(x)\)形式更简单的时候这么做，比如下面的例子。

例子

从指数分布\(x\sim\mathrm{Exp}(\mu)\)采样，可以得到\(F(x;\,\mu)=1-\exp\{-x/\mu\}\)，\(S(x;\,\mu)=\exp\{-x/\mu\}\)，\(S^{-1}(u;\,\mu)=-\mu\log(u)\)，于是

1
2
3
4
5
6
7
8
my.rexp <- function(n, mu) {
    u = runif(n)
    return (-mu * log(u))
}
set.seed(42)
s <- my.rexp(10000,1)
hist(s, prob=T)
curve(dexp(x,1), 0, 10, add=T)

优劣

优点

  • 直接从目标分布中采样
  • 不需要考虑其他参数，比如格点法的格点密度

缺点

  • \(F^{-1}(u)\)不一定有闭式解
  • 计算量可能比较大

利用现有关系

抽样方法

如果目标分布是可以由某个已知分布的样本转化得来，那么直接从那个已知分布里采样就好。最后再做相应的变换即可。

例子

  • 若\(X\sim \mathcal{N}(\mu,\sigma^2)\)，则\(Y=\exp\{X\}\sim \mathrm{Lognormal}(\mu,\sigma^2)\)
  • 若\(X\sim \mathcal{N}(0, 1)\)，则\(Y=X^2\sim \chi_1^2\)
  • 若\(X_\alpha\sim \mathrm{Gamma}(1,\alpha),\,X_\beta\sim \mathrm{Gamma}(1,\beta)\)，则\(Y=x_\alpha/(X_\alpha+X_\beta)\sim \mathrm{Beta}(a,b)\)
  • 若\(X\sim U(0, 1)\)，则\(Y=-\log X\sim \mathrm{Exp}(1)\)

事实上，逆密度法就是本法的一个特例（为什么？）

1
2
3
4
5
6
7
8
9
my.rbeta2 <- function(n, shape1, shape2) {
    x1 = rgamma(n, shape1)
    x2 = rgamma(n, shape2)
    return (x1 / (x1 + x2))
}
set.seed(42)
s <- my.rbeta2(10000, 2, 3)
hist(s, prob=T)
curve(dbeta(x,2,3), 0, 1, add=T)

优劣

优点

  • 与从目标分布中采样完全等价

缺点

  • 不一定有好用的关系
  • 计算\(\log,\sin\)等可能代价也比较高

拒绝采样法

其实他的思想也十分简单，也十分“蒙特卡洛”。正如一开始我们使用随机的方式求圆周率一样，他使用随机的方式求积分的结果——寻找一个能够“盖住”原分布的“分布”（这里加引号是因为，如果一个分布（积分为1）要能盖住另一个分布，那么它两只能几乎处处相等——显然这样的分布没有什么用，因此我们只需要那个“盖子”的积分为有限值就行）然后我们从这样的分布中抽样，然后剔除掉高于目标积分的区域即可。用数学语言写作如下，对于一个目标分布\(f(x)\)，我们寻找一个分布\(g(x)\)，满足\(f(x)\leq c\cdot g(x)\overset{\Delta}{=}h(x),\,\forall x\)。显然这里通常\(c>1\)。

抽样方法

  1. 从\(g(x)\)中采样得样本\(x_0\)
  2. 计算该样本的接受率\(r(x_0)=\frac{f(x_0)}{cg(x_0)}=\frac{f(x_0)}{h(x_0)}\leq 1\)
  3. 从均匀分布中抽样\(u\sim U(0,1)\)
    1. 若\(u\leq r(x_0)\)，则接受该样本，返回样本\(x_0\)
    2. 否则，拒绝该样本
  4. 重复上述步骤直到采样出所需的样本容量

这里我们证明这个方法是正确的，记事件\(I\)为样本被接受。那么

\[\begin{align} \mathbb{P}(I=1)&=\int \mathbb{P}(I=1\lvert X=x)g(x)\,\mathrm{d}x\\ &=\int r(x)g(x)\,\mathrm{d}x\\ &=\int f(x)/c\,\mathrm{d}x=1/c \end{align}\]

于是，被接受的样本中

\[\mathbb{P}(x\lvert I=1)=p(x,I=1)/p(I=1)=\frac{f(x)}{c}\cdot c=f(x)\]

需要注意的是，\(1/c\)表征着抽样的效率，因为平均抽取\(c\)个样本后才接受1个。但另一方面，由于需要满足“盖住”的需求，故\(c\)不可能太小。显然最优的\(c=\sup f(x)/g(x)\)

例子

设\(f(x)=2\phi(x) I(x\geq 0)=\sqrt{2/\pi}\exp\{-0.5 x^2\}I(x\geq 0)\)，即半边正态分布。我们选择\(g(x)=\exp\{-x\}I(x\geq 0)\)。

首先我们计算最小的\(c\)值，不难看出\(x=1\)时\(f(x)/g(x)\)取最大值\(\sqrt{2/\pi}\exp\{0.5\}\approx 1.315\)，故采样的效率约为\(1/c\approx 0.76\)。而此时接受率\(r(x)=\exp\{-0.5(x-1)^2\}\)。采样代码如下

当然，事实上更简单的方式是从正态分布抽出样本然后取绝对值。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
my.half.norm <- function(n) {
  c = sqrt(2/pi) * exp(0.5)
  # make it a little bit larger to reduce the probabilities that regenerate
  # because I like vectorize operation
  N = ceiling(n / (c - 0.1))
  samples = rexp(N, 1)
  ratios = exp(-0.5*(samples-1)^2)
  us = runif(N)
  samples = samples[us < ratios]
  if (length(samples) >= n) {
    return (samples[1:n])
  } else {
    return(c(samples, my.half.norm(n - length(samples))))
  }
}
set.seed(42)
s <- my.half.norm(10000)
hist(s, prob=T)
curve(2 * dnorm(x), 0, 20, add=T)

采样分布的构造

许多时候选择一个合适的采样分布并不容易。但对于\(\log f(x)\)严格上凸Concave的情形，我们知道其切线必定在整个曲线上方，也就起到了盖住的作用。而取了对数后的直线，也就是原始的指数。也就是说此时我们可以用指数分布盖住它。

这里，我们考察一个例子\(\mathrm{Gamma}(2,1)\)，其概率密度函数\(f(x)=x\exp\{-x\}\)，而\(\log f(x)=\log(x)-x\),\(\frac{\mathrm{d}^2\log f(x)}{\mathrm{d}x^2}=-1/x^2<0\)，故满足严格上凸的条件。同时我们知道其概率密度的极大值出现在\(x=1\)处，因此我们将其分为左右两边，分别寻找采样分布。

出于简单考虑，我们假设两条切线的切点分别在\(x_l=0.5\)和\(x_r=2\)处，立刻可以写出两条切线

\[\begin{align} \mathrm{l:}&\quad y-\log(0.5)+0.5=1\cdot (x-0.5)\\ \mathrm{r:}&\quad y-\log(2)+2=-0.5\cdot(x-2)\\ \end{align}\]

做出图如下

1
2
3
curve((log(x)-x), 0, 6, ylim=c(-4,0))
abline(a=log(0.5)-1, b=1, col='red')
abline(a=log(2)-1, b=-0.5, col='red')

转化为原始形式，指数分布有\(\lambda\exp\{\lambda x\}\)的形式，则

\[\begin{align} h_l(x)&=c_lg_l(x)=0.5\cdot\exp\{(x-1)\}\mathrm{1}(x\leq 1)\\ h_r(x)&=c_rg_r(x)=4\exp\{-1.5\}\cdot0.5\exp\{-0.5(x-1)\}\mathrm{1}(x> 1) \end{align}\]

做出图如下

1
2
3
curve(dgamma(x,2,1), 0, 6, ylim=c(0,0.5))
curve(0.5 * exp(x-1), 0, 2, col='red', add=T)
curve(4 * exp(-1.5) * 0.5 * exp(-0.5*(x-1)), 0, 6, col='red', add=T)

此时我们的采样分布为

\[g(x)=\frac{c_l}{c_l+c_r}g_l(x)+\frac{c_r}{c_l+c_r}\]

因此我们可以先“掷硬币”——\(u\sim U(0,1)\)，然后根据是否小于\(c_l/(c_l+c_r)\)决定从哪个分布采样，剩下的事情就和之前一样了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
my.gamma21 <- function(n) {
  c_l = 0.5
  c_r = 4*exp(-1.5)
  N = 2 * n
  left.cnt  = sum(runif(N) < c_l / (c_l + c_r))
  left.samples = 1 - rexp(left.cnt, 1)
  right.samples = 1 + rexp(N - left.cnt, 0.5)
  left.ratios = dgamma(left.samples,2,1) / (0.5 * exp(left.samples-1))
  right.ratios = dgamma(right.samples,2,1) / (4 * exp(-1.5) * 0.5 * exp(-0.5*(right.samples-1)))
  left.samples = left.samples[runif(left.cnt) < left.ratios]
  right.samples = right.samples[runif(N - left.cnt) < right.ratios]
  samples = c(left.samples, right.samples)
  if (length(samples) >= n) {
    return (samples[1:n])
  } else {
    return(c(samples, my.gamma21(n - length(samples))))
  }
}
set.seed(42)
s <- my.gamma21(10000)
hist(s, prob=T, breaks=30)
curve(dgamma(x,2,1),add=T)

优劣

有一点必须要提前说：敏锐的读者可能已经发现，我们其实并不需要\(f(x)\)为一个分布，就算其积分不是1也丝毫没有关系。这有什么意义呢？还记得我们的后验分布\(p(\theta\lvert y)\propto p(y\lvert \theta)p(\theta)\)么，我们不需要计算这个归一化因子，直接从这个形式就能从后验分布里抽样了，减少了许多计算上的困扰。

优点

  • 样本服从目标分布
  • 很多目标分布都能使用这一方法
  • 采样分布也有很多选择的空间，相当一部分目标分布已经有比较成熟的采样分布的选择方式了

缺点

  • 当目标分布形式不好，或采样分布选取不当时可能导致采样效率低下，即\(c\)很大。在高维情况下尤为明显
  • 有些时候，采样分布的选择还是需要构思的

重要性抽样

读者在上一节的拒绝采样中可能会有这样的疑惑——我为什么一定要让接受率小于1呢？如果我让“接受率”可以大于1，相当于接受了大于1个样本；而我在小于1的时候也不通过从均匀分布抽样的方式判断是否接受，而相当于接受了零点几个样本，其本质上对于你做相关的计算(如积分)应该没有太大的影响。这样的想法对不对呢，我们用数学的语言阐述一下。

我们有目标分布\(f(x)\)和抽样分布\(g(x)\)，此时我们不要求其归一化。每当我们从\(g(x)\)中抽样得到\(x_i\)，我们计算其重要性\(w_i=f(x_i)/g(x_i)\)，并将元组\((x_i,\,w_i)\)叫做带权样本。

此时我们想计算\(h(x)\)在\(x\sim f(x)\)下的期望，则

\[\mathbb{E}_fh(X)=\int h(x)f(x)\,\mathrm{d}x=\int h(x)w(x)g(x)\,\mathrm{d}x=\mathbb{E}_g\left[h(x)w(x)\right]\]

也就是说，对于\(g(x)\)中独立抽出的一组样本\(x_1,\dots,x_n\)，我们使用\(\frac{1}{n}\sum_{i=1}^{n}h(x_i)w(x_i)\)就可以作为\(\mathbb{E}_fh(X)\)的一个无偏估计。

抽样方法

从\(g(x)\)中抽样，然后计算权重后成为带权样本，然后使用其进行计算。

例子

我们的目标分布是标准正态分布被截到区间\(\left[0,1\right]\)的一个分布，也就是\(f(x)=\phi(x)/\int_0^1\phi(x)\,\mathrm{d}x\ 0\leq x\leq 1\)。但我们现在感兴趣的是其均值。

当然，一个简单的做法就是我们从标准正态分布抽取样本，然后将不落在区间\(\left[0,1\right]\)的样本丢弃，但这样带来的弊端是，如果样本落在我们所感兴趣的区间的概率很小，那么有相当多的采样都是被浪费了；为了抽取1000个样本，我们可能进行了上万次的抽样，这是十分高昂的代价。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
my.bruteforce <- function(n) {
    remaining = n
    samples = numeric(0)
    while (remaining != 0) {
        tmps = rnorm(2 * remaining)
        tmps = tmps[(0 < tmps) & (tmps < 1)]
        if (length(tmps) > remaining) {
            samples = c(samples, tmps[1:remaining])
            remaining = 0
        } else {
            samples = c(samples, tmps)
            remaining = remaining - length(tmps)
        }
    }
    return(mean(samples))
}
set.seed(42)
my.bruteforce(10000) # 0.4581021

于是我们采用重要性抽样法，我们选取采样分布为\(U(0,1)\)，此时权重\(w(x)=f(x)/g(x)=\phi(x)/\int_0^1\phi(x)\,\mathrm{d}x / 1\)

1
2
3
4
5
6
7
my.imp <- function(n) {
    normalize = pnorm(1) - 0.5
    samples = runif(n)
    imps = (dnorm(samples) / normalize) / 1
    return(mean(samples * imps))
}
my.imp(10000) # 0.462356

我们再算算真值

1
2
integrate(function(x){x * dnorm(x)}, 0, 1)$value / integrate(dnorm, 0, 1)$value
# 0.4598622

从拒绝采样的算法中可以看到，没有样本会被拒绝，因此重要性检验从采出样本的效率上看是不能再好了。

然而，我们希望采出样本的权重的波动不那么大。直观的想法就是，如果不同样本的权重差别都很大，特别是如果有少量样本有着极高的权重的话，那我们在我们的一次实验中，如果我们的实验没有抽到这些样本可能就导致我们的结果有较大的偏差。因此我们希望我们的权重的波动，也就是方差越小越好。

最理想的情况当然是\(f(x)\)和\(g(x)\)成比例了，此时方差是0，大家的权重都是1。但这几乎是不太可行的，因此我们通常会选用一些看起来比较相似的分布——该有峰的地方大家都有峰，该趋于0的时候大家都趋于0。一个熟悉的例子就是正态分布和\(t\)分布了。

轻尾和重尾

轻尾和重尾都是相对而言的，简单说当趋向尾部（通常是正负无穷远）的时候，值较大的就称为重尾，反之称为轻尾。

1
2
curve(dnorm, -10, 10, col='blue')
curve(dt(x,1), -10, 10, col='green', add=T)

可以看到绿色的\(t\)分布相对于正态分布是重尾的。

在重要性抽样里，我们只能用更加重尾的分布作为轻尾分布的采样分布，而绝不能反过来。因为一旦采样分布更加轻尾，就意味着当\(x\)很大时，其样本的权重很大。此时，如果我们在某次实验中多抽中或少抽中了几个很大的\(x\)，就会导致计算的结果发生剧变。一般的，由于抽到极端值的概率相对较小，因此我们的估计绝大多数都会是少抽到了这些离群值，而有少数几次实验多抽到了几次离群值。而由于最后我们的估计数学上是无偏的，因此不难想象——我们的结果大多都是有小的偏差（比如比真值略小），而有某几次实验的结果远远高出真值。

Talk is cheap, just show the code.

使用\(t_3\)作为采样分布，求标准正态分布的均值（真值0）和方差（真值1）
1
2
3
4
5
6
7
8
9
10
11
12
13
# target distribution: normal distribution
# proposal distribution: t distribution
# it will work fine
set.seed(42)
f <- function (n) {
  samples <- rt(n, 3)
  imps <- dnorm(samples) / dt(samples, 3)
  m1 <- mean (samples * imps)
  m2 <- mean (samples * samples * imps)
  v <- m2 - m1 ^2
  print(paste('mu: ', toString(m1),', sigma : ', toString(v)))
}
f(1000) # "mu:  -0.0200670798678871 , sigma :  0.999501342932979"

可以看到结果很正常，而权重的分布中，有一些集中在0附近，也就是尾部的样本。但由于他们权重太小，翻不起浪，对结果也没太大影响。而与之平衡的，就是有部分样本权重大于1，也就是集中在0附近的样本。可以看到这些样本的最大权重也都不超过1.2（事实上，理论上权重的极大值也就大概为1.12）。因此得到的结果就比较稳定。

使用标准正态分布作为采样分布，求\(t_3\)的均值（真值0）和方差（真值3）
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# target distribution: t distribution
# proposal distribution: normal distribution
# it will be biased
set.seed(42)
g <- function (n) {
  samples = rnorm (n, 0, 1)
  imps = dt(samples ,3) / dnorm(samples)
  m1 = mean (samples * imps)
  m2 = mean (samples * samples * imps)
  v = m2 - m1 ^2
  print(paste('mu: ',toString(m1),', sigma : ', toString(v)))
  imps[imps > quantile (imps, .95)] = quantile (imps, .95)
  hist(imps, breaks = 25)
}
g(1000) # "mu:  0.00495172600802457 , sigma :  1.72064173002289"
g(10000) # "mu:  0.0538152854406483 , sigma :  1.92538642121928"
g(100000) # "mu:  0.00574719665677275 , sigma :  1.5573022759644"
g(10000000) # "mu:  0.00492975029866404 , sigma :  2.03582543863898"

可以看到我们给出的方差普遍偏低，随着样本量的加大虽然有所上升但也不太稳定，甚至还碰巧出现了大幅下跌。原因正是上面讲的，我们大多数实验都会有小幅偏差，而有少量实验会大幅偏差。直方图也说明了确有一部分样本有着极高的权重。

1
2
3
4
5
6
7
8
9
10
g2 <- Vectorize(function (n) {
  samples = rnorm (n, 0, 1)
  imps = dt(samples ,3) / dnorm(samples)
  m1 = mean (samples * imps)
  m2 = mean (samples * samples * imps)
  v = m2 - m1 ^2
  return(v)
})
set.seed(42)
plot(x=1:100, y=g2(5001:5100), type='l')

而通过重复跑多次实验，也确实能看到这样的情况：极个别跳到特别高，而绝大多数都是特别低的。

有效样本量

\[\begin{align} \mathrm{Var}(\widehat{h})&=\frac{1}{n}\mathrm{Var}_g\left[h_xg_x\right]\\ &=\frac{1}{n}\left[\mathbb{E}_g(h_x^2w_x^2)-(\mathbb{E}_gh_xw_x)^2\right]\\ \mathrm{Assume}\ h_x\perp w_x\rightarrow&=\frac{1}{n}\left[\mathbb{E}_g(h_x^2)\cdot \mathbb{E}_gw_x^2)-(\mathbb{E}_gh_x)^2(\mathbb{E}_gw_x)^2\right]\\ \mathrm{Notice that}\ \mathbb{E}_gw_x=\mathbb{E}_g\frac{f(x)}{g(x)}=1\rightarrow&=\frac{1}{n}\left[\mathbb{E}_g(h_x^2)\cdot \mathbb{E}_gw_x^2 -(\mathbb{E}_gh_x)^2\right]\\ \end{align}\]

\[n_{\mathrm{eff}}=n\cdot\frac{\mathrm{Var}(\tilde{h})}{\mathrm{Var}(\widehat{h})}=n\frac{\mathbb{E}_fh_x^2-(\mathbb{E}_fh_x)^2}{\mathbb{E}_g(h_x^2)\cdot \mathbb{E}_g(w_x)^2 -(\mathbb{E}_gh_x)^2}\approx \frac{n}{\mathbb{E}_gw_x^2}\]

其中最后一个约等号假定了\(\mathbb{E}_fh_x\approx \mathbb{E}_gh_x\approx 0,\,\mathbb{E}_fh_x^2\approx \mathbb{E}_gh_x^2\)

由此可以看到，如果分母\(\mathbb{E}_gw_x^2\)较大，则我们的样本效能会降低。

重要性重抽样

其抽样方法如下

  1. 从\(g(x)\)中采样出一系列样本\(\{x_1,\ldots,x_m\}\)，并计算对应的权重\(w(x_i)\)
  2. 从上述的采样集合中采样，每个样本被采样的概率正比于归一化后的权重\(\tilde{w}(x_i)=w(x_i)/\sum_{i=1}^{m}w(x_i)\)
  3. 如果抽得的样本已经被抽取过，则丢弃
  4. 无放回地重复步骤2-3直至抽取到足够的样本量

其一个代码实现如下——使用均匀分布作为先验，目标分布为\(\left[0,1\right]\)上的正态分布

1
2
3
4
5
6
7
8
9
10
11
12
my.rbeta2 <- function(n, shape1, shape2) {
    samples = runif(n)
    imps = dbeta(samples, shape1, shape2) / 1
    idxs = sample(1:n, size=n, replace=T, prob=imps)
    results = samples[idxs]
    return(results[1:n])
}
set.seed(42)
s <- my.rbeta2(10000, 2, 3)
mean(s) # 0.3984761
hist(s, prob=T)
curve(dbeta(x,2,3), 0, 1, add=T)
TAGS
 LATESTS