# 马尔可夫链的收敛性# 马尔可夫链简介马尔可夫链是一个随机变量序列 θ ( t ) \theta^{(t)} θ ( t )
p ( θ ( t ) ∣ θ ( t − 1 ) , ⋯ , θ ( 0 ) ) = p ( θ ( t ) ∣ θ ( t − 1 ) ) p(\theta^{(t)}\mid \theta^{(t-1)},\cdots, \theta^{(0)})=p(\theta^{(t)}\mid \theta^{(t-1)})
p ( θ ( t ) ∣ θ ( t − 1 ) , ⋯ , θ ( 0 ) ) = p ( θ ( t ) ∣ θ ( t − 1 ) )
其中，p p p 转移分布 。
马尔可夫链的状态空间 是其所有可能状态的集合。对于有限状态的马尔可夫链，状态转移可以用状态转移矩阵 P P P p i j p_{ij} p ij i i i j j j
当前状态向量 ：π i ( t ) = P ( θ ( t ) = i ) \pi_i^{(t)}=P(\theta^{(t)}=i) π i ( t ) = P ( θ ( t ) = i ) 状态转移方程 ：π ( t ) = π ( t − 1 ) P = π ( 0 ) P t \pi^{(t)}=\pi^{(t-1)}P=\pi^{(0)}P^t π ( t ) = π ( t − 1 ) P = π ( 0 ) P t
# 平稳分布与收敛条件平稳分布 （Stationary Distribution）是指马尔可夫链在经过足够长的时间后，其状态分布不再变化的分布，满足 π = π P \pi = \pi P π = π P
为了确保马尔可夫链能够收敛到唯一的平稳分布，并让该平稳分布成为极限分布，需要满足以下条件：
不可约性 （Irreducibility）：对于状态空间中的任意两个状态 i i i j j j t i j t_{ij} t ij i i i j j j P ( θ ( t i j ) = j ∣ θ ( 0 ) = i ) > 0 P(\theta^{(t_{ij})}=j\mid\theta^{(0)}=i)>0 P ( θ ( t ij ) = j ∣ θ ( 0 ) = i ) > 0
非周期性 （Aperiodicity）：状态 i i i k i k_i k i i i i k i = g c d { t : P ( θ ( t ) = i ∣ θ ( 0 ) = i ) > 0 } k_i=gcd\{t:P(\theta^{(t)}=i\mid\theta^{(0)}=i)>0\} k i = g c d { t : P ( θ ( t ) = i ∣ θ ( 0 ) = i ) > 0 } k i = 1 k_i=1 k i = 1 i i i
常返性 （Recurrence）：
设 T i T_i T i i i i P ( T i < ∞ ) = 1 P(T_i < \infty) = 1 P ( T i < ∞ ) = 1 i i i 常返的 。
若常返态 i i i E [ T i ] E[T_i] E [ T i ] 正常返的 ；否则为零常返的 。
对不可约马尔可夫链，如果存在一个正常返（或零常返）的状态，则所有状态都正常返（或零常返）。
若马尔可夫链的状态 i i i lim t → ∞ π i ( t ) = 1 / E [ T i ] \lim_{t \to \infty} \pi_i^{(t)} = 1/E[T_i] lim t → ∞ π i ( t ) = 1/ E [ T i ]
结论 ：在有限状态空间中：
不可约马尔可夫链具有唯一的平稳分布 π \pi π
如果不可约马尔可夫链同时也是非周期的，那么对于任意初始状态 π ( 0 ) \pi^{(0)} π ( 0 ) π \pi π lim t → ∞ π ( t ) = π \lim_{t \to \infty}\pi^{(t)}=\pi lim t → ∞ π ( t ) = π 极限分布 。
如果不可约非周期马尔可夫链是正常返的，则它存在唯一的平稳分布 π \pi π π i = 1 / E [ T i ] \pi_i=1/E[T_i] π i = 1/ E [ T i ]
遍历性 （Ergodicity）：如果一个马尔可夫链是不可约、非周期和正常返的，则称其为遍历的 。遍历性是马尔可夫链收敛到唯一平稳分布并使其成为极限分布的充分条件。
对于连续状态空间 的马尔可夫链，我们通常讨论其在特定集合上的概率。
示例 ：AR(1) 模型
θ ( t ) = μ + ρ ( θ ( t − 1 ) − μ ) + ϵ t , ϵ t ∼ N ( 0 , σ 2 ) \theta^{(t)}=\mu+\rho(\theta^{(t-1)}-\mu)+\epsilon_t,~\epsilon_t\sim N(0,\sigma^2)
θ ( t ) = μ + ρ ( θ ( t − 1 ) − μ ) + ϵ t , ϵ t ∼ N ( 0 , σ 2 )
即：
θ ( t ) ∣ θ ( t − 1 ) ∼ N ( μ + ρ ( θ ( t − 1 ) − μ ) , σ 2 ) \theta^{(t)}\mid\theta^{(t-1)}\sim N(\mu+\rho(\theta^{(t-1)}-\mu),\sigma^2)
θ ( t ) ∣ θ ( t − 1 ) ∼ N ( μ + ρ ( θ ( t − 1 ) − μ ) , σ 2 )
当 ∣ ρ ∣ < 1 |\rho|<1 ∣ ρ ∣ < 1
θ ( t ) ∼ N ( μ , σ 2 / ( 1 − ρ 2 ) ) \theta^{(t)}\sim N(\mu,\sigma^2/(1-\rho^2))
θ ( t ) ∼ N ( μ , σ 2 / ( 1 − ρ 2 ))
# MCMC 的有效性验证MCMC （Markov Chain Monte Carlo）方法的基本思想是：构造一个马尔可夫链，使其平稳分布为我们想要采样的目标分布 p ( θ ∣ y ) p(\theta \mid y) p ( θ ∣ y )
细致平衡条件 （Detailed Balance Condition, DBC）：p ( ϕ ∣ θ ) p(\phi \mid \theta) p ( ϕ ∣ θ )
π ( θ ) p ( ϕ ∣ θ ) = π ( ϕ ) p ( θ ∣ ϕ ) \pi(\theta)p(\phi\mid\theta)=\pi(\phi)p(\theta\mid\phi)
π ( θ ) p ( ϕ ∣ θ ) = π ( ϕ ) p ( θ ∣ ϕ )
则 π ( θ ) \pi(\theta) π ( θ ) p ( θ ∣ y ) p(\theta \mid y) p ( θ ∣ y )
# 马尔可夫链的构造与诊断# 初始值选择初始值的选择会影响 MCMC 采样的收敛速度。虽然理论上马尔可夫链最终会收敛到平稳分布，但如果初始值离目标分布很远，收敛会非常慢。
合理初值 ：任何不介意包含在样本中的点都是一个好的起点。
获得好初值的方法 ：
多点老化 （Burn-in）：从多个不同的起始点开始采样，并丢弃开始的若干次迭代（例如，丢弃前 X 次迭代）。从最大似然估计（MLE）开始 ：arg max θ p ( y ∣ θ ) \arg\max_\theta p(y \mid \theta) arg max θ p ( y ∣ θ ) 从最大后验概率（MAP）开始 ：arg max θ p ( θ ∣ y ) \arg\max_\theta p(\theta \mid y) arg max θ p ( θ ∣ y ) 选择多个起始点 （例如 5-10 个）运行多条链，以确保覆盖整个样本空间。
# 样本大小与收敛诊断为了确保 MCMC 采样的有效性，我们需要判断链是否已收敛以及需要多少有效样本。
有效样本量 ：通常 100-2000 个有效样本就足够了。如果需要考察极值，可能需要更多。
# Gelman-Rubin 势尺度缩减因子（Potential Scale Reduction Factor）Gelman-Rubin 诊断用于判断多条马尔可夫链是否已收敛。其核心思想是比较链内方差与链间方差。
步骤 ：
选择 m m m m m m n n n ψ i j \psi_{ij} ψ ij
计算链内方差 W W W B B B
单链均值：ψ ˉ . j = ∑ i = 1 n ψ i j / n \bar{\psi}_{.j}=\sum_{i=1}^n\psi_{ij}/n ψ ˉ . j = ∑ i = 1 n ψ ij / n
总体均值：ψ ˉ . . = ∑ j = 1 m ψ ˉ . j / m \bar{\psi}_{..}=\sum_{j=1}^m\bar{\psi}_{.j}/m ψ ˉ .. = ∑ j = 1 m ψ ˉ . j / m
单链方差：s j 2 = ∑ i = 1 n ( ψ i j − ψ ˉ . j ) 2 / ( n − 1 ) s_j^2=\sum_{i=1}^n(\psi_{ij}-\bar{\psi}_{.j})^2/(n-1) s j 2 = ∑ i = 1 n ( ψ ij − ψ ˉ . j ) 2 / ( n − 1 )
链间方差：B = n m − 1 ∑ j = 1 m ( ψ ˉ . j − ψ ˉ . . ) 2 B=\frac{n}{m-1}\sum_{j=1}^m (\bar{\psi}_{.j}-\bar{\psi}_{..})^2 B = m − 1 n ∑ j = 1 m ( ψ ˉ . j − ψ ˉ .. ) 2
链内方差：W = 1 m ∑ j = 1 m s j 2 W=\frac{1}{m}\sum_{j=1}^ms_j^2 W = m 1 ∑ j = 1 m s j 2
计算估计量的边缘后验方差 V a r ^ + ( ψ ∣ y ) = n − 1 n W + 1 n B \widehat{Var}^+(\psi \mid y) = \frac{n-1}{n}W+\frac{1}{n}B Va r + ( ψ ∣ y ) = n n − 1 W + n 1 B
计算势尺度缩减因子 R ^ ψ = V a r ^ + ( ψ ∣ y ) W \hat{R}_\psi = \sqrt{\frac{\widehat{Var}^+(\psi \mid y)}{W}} R ^ ψ = W Va r + ( ψ ∣ y )
判断 ：如果 R ^ ψ \hat{R}_\psi R ^ ψ
注意事项 ：
多峰分布 ：很小的 R ^ ψ \hat{R}_\psi R ^ ψ 适用性 ：该方法基于均值和方差，更适用于近似正态分布的变量。对于非正态分布的变量，进行变换后再诊断效果更好（如对正值变量取对数）。
# Geweke 检验Geweke 检验通过比较马尔可夫链首尾部分的均值来诊断收敛性。如果检验结果表明首尾两部分来自不同的分布，则可能需要丢弃更多的前段样本，以确保链的剩余部分已收敛。
# 有效样本量MCMC 估计量是马尔可夫链产生的样本的均值。根据强大数定律（SLLN） ，在正则条件下，样本均值 h ^ T = 1 T ∑ t = 1 T h ( θ ( t ) ) \hat{h}_T = \frac{1}{T}\sum_{t=1}^T h(\theta^{(t)}) h ^ T = T 1 ∑ t = 1 T h ( θ ( t ) ) E [ h ( θ ) ∣ y ] E[h(\theta) \mid y] E [ h ( θ ) ∣ y ]
根据中心极限定理（CLT） ，在强正则条件下，样本均值近似服从正态分布：
h ^ T ∼ N ( E [ h ( θ ) ∣ y ] , σ 2 / T ) \hat{h}_T\sim N(E[h(\theta)\mid y],~\sigma^2/T)
h ^ T ∼ N ( E [ h ( θ ) ∣ y ] , σ 2 / T )
其中，σ 2 \sigma^2 σ 2 σ 2 = V a r ( h ( θ ) ∣ y ) ( 1 + 2 ∑ k = 1 ∞ ρ k ) \sigma^2=Var(h(\theta)\mid y)\left(1+2\sum_{k=1}^\infty\rho_k\right) σ 2 = Va r ( h ( θ ) ∣ y ) ( 1 + 2 ∑ k = 1 ∞ ρ k ) ρ k \rho_k ρ k k k k
有效样本量 （Effective Sample Size, n e f f n_{eff} n e ff
n e f f = m n 1 + 2 ∑ t = 1 ∞ ρ t n_{eff}=\frac{mn}{1+2\sum_{t=1}^\infty\rho_t}
n e ff = 1 + 2 ∑ t = 1 ∞ ρ t mn
通常用估计的自相关系数 ρ ^ t \hat{\rho}_t ρ ^ t
n e f f ≈ m n 1 + 2 ∑ t = 1 T ρ ^ t n_{eff}\approx\frac{mn}{1+2\sum_{t=1}^T\hat{\rho}_t}
n e ff ≈ 1 + 2 ∑ t = 1 T ρ ^ t mn
其中，ρ ^ t = 1 − V t 2 v a r ^ + \hat{\rho}_t=1-\frac{V_t}{2\widehat{var}^+} ρ ^ t = 1 − 2 v a r + V t V t = 1 m ( n − t ) ∑ j = 1 m ∑ i = t + 1 n ( ψ i , j − ψ i − t , j ) 2 V_t=\frac{1}{m(n-t)}\sum_{j=1}^m\sum_{i=t+1}^n(\psi_{i,j}-\psi_{i-t,j})^2 V t = m ( n − t ) 1 ∑ j = 1 m ∑ i = t + 1 n ( ψ i , j − ψ i − t , j ) 2
