Chapter 6 MCMC方法

6.1 介绍

前面介绍蒙特卡洛方法估计积分 $$\int_Ag(t)\text{d}t,$$ 是把此积分表示成一个对某个概率密度 $f(t)$ 下的期望，从而积分估计问题转化成从目标概率密度 $f(t)$ 中产生随机样本。

在MCMC方法中，首先建立一个Markov链，使得 $f(t)$ 为其平稳分布，则可以运行此Markov链充分长时间直至收敛到平稳分布，那么从目标分布 $f(t)$ 中产生随机样本，就是从达到平稳状态的Markov链中产生其样本路径。

一个好的链应该具有快速混合性质——从任意位置出发很快达到平稳分布。

6.1.1 贝叶斯推断中的积分问题

从贝叶斯的观点来看，模型中的观测变量和参数都是随机变量。因此，样本 $\pmb x=(x_1,\dots,x_n)$ 和参数 $\theta$ 的联合分布可以表示为 $$f_{\pmb x,\theta}(\pmb x,\theta)=f_{\pmb x|\theta}(x_1,\dots,x_n)\pi(\theta).$$ 从而根据Bayes定理，可以通过样本 $\pmb x$ 的信息对 $\theta$ 的分布进行更新： $$f_{\theta|\pmb x}(\theta\mid\pmb x)=\dfrac{f_{\pmb x|\theta}(\pmb x)\pi(\theta)}{\int_\Theta f_{\pmb x|\theta}(\pmb x)\pi(\theta)\text{d}\theta}.$$ 则在后验分布下，$g(\theta)$ 的期望为 $$\mathbb E[g(\theta\mid\pmb x)]=\int_\Theta g(\theta)f_{\theta|\pmb x}(\theta\mid\pmb x)\text{d}\theta=\dfrac{\int_\Theta g(\theta)f_{\pmb x|\theta}(\pmb x)\pi(\theta)\text{d}\theta}{\int_\Theta f_{\pmb x|\theta}(\pmb x)\pi(\theta)\text{d}\theta}.$$ 在许多情况下，得到直接积分的结果是十分困难的，需要通过数值计算。

6.1.2 MCMC积分

积分 $\displaystyle\mathbb E[g(\theta\mid\pmb x)]=\int_\Theta g(\theta)f_{\theta|\pmb x}(\theta\mid\pmb x)\text{d}\theta$ 的蒙特卡洛估计为样本均值 $$\overline{g}=\dfrac{1}{m}\sum_{i=1}^m g(\theta_i),$$ 其中 $\theta_i$ 是从分布 $f_{\theta|\pmb x}(\theta\mid\pmb x)$ 中抽取的样本。当 $\theta_1,\dots,\theta_m$ 独立时，则可以根据大数律知当样本量 $m$ 趋于无穷时，$\overline{g}$ 收敛到 $\mathbb E[g(\theta\mid\pmb x)]$。

MCMC方法的理论依据在于下述的极限定理：

定理1. 设 $\{X_t,t\geq0\}$ 为一不可约、非周期的状态空间为 $\Omega$ 的Markov链，$\pi$ 为平稳分布，$\pi_0$ 为初始分布，则 $$\pi_t\rightarrow\pi,\quad t\rightarrow\infty.$$

这一定理揭示了当Markov链运行充分长时间后，链的状态分布即平稳分布，且与初始分布无关。

定理2. 若 $X_1,X_2,\dots$ 为一遍历的平稳分布为 $\pi$ 的Markov链值，则 $X_n$ 依分布收敛到分布为 $\pi$ 的随机变量 $X$，且对任意函数 $f$，当 $\mathbb E_\pi|f(X)|$ 存在时，且 $n\rightarrow\infty$ 时，则 $$\overline{f}_n=\dfrac{1}{n}\sum_{t=1}^nf(X_t)\rightarrow\mathbb E_\pi f(X),\quad\text{a.s.}$$

下面计算遍历均值 $\overline{f}_n$ 的方差。记 $f^t=f(X_t)$，$\gamma_k=\text{Cov}(f^t,f^{t+k})$，则 $f^t$ 的方差为 $\sigma^2=\gamma_0$，$k$ 步相关系数 $\rho_k=\gamma_k/\sigma^2$。从而 $$\begin{align*} \tau_n^2&=n\text{Var}_\pi(\overline{f}_n)=n\text{Var}_\pi\left(\dfrac{1}{n}\sum_{t=1}^nf^t\right)\\ &=\dfrac{1}{n}\sum_{t=1}^n\text{Var}_\pi(f^t)+\dfrac{2}{n}\sum_{k=1}^{n-t}\text{Cov}(f^t,f^{t+k})=\sigma^2\left[1+2\sum_{k=1}^{n-1}\dfrac{n-k}{n}\rho_k\right]. \end{align*}$$

可以证明 $$\tau_n^2\rightarrow\tau^2=\sigma^2\left[1+2\sum_{k=1}^{\infty}\rho_k\right].$$ 因此遍历均值 $\overline{f}_n$ 的方差为 $$\text{Var}_\pi(\overline{f}_n)=\dfrac{\sigma^2}{n}\left[1+2\sum_{k=1}^{n-1}\dfrac{n-k}{n}\rho_k\right].$$

定理3. 若一个链是一致的几何遍历（转移概率以速度 $\lambda^t,0<\lambda<1$ 收敛）的，以及 $f$ 相对于平稳分布 $\pi$ 是平方可积的，则有 $$\sqrt{n}\cdot\dfrac{\overline{f}_n-\mathbb E_\pi f(X)}{\tau}\rightarrow N(0,1),\quad n\rightarrow\infty.$$

这个定理对构造平稳分布 $\pi$ 的参数的置信区间提供了理论依据。

综合上述定理：

• 不可约+非周期+正常返 $\Rightarrow$ 唯一的平稳分布；
• LLN：Markov链的实值函数的遍历均值以概率1收敛到极限分布下的均值；
• CLT：遍历均值作合适的变化依分布收敛到标准正态分布。

因此，我们可以建立一个以 $\pi$ 为平稳分布的Markov链，则在运行此链足够长时间后，该Markov链就会达到平稳状态，此时Markov链的值相当于从分布 $\pi$ 中抽取样本。

MCMC算法的优点在于适用于广泛且困难的问题，此外问题维数的增加通常不会降低其收敛速度或使其更复杂。下面介绍两种重要的MCMC算法。

6.2 Metropolis-Hastings算法

M-H算法是一类常用的构造Markov链的方法。MCMC抽样方法的核心问题是确定满足要求（不可约、非周期、正常返且平稳分布为目标抽样分布）的一个马尔可夫链，因此确定从当前值转移到下一个值的规则十分重要。M-H抽样方法通过如下方式产生Markov链 $\{X_0,X_1,X_2,\dots\}$：

• 构造合适的提议分布 $g(\cdot\mid X_t)$.
• 从某个分布 $g$ 中产生 $X_0$.
• 重复（直到Markov链达到平稳状态）
  • 从 $g(\cdot\mid X_t)$ 中产生 $Y$.
  • 从 $U(0,1)$ 中产生 $U$.
  • 计算接受概率 $$\alpha(X_t,Y)=\min\left\{1,\dfrac{f(Y)g(X_t\mid Y)}{f(X_t)g(Y\mid X_t)}\right\}.$$
  • 如果 $U\leq\alpha(X_t,Y)$ 则接受 $Y$ 并令 $X_{t+1}=Y$，否则 $X_{t+1}=X_t$.
  • 增加 $t$，返回到循环内的第一步。

上述过程是有理论保证的。事实上，记 $q_{ij}=g(x_{t+1}=j\mid x_t=i)$，首先由M-H抽样方法产生的序列 $\{x_0,x_1,\dots\}$ 为一马尔可夫链，其转移概率为 $$p_{ij}=q_{ij}\alpha(i,j)+\delta_i(j)(1-r(i)),$$ 其中 $r(i)=\sum_jq_{ij}\alpha(i,j)$，$\alpha_{i,j}=\min\{1,f_jq_{ji}/f_iq_{ij}\}$，$\delta_i(j)=1$ 当且仅当 $i=j$。所以 $$f_ip_{ij}=f_iq_{ij}\min\left\{1,\dfrac{f_jq_{ji}}{f_iq_{ij}}\right\}=f_jq_{ji}\min\left\{1,\dfrac{f_iq_{ij}}{f_jq_{ji}}\right\}=f_jp_{ji},\quad i\neq j.$$ 而 $i=j$ 时，$f_i\delta_i(j)(1-r(i))=f_j\delta_j(i)(1-r(j))$，因此细致平稳方程满足，$f$ 即该链的平稳分布。进一步可以验证该链是不可约、正常返、非周期的。此外，如果状态空间是连续的，证明类似，只需将转移概率修改为转移核即可。

注意，提议分布 $g$ 的选择要使得产生的马尔可夫链满足不可约、正常返、非周期且具有平稳分布 $f$ 等正则化条件外，还应满足：

• 提议分布的支撑集包含目标分布的支撑集；
• 容易从提议分布中抽样，常取为已知的分布，如正态分布或 $t$ 分布等；
• 提议分布应使接受概率容易计算；
• 提议分布的尾部要比目标分布的尾部厚。

此外，Robert和Casella（1999）指出，接受概率并非越大越好，因为这可能导致较慢的收敛性；Gelman等（1996）建议当参数维数为1时，接受概率应略小于0.5，当参数维数大于5时，接受概率应降至0.25左右。

例如，使用M-H抽样方法从瑞利（Rayleigh）分布中抽样，其密度函数为 $$f(x)=\dfrac{x}{\sigma^2}\text{e}^{-\frac{x^2}{2\sigma^2}},\quad x\geq0,\sigma>0.$$ 瑞利分布通常被用来建模具有快速衰减的生存时间，因为风险率 $f(x)/S(x)=x/\sigma^2$ 是线性增长的。

使用M-H抽样方法如下：

• 令 $g(\cdot\mid X)$ 为 $\chi^2(X)$；
• 从 $\chi^2(1)$ 中产生 $X_1$；
• 对 $i=2,\dots,N$，重复：
  • 从 $\chi^2(X_{i-1})$ 中产生 $Y$；
  • 产生 $U\sim U(0,1)$；
  • 计算 $$r(X_{i-1},Y)=\dfrac{f(Y)g(X_{i-1}\mid Y)}{f(X_{i-1})g(Y\mid X_{i-1})},$$
  • 若 $U\leq r(X_{i-1},Y)$，则接受 $Y$，令 $X_i=Y$；否则令 $X_i=X_{i-1}$。
  • 增加 $i$，返回到循环内的第一步。

f <- function(x, sigma) {
  if (any(x < 0)) return(0)
  return(x / sigma^2 * exp(-x^2 / (2 * sigma^2)))
}

m <- 10000
sigma <- 4
x <- numeric(m)
x[1] <- rchisq(1, 1)
k <- 0
u <- runif(m)
for (i in 2:m) {
  y <- rchisq(1, x[i-1])
  r <- f(y, sigma) * dchisq(x[i-1], y) / f(x[i-1], sigma) / dchisq(y, x[i-1])
  if (u[i] <= r) { x[i] <- y }
  else { x[i] <- x[i-1]; k <- k + 1 }
}
print(k / m)

## [1] 0.4016

可以看到约有40.16%的候选点被拒绝了，可见算法效率不是特别高。

index <- 5000:5500
plot(index, x[index], type = "l", main = "", ylab = "x")

b <- 2000 # burnin period
y <- x[(b+1):m]
a <- ppoints(100)
QR <- sigma * sqrt(-2 * log(1 - a))
Q <- quantile(x, a)

par(mfrow = c(1,2))
hist(y, breaks = "scott", xlab = "sample", freq = FALSE)
lines(QR, f(QR, sigma))
qqplot(QR, Q, xlab = "Rayleigh Quantiles", ylab = "Sample Quantiles")
abline(0, 1, col = 2, lwd = 2)

从分布拟合的角度看，预烧期2000之后的样本基本已经达到平稳分布，拟合效果较好。

根据提议分布 $g$ 的不同选择，M-H抽样方法衍生出了几个不同的变种。

6.2.1 Metropolis抽样方法

在Metropolis抽样方法中，提议分布是对称的，即 $g(\cdot\mid X_n)$ 满足 $g(X\mid Y)=g(Y\mid X)$，因此接受概率为 $$\alpha(X_t,Y)=\min\left\{1,\dfrac{f(Y)}{f(X_t)}\right\}.$$

6.2.2 随机游动Metropolis抽样方法

随机游动Metropolis抽样方法是Metropolis方法的一个应用例子。假设候选点 $Y$ 从一个对称的提议分布 $g(Y\mid X_t)=g(\vert X_t-Y\vert)$ 中产生，则在每一次迭代中，从 $g(\cdot)$ 中产生一个随机增量 $Z$，然后 $Y=X_t+Z$。比如增量 $Z$ 可以从正态分布中产生，此时 $Y\mid X_t\sim N(X_t,\sigma^2)$，$\sigma^2>0$。

例如，使用提议分布 $N(X_t,\sigma^2)$ 和随机游动Metropolis算法产生自由度为 $\nu$ 的 $t$ 分布随机数，并对不同方差 $\sigma^2$ 重复此过程。记 $f(x)$ 为 $t_\nu$ 的密度函数，则接受概率为 $$\alpha(X_t,Y)=\min\left\{1,\dfrac{f(Y)}{f(X_t)}\right\}.$$

rw.Metropolis <- function(n, sigma, x0, N) {
  x <- numeric(N)
  x[1] <- x0
  u <- runif(N)
  k <- 0
  for (i in 2:N) {
    y <- rnorm(1, x[i-1], sigma)
    if (u[i] <= dt(y, n) / dt(x[i-1], n)) { x[i] <- y }
    else { x[i] <- x[i-1]; k <- k + 1 }
  }
  return(list(x = x, k = k))
}

n <- 4
N <- 2000
sigma <- c(0.05, 0.5, 2, 16)
x0 <- 25
rw1 <- rw.Metropolis(n, sigma[1], x0, N)
rw2 <- rw.Metropolis(n, sigma[2], x0, N)
rw3 <- rw.Metropolis(n, sigma[3], x0, N)
rw4 <- rw.Metropolis(n, sigma[4], x0, N)
print(c(rw1$k, rw2$k, rw3$k, rw4$k) / N)

## [1] 0.0035 0.0945 0.4605 0.8985

可以看到只有第三个链的拒绝率在区间 $[0.15,0.5]$ 内。

par(mfrow = c(2, 2))
refline <- qt(c(0.025, 0.975), n)
rw <- cbind(rw1$x, rw2$x, rw3$x,  rw4$x)
for (j in 1:4) {
  plot(rw[,j], type = "l", xlab = bquote(sigma == .(round(sigma[j],3))),
       ylab = "X", ylim=range(rw[,j]))
  abline(h = refline)
}

par(mfrow = c(1,1))

可以看到，当 $\sigma=0.05$ 时增量太小，几乎每个候选点都被接受了，链在迭代2000次后还没有收敛；当 $\sigma=0.5$ 时，链的收敛较慢；当 $\sigma=2$ 时，链的收敛很快；当 $\sigma=16$ 时，接受的概率太小，链虽然收敛了，但是效率很低。

a <- c(0.05, seq(0.1, 0.9, 0.1), 0.95)
Q <- qt(a, n)
mc <- rw[501:N,]
Qrw <- apply(mc, 2, function(x) quantile(x, a))
QQ <- data.frame(round(cbind(Q, Qrw), 3))
names(QQ) <- c('True', 'sigma=0.05', 'sigma=0.5', 'sigma=2', 'sigma=16')
knitr::kable(QQ)

	True	sigma=0.05	sigma=0.5	sigma=2	sigma=16
5%	-2.132	21.639	-1.404	-2.144	-2.376
10%	-1.533	21.848	-0.945	-1.521	-1.370
20%	-0.941	22.132	-0.399	-0.860	-0.694
30%	-0.569	22.420	0.128	-0.535	-0.482
40%	-0.271	22.587	0.596	-0.267	-0.065
50%	0.000	22.790	1.160	0.005	0.272
60%	0.271	23.207	2.157	0.232	0.494
70%	0.569	23.510	8.362	0.535	0.923
80%	0.941	23.766	16.147	0.959	1.378
90%	1.533	24.073	20.063	1.650	2.168
95%	2.132	24.718	26.414	2.377	2.936

可以看到第三个链的分位数和真实分布的分位数也是十分接近的。

6.2.3 独立抽样方法

独立抽样是M-H抽样方法的另一个特殊情形，其中的提议分布不依赖于链的前一步状态值，因此 $g(Y\mid X_t)=g(Y)$，此时接受概率为 $$\alpha(X_t,Y)=\min\left\{1,\dfrac{f(Y)g(X_t)}{f(X_t)g(Y)}\right\}.$$

例如，假设从一个正态混合分布 $pN(\mu_1,\sigma_1^2)+(1-p)N(\mu_2,\sigma_2^2)$ 中观测到一个样本 $\pmb z=(z_1,\dots,z_n)$，利用独立抽样方法求 $p$ 的估计。在 $p$ 无先验信息可以利用时，考虑先验分布 $\pi(p)$ 为 $U(0,1)$，则 $p$ 的后验分布为 $$\pi(p\mid\pmb z)\propto\prod_{j=1}^n[pf_1(z_j\mid p)+(1-p)f_2(z_j\mid p)].$$

其中 $f_1,f_2$ 分别为两个正态的密度。考虑提议分布为贝塔分布 $Be(a,b)$，其密度函数为 $g(y)\propto y^{a-1}(1-y)^{b-1}$，接受概率为 $$\alpha(X_t,Y)=\min\left\{1,\dfrac{\pi(Y\mid\pmb z)g(X_t)}{\pi(X_t\mid\pmb z)g(Y)}\right\}.$$

set.seed(12345)
m <- 10000
xt <- numeric(m)
a <- 1
b <- 1
p <- 0.2
n <- 100
mu <- c(0, 5)
sigma <- c(1, 1)
i <- sample(1:2, size = n, replace = TRUE, prob = c(p, 1-p))
x <- rnorm(n, mu[i], sigma[i])

u <- runif(m)
y <- rbeta(m, a, b)
xt[1] <- 0.5
for (i in 2:m) {
  fy <- y[i] * dnorm(x, mu[1], sigma[1]) + (1 - y[i]) * dnorm(x, mu[2], sigma[2])
  fx <- xt[i-1] * dnorm(x, mu[1], sigma[1]) + (1 - xt[i-1]) * dnorm(x, mu[2], sigma[2])
  r <- prod(fy / fx) * dbeta(xt[i-1], a, b) / dbeta(y[i], a, b)
  if (u[i] <= r) { xt[i] <- y[i] }
  else { xt[i] <- xt[i-1] }
}

par(mfrow = c(1, 2))
plot(xt, type = "l", ylab = "p")
abline(h = p, v = 101, lty = 3)
hist(xt[-(1:1000)], main = "", xlab = "p", prob = TRUE)
z <- seq(min(xt), max(xt), length = 100)
lines(z, dnorm(z, mean(xt), sd(xt)))
abline(v = mean(xt), col = 2, lwd = 2)

可以看到 $p$ 的后验期望接近0.2（真值）。

6.3 Gibbs算法

Gibbs抽样方法特别适合于目标分布是多元的场合，其最令人感兴趣的方面是为了产生不可约、非周期、并以高维空间的目标分布作为其平稳分布的马尔可夫链，只需要从一些一元分布（全条件分布）中进行抽样就可以了。即将从多元目标分布中抽样转化为从一元目标分布抽样，这是Gibbs抽样的重要性所在。

令 $\pmb X=(X_1,\dots,X_p)\in\mathbb R^p$ 为随机变量，其联合分布 $f(\pmb x)$ 为目标分布。记 $$\pmb X_{-i}=(X_1,\dots,X_{i-1},X_{i+1},\dots,X_p)\in\mathbb R^{p-1},$$ 并记 $X_i\mid\pmb X_{-i}$ 的满（全）条件密度为 $f(x_i\mid\pmb x_{-i})$ $(i=1,\dots,p)$。Gibbs抽样方法是从这 $p$ 个条件分布中产生候选点，以解决直接从 $f$ 中进行抽样的困难。Gibbs抽样从一个初始值 $\pmb x^{(0)}=(x_1^{(0)},\dots,x_p^{(0)})$ 出发，按下述过程迭代更新分量： $$\begin{align*} x_1^{(t)}&\sim f(x_1\mid x_2^{(t-1)},x_3^{(t-1)},\dots,x_p^{(t-1)}),\\ x_2^{(t)}&\sim f(x_2\mid x_1^{(t)},x_3^{(t-1)},\dots,x_p^{(t-1)}),\\ &\vdots\\ x_i^{(t)}&\sim f(x_i\mid x_1^{(t)},\dots,x_{i-1}^{(t)},x_{i+1}^{(t-1)},\dots,x_p^{(t-1)}),\\ &\vdots\\ x_p^{(t)}&\sim f(x_p\mid x_1^{(t)},\dots,x_{p-1}^{(t)}). \end{align*}$$

有时候对所有分量完成一次更新比较复杂，可以采用随机化的方式，每次仅更新一个分量（随机选择），而其它分量保持不变。注意，这里无需计算每一个候选点的接受概率 $\alpha$，全部接受。

例如，使用Gibbs抽样方法产生二元正态分布 $N(\mu_1,\mu_2,\sigma_1^2,\sigma_2^2,\rho)$ 的随机数。注意到 $$\begin{align*} X_1\mid X_2=x_2&\sim N\left(\mu_1+\rho\dfrac{\sigma_1}{\sigma_2}(x_2-\mu_2),(1-\rho^2)\sigma_1^2\right),\\ X_2\mid X_1=x_1&\sim N\left(\mu_2+\rho\dfrac{\sigma_2}{\sigma_1}(x_1-\mu_1),(1-\rho^2)\sigma_2^2\right). \end{align*}$$

于是通过如下代码就可以实现目的。

m <- 10000
burnin <- 1000
X <- matrix(0, m, 2)
rho <- 0.8
mu <- c(0, 2)
sigma <- c(1, 0.5)
s <- sqrt(1 - rho^2) * sigma
X[1,] <- mu
for (i in 2:m) {
  x2 <- X[i-1, 2]
  m1 <- mu[1] + rho * (x2 - mu[2]) * sigma[1] / sigma[2]
  X[i,1] <- rnorm(1, m1, s[1])
  x1 <- X[i, 1]
  m2 <- mu[2] + rho * (x1 - mu[1]) * sigma[2] / sigma[1]
  X[i,2] <- rnorm(1, m2, s[2])
}
x <- X[-(1:burnin),]
plot(x[,1], x[,2], xlab = bquote(X[1]), ylab = bquote(X[2]),
     main = "", cex = 0.5, ylim = range(x[,2]))

colMeans(x)

## [1] 0.004515394 2.005059651

cov(x)

##           [,1]     [,2]
## [1,] 0.9821537 0.396883
## [2,] 0.3968830 0.248386

cor(x)

##           [,1]      [,2]
## [1,] 1.0000000 0.8035431
## [2,] 0.8035431 1.0000000

可以看到，模拟的数值与真值十分接近。

6.4 收敛监视

较主观的监视收敛的方法有

• 监视蒙特卡洛误差：较小的MC误差表明在计算感兴趣的量时精度较高，即MC误差越小表明马尔可夫链的收敛性越好。
• 样本路径图：将马尔可夫链生成的值对迭代次数作图，如果所有的值都在一个区域里没有明显的周期性或趋势性，则认为收敛性已经达到。
• 累积均值图：将马尔可夫链的累积均值对迭代次数作图，如果累积均值在经过一些迭代后基本稳定，则认为收敛性已经达到。
• 自相关函数（ACF）图：将马尔可夫链生成样本的自相关函数对迭代次数作图，较低或较高的自相关性分别表明了收敛性的好与坏。

除此之外，下面介绍一种基于统计量的方法。

6.4.1 Gelman-Rubin方法

假设有 $k$ 个链，每个链有 $n$ 个样本，感兴趣的量记为 $\phi$，其在目标分布下有期望 $\mu$ 和方差 $\sigma^2$。记 $\phi_{jt}$ 表示到链 $j$ 的第 $t$ 个样本时 $\phi$ 的值，那么在混合样本中，$\mu$ 的一个无偏估计为 $\hat{\mu}=\overline{\phi}_{\cdot\cdot}$，而链之间的方差 $B_n$ 和链内的方差 $W_n$ 分别为 $$B_n=\dfrac{n}{k-1}\sum_{j=1}^k(\overline{\phi}_{j\cdot}-\overline{\phi}_{\cdot\cdot})^2,\quad W_n=\dfrac{1}{k(n-1)}\sum_{j=1}^k\sum_{t=1}^n(\phi_{jt}-\overline{\phi}_{j\cdot})^2.$$

从而可以使用 $B_n$ 和 $W_n$ 的加权估计 $\sigma^2$： $$\hat{\sigma}^2=\dfrac{n-1}{n}W_n+\dfrac{1}{n}B_n.$$

如果初始值是从目标分布中抽取的，那么 $\hat{\sigma}^2$ 是 $\sigma^2$ 的无偏估计，但如果初始值过度分散，则会高估 $\sigma^2$。

一个自然的想法就是令 $R=\hat{\sigma}^2/\sigma^2$，则称 $\sqrt{R}$ 为尺度缩减因子，它的估计为 $$\hat{R}=\dfrac{\hat{\sigma}^2}{W_n},$$ 称 $\sqrt{\hat{R}}$ 为潜在尺度缩减因子。当链达到收敛并且产生的数据很大时，$\hat{R}$ 应该趋于1，一般而言，当 $\hat{R}<1.1$ 或 $1.2$ 时就认为已达到收敛。

例如，目标分布为 $N(0,1)$，提议分布为 $N(X_t,\sigma^2)$，$\phi_{jt}$ 表示第 $j$ 个链前 $t$ 个样本的平均。

# 计算潜在尺度缩减因子
Gelman.Rubin <- function(phi) {
  phi <- as.matrix(phi)
  k <- nrow(phi); n <- ncol(phi)
  phi.means <- rowMeans(phi)
  B <- n * var(phi.means)
  phi.w <- apply(phi, 1, var)
  W <- mean(phi.w)
  v.hat <- W * (n - 1) / n + B / n
  r.hat <- v.hat / W
  return(r.hat)
}

# 采用M-H算法生成链
normal.chain <- function(sigma, N, X1) {
  x <- rep(0, N)
  x[1] <- X1
  u <- runif(N)
  for (i in 2:N) {
    xt <- x[i - 1]
    y <- rnorm(1, xt, sigma)
    r1 <- dnorm(y, 0, 1) * dnorm(xt, y, sigma)
    r2 <- dnorm(xt, 0, 1) * dnorm(y, xt, sigma)
    if (u[i] <= r1 / r2) { x[i] <- y} else { x[i] <- xt }
  }
  return(x)
}

set.seed(12345)
sigma <- 1
k <- 4 # numebr of chains
n <- 15000 # length of chains
b <- 1000 # burnin period
x0 <- c(-10, -5, 5, 10) # initial value
X <- matrix(nrow = k, ncol = n)
for (i in 1:k) {
  X[i, ] <- normal.chain(sigma, n, x0[i])
}

# 绘制累积均值(phi)图
phi <- t(apply(X, 1, cumsum))
for (i in 1:nrow(phi)) {
  phi[i,] <- phi[i,] / (1:ncol(phi))
}
for (i in 1:k) {
  if (i == 1) {
    plot((b+1):n, phi[i, (b+1):n], ylim = c(-0.2, 0.2),
         type = "l", xlab = 'Index', ylab = bquote(phi))
  } else { lines(phi[i, (b+1):n], col = i) }
}

# 绘制潜在尺度缩减因子图
rhat <- rep(0, n)
for (j in (b+1):n) {
  rhat[j] <- Gelman.Rubin(phi[, 1:j])
}
plot(rhat[(b+1):n], type = "l", xlab = "", ylab = "R")
abline(h = 1.1, lty = 2)