马尔科夫链蒙特卡罗方法

马尔科夫链蒙特卡罗 ( M C M C ) (MCMC) (MCMC) 方法是一组从后验分布（未知其标准化常量）中渐进地产生样本的非常有效的方法，它通过形成样本平均值得到关于模型参数的推论。 ( M C M C ) (MCMC) (MCMC) 方法可以分析复杂的问题，而不再要求用户为能使用解析方法而强行把问题归在某一简单的框架之下。

1. 吉布斯采样器

吉布斯采样器是一种非常通用的 ( M C M C ) (MCMC) (MCMC) 方法，吉布斯采样器和一种更通用的 M e t r o p o l i s − H a s t i n g s Metropolis-Hastings Metropolis−Hastings 算法共同形成了多种情况的分析基础。

设用 f ( θ ) f(\theta) f(θ) 表示后验分布，我们希望从中抽取样本； θ \theta θ 是 d d d 维参数向量 ( θ 1 , . . . , θ d ) T (\theta_1,...,\theta_d)^T (θ1​,...,θd​)T。很可能无法确切的获知 f ( θ ) f(\theta) f(θ) ，但已知函数 g ( θ ) g(\theta) g(θ) ，而 f ( θ ) = g ( θ ) / ∫ g ( θ ′ ) d θ ′ f(\theta)=g(\theta)/\int g(\theta')d\theta' f(θ)=g(θ)/∫g(θ′)dθ′ 。例如似然函数和先验密度的后验密度函数，可能未知标准化常量 p ( D ) = ∫ p ( D ∣ θ ) p ( θ ) p(D)=\int p(D|\theta)p(\theta) p(D)=∫p(D∣θ)p(θ) ，这时认为 g ( θ ) = p ( D ∣ θ ) p ( θ ) g(\theta)=p(D|\theta)p(\theta) g(θ)=p(D∣θ)p(θ)。

吉布斯采样算法

在吉布斯采样算法中，给定了所有先前值 θ 0 , θ 1 , . . . , θ t − 1 , θ t \theta^0,\theta^1,...,\theta^{t-1},\theta^t θ0,θ1,...,θt−1,θt 的分布仅取决于 θ t − 1 \theta^{t-1} θt−1，这体现为马尔可夫性，由此产生的序列称为马尔科夫序列。

• 从先验/后验分布中选取 θ \theta θ 的任意初始值 θ 0 = ( θ 1 0 , . . . , θ d 0 ) T \theta_0=(\theta^0_1,...,\theta^0_d)^T θ0​=(θ10​,...,θd0​)T。
• 在迭代的第 t + 1 t+1 t+1 步：
• 从 f ( θ 1 ∣ θ 2 t , . . . , θ d t ) f(\theta_1|\theta_2^t,...,\theta_d^t) f(θ1​∣θ2t​,...,θdt​) 中抽取样本 θ 1 t + 1 \theta_1^{t+1} θ1t+1​；
• 从 f ( θ 2 ∣ θ 1 t , θ 3 t , . . . , θ d t ) f(\theta_2|\theta_1^t,\theta_3^t,...,\theta_d^t) f(θ2​∣θ1t​,θ3t​,...,θdt​) 中抽取样本 θ 2 t + 1 \theta_2^{t+1} θ2t+1​；
• 从 f ( θ 3 ∣ θ 1 t , θ 2 t , θ 4 t . . . , θ d t ) f(\theta_3|\theta_1^t,\theta_2^t,\theta_4^t...,\theta_d^t) f(θ3​∣θ1t​,θ2t​,θ4t​...,θdt​) 中抽取样本 θ 3 t + 1 \theta_3^{t+1} θ3t+1​；
• 对所有的变量继续这种操作，直到最后从 f ( θ d ∣ θ 1 t + 1 , . . . , θ d − 1 t + 1 ) f(\theta_d|\theta_1^{t+1},...,\theta_{d-1}^{t+1}) f(θd​∣θ1t+1​,...,θd−1t+1​)；
• θ t + 1 = ( θ 1 t + 1 , . . . , θ d t + 1 ) T \theta^{t+1}=(\theta_1^{t+1},...,\theta_d^{t+1})^T θt+1=(θ1t+1​,...,θdt+1​)T。
• 经过大量的迭代，向量 θ t \theta^t θt 就如同从联合密度 f ( θ ) f(\theta) f(θ) 中随机抽取一样。

经过足够步骤的迭代之后，样本 { θ t } \lbrace \theta^t \rbrace {θt} 将成为来自后验分布 f ( θ ) f(\theta) f(θ) 的相依样本。函数 h ( θ ) h(\theta) h(θ) 的期望值可用下式近似：
E [ h ( θ ) ] ≈ 1 N − M ∑ t = M + 1 N h ( θ t ) E[h(\theta)]\approx \frac{1}{N-M}\sum_{t=M+1}^Nh(\theta^t) E[h(θ)]≈N−M1​t=M+1∑N​h(θt)
其中， N N N 是吉布斯采样器的迭代次数， M M M 是老化期间被摒弃的初始样本的个数。

• 举例如何得出多元正态分布的条件分布

  设 X X X 是取自多元正态分布 N ( X ∣ μ , Σ ) N(X|\mu,\Sigma) N(X∣μ,Σ) 的 d d d 维向量，我们希望确定条件分布 X 1 ∣ X 2 X_1|X_2 X1​∣X2​ ，其中 X 1 X_1 X1​ 是 X X X 的前 d 1 d_1 d1​ 个分量组成的向量， X 2 X_2 X2​ 是 X X X 的后 ( d − d 1 ) (d-d_1) (d−d1​) 个分量组成的向量。

  把均值向量和协方差矩阵以分块的形式写成
  μ = ( μ 1 μ 2 ) , Σ = ( Σ 1 , 1 Σ 1 , 2 Σ 2 , 1 Σ 2 , 2 ) \mu=\begin{pmatrix} \mu_1\\\mu_2 \end{pmatrix},\Sigma=\begin{pmatrix} \Sigma_{1,1}&\Sigma_{1,2} \\ \Sigma_{2,1}&\Sigma_{2,2}\end{pmatrix} μ=(μ1​μ2​​),Σ=(Σ1,1​Σ2,1​​Σ1,2​Σ2,2​​)
  其中， μ 1 \mu_1 μ1​为 d 1 d_1 d1​ 维向量， μ 2 \mu_2 μ2​ 为 ( d − d 1 ) (d-d_1) (d−d1​) 维向量， Σ 1 , 1 \Sigma_{1,1} Σ1,1​ 为 d 1 ∗ d 1 d_1*d_1 d1​∗d1​ 维矩阵， Σ 1 , 1 \Sigma_{1,1} Σ1,1​ 为 d 1 ∗ （ d − d 1 ) d_1*（d-d_1) d1​∗（d−d1​) 维矩阵， Σ 2 , 1 = Σ 1 , 2 T \Sigma_{2,1}=\Sigma_{1,2}^T Σ2,1​=Σ1,2T​ ， Σ 2 , 2 \Sigma_{2,2} Σ2,2​ 为 ( d − d 1 ) ∗ ( d − d 1 ) (d-d_1)*(d-d_1) (d−d1​)∗(d−d1​) 维矩阵。

  可得
  p ( x 1 ∣ x 2 ) ∝ p ( x 1 , x 2 ) ∝ e x p [ − 0.5 ( x 1 − μ 1 ) T Σ c , 1 − 1 ( x 1 − μ 1 ) − 0.5 ( x 2 − μ 2 ) T Σ c , 2 − 1 ( x 2 − μ 2 ) − ( x 1 − μ 1 ) T Σ c , 1 − 1 Σ 1 , 2 Σ 2 , 2 − 1 ( x 2 − μ 2 ) ] p(x_1|x_2)\propto p(x_1,x_2)\\ \propto exp[-0.5(x_1-\mu_1)^T\Sigma_{c,1}^{-1}(x_1-\mu_1)-0.5(x_2-\mu_2)^T\Sigma_{c,2}^{-1}(x_2-\mu_2)-(x_1-\mu_1)^T\Sigma_{c,1}^{-1}\Sigma_{1,2}\Sigma_{2,2}^{-1}(x_2-\mu_2)] p(x1​∣x2​)∝p(x1​,x2​)∝exp[−0.5(x1​−μ1​)TΣc,1−1​(x1​−μ1​)−0.5(x2​−μ2​)TΣc,2−1​(x2​−μ2​)−(x1​−μ1​)TΣc,1−1​Σ1,2​Σ2,2−1​(x2​−μ2​)]
  其中
  Σ c , 1 = Σ 1 , 1 − Σ 1 , 2 Σ 2 , 2 − 1 Σ 2 , 1 Σ c , 2 = Σ 2 , 2 − Σ 2 , 1 Σ 1 , 1 − 1 Σ 1 , 2 \Sigma_{c,1}=\Sigma_{1,1}-\Sigma_{1,2}\Sigma_{2,2}^{-1}\Sigma_{2,1}\\\Sigma_{c,2}=\Sigma_{2,2}-\Sigma_{2,1}\Sigma_{1,1}^{-1}\Sigma_{1,2} Σc,1​=Σ1,1​−Σ1,2​Σ2,2−1​Σ2,1​Σc,2​=Σ2,2​−Σ2,1​Σ1,1−1​Σ1,2​
  忽略不随 x 1 x_1 x1​ 变化的相乘因子，得到
  p ( x 1 ∣ x 2 ) ∝ e x p [ − 0.5 ( x 1 − μ c , 1 ) T Σ c , 1 − 1 ( x 1 − μ c , 1 ) ] ∝ N ( x 1 ; μ c , 1 , Σ c , 1 ) p(x_1|x_2)\propto exp[-0.5(x_1-\mu_{c,1})^T\Sigma_{c,1}^{-1}(x_1-\mu_{c,1})]\propto N(x_1;\mu_{c,1},\Sigma_{c,1}) p(x1​∣x2​)∝exp[−0.5(x1​−μc,1​)TΣc,1−1​(x1​−μc,1​)]∝N(x1​;μc,1​,Σc,1​)

2. M e t r o p o l i s − H a s t i n g s Metropolis-Hastings Metropolis−Hastings 算法

• 在先验/后验分布中的支集中，为 θ \theta θ 任选一个初始值 θ 0 \theta^0 θ0。
• 在算法的第 t + 1 t+1 t+1 次迭代，进行以下操作：
• 从建议分布 q ( θ ∣ θ t ) q(\theta|\theta^t) q(θ∣θt) 中取一个样本 θ \theta θ。
• 从[0，1]区间的均匀分布中取一个样本 u u u。
• 利用下式计算： α ( θ t , θ ) = m i n ( 1 , g ( θ ) q ( θ t ∣ θ ) g ( θ t ) q ( θ ∣ θ t ) ) \alpha(\theta^t,\theta)=min(1,\frac{g(\theta)q(\theta^t|\theta)}{g(\theta^t)q(\theta|\theta^t)}) α(θt,θ)=min(1,g(θt)q(θ∣θt)g(θ)q(θt∣θ)​)
• 如果 u ⩽ α ( θ t ∣ θ ) u \leqslant \alpha(\theta^t|\theta) u⩽α(θt∣θ) ，则接受 θ \theta θ ，并令 θ t + 1 = θ t \theta^{t+1}=\theta^t θt+1=θt
• 经大量迭代所采集到的样本 θ t \theta^t θt ，就如同随机采样于密度函数 f ( θ ) f(\theta) f(θ)。

M e t r o p o l i s − H a s t i n g s Metropolis-Hastings Metropolis−Hastings 过程产生了不同于吉布斯抽样的马尔科夫链，但它们有着相同的分布限制 g ( θ ) / ∫ g ( θ ′ ) d θ ′ g(\theta)/\int g(\theta')d\theta' g(θ)/∫g(θ′)dθ′ 。建议分布可以采用任何实用的形式，但平稳分布依旧是 g ( θ ) / ∫ g ( θ ′ ) d θ ′ g(\theta)/\int g(\theta')d\theta' g(θ)/∫g(θ′)dθ′ 。然而，假设 q ( X ∣ Y ) q(X|Y) q(X∣Y) 是多元正态分布， Y Y Y是其均值， Σ \Sigma Σ 是其固定协方差矩阵。如果， Σ \Sigma Σ 的尺度过小，会导致高接受率，相容性变差。

还有一种单变量的 M e t r o p o l i s − H a s t i n g s Metropolis-Hastings Metropolis−Hastings 算法，如果样本总是被接受，则算法等同于吉布斯采样，常常被称为吉布斯的 M e t r o p o l i s Metropolis Metropolis 算法。

3. 可逆跳跃马尔科夫链蒙特卡罗方法

M C M C MCMC MCMC 方法可以拓展到模型阶次未知的情形，这时，使用可逆跳转马尔科夫链蒙特卡洛 ( R J M C M C ) (RJMCMC) (RJMCMC) 过程。正如标准的 M C M C MCMC MCMC 技术，模型维度不变时的参数，更新， R J M C M C RJMCMC RJMCMC 方法允许通过引入新参数或去掉现有参数，改变模型的维数。

假定打算从 k k k 维状态空间下的模型 M k M_k Mk​ 运动到 k ′ k' k′ 维状态空间下的模型 M k ′ M_{k'} Mk′​ ，其中 k ′ = k + d k'=k+d k′=k+d ， d d d 是大于等于1的整数。这个维数改变的方式是在当前状态 θ k \theta_k θk​ 之外抽取一个连续随机变量的 d d d 维向量 u d u_d ud​ ，再通过双映射 θ k ′ = θ ′ = ( θ k , u d ) \theta_{k'}=\theta'=(\theta_k,u_d) θk′​=θ′=(θk​,ud​) 定义一个新的状态 θ k ′ \theta_{k'} θk′​ 。通过指定运动接受概率也可使精细平衡满足
α ( θ k , θ k ′ ) = m i n ( 1 , p ( θ k ′ ∣ D ) q k ′ , k p ( θ k ∣ D ) q k , k ′ p ( u d ) ∣ ∂ θ k ′ ∂ ( θ k , u d ) ∣ ) \alpha(\theta_k,\theta_{k'})=min(1,\frac{p(\theta_{k'}|D)q_{k',k}}{p(\theta_k|D)q_{k,k'}p(u_d)}|\frac{\partial \theta_{k'}}{\partial(\theta_k,u_d)}|) α(θk​,θk′​)=min(1,p(θk​∣D)qk,k′​p(ud​)p(θk′​∣D)qk′,k​​∣∂(θk​,ud​)∂θk′​​∣)
其中， p ( θ ∣ D ) p(\theta|D) p(θ∣D) 表示数据 D D D 下参数 θ \theta θ 的后验分布； q k , k ′ , q k ′ , k q_{k,k'},q_{k',k} qk,k′​,qk′,k​ 分别是模型从 M k M_k Mk​ 运动到 M k ′ M_{k'} Mk′​ 以及模型从 M k ‘ M_{k‘} Mk‘​ 运动到 M k M_k Mk​ 的选择概率； p ( u d ) p(u_d) p(ud​) 表示随机变量 u d u_d ud​ 的概率密度函数； ∣ ∂ θ k ′ ∂ ( θ k , u d ) ∣ |\frac{\partial \theta_{k'}}{\partial(\theta_k,u_d)}| ∣∂(θk​,ud​)∂θk′​​∣ 是变量 ( θ k , u d ) (\theta_k,u_d) (θk​,ud​) 向 θ k ′ \theta_{k'} θk′​ 变换的雅可比行列式。

可以抽象为两个运动模型：生长运动和消亡运动。由于 k ′ = k + d k'=k+d k′=k+d ， d d d 是大于等于1的整数，所以，生长运动的选择概率为 q k , k ′ q_{k,k'} qk,k′​，消亡运动的选择概率为 q k ′ , k q_{k',k} qk′,k​。

4. 切片采样

切片采样是在概率密度函数曲线或与概率密度函数成比例的函数曲线下的区域进行移植采样来获取分布的样本，是用到增广数据的吉布斯采样的一个特例，这里就不再赘述。