统计计算五|MCMC( Markov Chain Monte Carlo)

系列文章目录

统计计算一|非线性方程的求解
统计计算二|EM算法(Expectation-Maximization Algorithm,期望最大化算法)
统计计算三|Cases for EM
统计计算四|蒙特卡罗方法(Monte Carlo Method)

文章目录

  • 系列文章目录
  • 一、基本概念
    • (一)马尔科夫链
      • 1、定义
      • 2、性质
      • 3、常返
      • 4、平稳分布
    • (二)MCMC原理
      • 1、核心思想
      • 2、连续状态
      • 3、MCMC估计期望的步骤
  • 二、满条件分布
      • 1、定义
      • 2、考虑MCMC中的应用
      • 3、伽玛分布
  • 三、Metropolis–Hastings 算法
  • 四、Gibbs 算法
    • (一)基本概念
    • (二)采样过程
    • (三)延展
    • (四)示例
      • 1、简单正态例子:
      • 2、多项分布示例:
      • 3、分类消费者例子
  • 五、实施
    • (一)混合和收敛
    • (二)相关术语


一、基本概念

(一)马尔科夫链

1、定义

考虑在每个时间段有一个值的随机过程,令 X n X_n Xn表示它在时间段 n n n的值,如果:
P { X n + 1 = j ∣ X n = i , X n − 1 = i n − 1 , . . . , X 1 = i 1 } = P { X n + 1 = j ∣ X n = i } = P i j \begin{aligned} &P\{X_{n+1}=j|X_{n}=i,X_{n-1}=i_{n-1},...,X_1=i_1\} \\ =&P\{X_{n+1}=j|X_{n}=i \} \\ =&P_{ij} \end{aligned} ==P{Xn+1=jXn=i,Xn1=in1,...,X1=i1}P{Xn+1=jXn=i}Pij
这样的随机过程称为马尔科夫链。 P P P为一步转移概率 P i j P_{ij} Pij的矩阵,n次转移后的矩阵为 P n P^n Pn

对随机变量序列 { X 0 , X 1 , X 2 , . . . } \{X_0,X_1,X_2,...\} {X0,X1,X2,...},在任一时刻 n ≥ 0 n\geq 0 n0,序列中下一时刻 n + 1 n+1 n+1 X n + 1 X_{n+1} Xn+1由条件分布 P ( x ∣ X n ) P(x|X_n) P(xXn)产生,它只依赖于时刻 n n n的当前状态,而与时刻 n n n以前的历史状态 { X 0 , . . . , X n − 1 } \{X_0,...,X_{n-1}\} {X0,...,Xn1}无关。满足这样条件的随机变量序列称为Markov链。

若转移概率不随n改变,则称此链为时间齐性的,否则为时间非齐性。

2、性质

  • 不可约:如果从任一状态 i i i 经有限步后都可到达任一状态 j j j,即对于任两个状态 i i i, j j j,都存在 m > 0 m > 0 m>0 使得 p ( X ( m + n ) = j ∣ X ( n ) = i ) > 0 p(X^{(m+n)} = j|X^{(n)} = i) > 0 p(X(m+n)=jX(n)=i)>0(联通的),则称这一条马氏链是不可约的

  • 常返:称能以概率1回来的状态是常返的:
    f i i = ∑ i = 1 ∞ f i i ( n ) = 1 f_{ii}=\sum_{i=1}^∞f_{ii}^{(n)}=1 fii=i=1fii(n)=1
    其中 f i j ( n ) = p ( X ( n ) = j , X k ≠ j , k = 1 , 2 , . . . , n − 1 ∣ X 0 = i ) f_{ij}^{(n)}=p(X^{(n)}=j,X_k\neq j,k=1,2,...,n-1|X_0=i) fij(n)=p(X(n)=j,Xk=j,k=1,2,...,n1∣X0=i)为马氏链在0时从状态 i i i 出发,经过n步转移后,首次到达状态 j j j 的概率。若返回概率小于1,即 f i i < 1 f_{ii}<1 fii<1,则为非常返态。

  • 正常返和零常返:一个状态的平均返回时间 μ i \mu_i μi
    μ i = ∑ n = 1 ∞ n f i i ( n ) \mu_i=\sum_{n=1}^∞nf_{ii}^{(n)} μi=n=1nfii(n)
    μ i < ∞ \mu_i<∞ μi<,则为正常返,反之 μ i = ∞ \mu_i=∞ μi=为零常返。如果状态空间有限,其常返状态都是非零常返

  • 周期性:

    • 如果马尔可夫链只能以一定的规则间隔访问状态空间的某些部分,则它是周期性的。如果由状态 j j j 经非 d d d 整数倍步到达 j j j 的概率为 0,则称状态 j j j 具有周期 d d d,周期可定义为集合 { n : n ≥ 0 , p i i ( n ) > 0 } \{n : n ≥ 0, p^{(n)}_{ii} > 0\} {n:n0,pii(n)>0} 的最大公约数,即所有返回可能的次数的最大公约数。
    • 如果一条马氏链的每一个状态的周期都为 1, 则称此链为非周期的。
  • 遍历的:如果一条马氏链是不可约、非周期,且其所有状态都是非零常返的,则称之为遍历的
    在这里插入图片描述

3、常返

i i i 是常返态,则:

  • i i i 是零常返态的充要条件为: l i m n → ∞ p i i ( n ) = 0 lim_{n\rightarrow ∞}p_{ii}^{(n)}=0 limnpii(n)=0
  • i i i 是正常返态的充要条件为: l i m n → ∞ p i i ( n ) > 0 lim_{n\rightarrow ∞}p_{ii}^{(n)}>0 limnpii(n)>0
    在这里插入图片描述

4、平稳分布

平稳分布定义:若一离散分布 π π π 满足 π T P = π T π^TP = π^T πTP=πT,则称之为转移概率矩阵为 P P P 的马氏链的平稳分布。

  • (充分条件)如果一条时间齐性的马氏链满足:
    π i p i j = π j p j i , ∀ i , j ∈ S , \pi_ip_{ij}=\pi_jp_{ji},∀i, j ∈ S, πipij=πjpji,i,jS,
    π \pi π是此链的平稳分布,且称此链为可逆的,上述方程也称为细致平衡。

如果一个转移概率阵为 P P P,平稳分布为 π \pi π的马氏链是不可约的,正常返的,且非周期的(遍历),则 π \pi π唯一,且满足:
lim ⁡ n → ∞ P [ X n = j ] = π j \lim_{n\rightarrow ∞}P[X_n=j]=\pi_j nlimP[Xn=j]=πj
也就是极限分布即为平稳分布。其中 π j \pi_j πj π \pi π的第 j j j个元素,且满足如下方程组:
π j ≥ 0 , ∑ i ∈ S π i = 1 , 且 π j = ∑ i ∈ S π i p i j , ∀ j ∈ S \pi_j\geq 0,\sum_{i\in S}\pi_i=1,且\pi_j=\sum_{i\in S}\pi_ip_{ij},∀j ∈ S πj0,iSπi=1,πj=iSπipij,jS

推论:如果 X 1 , X 2 , . . . X_1,X_2,... X1,X2,...是一不可约、正常返的、非周期的平稳分布为 π \pi π的马氏链值,则 X ( n ) X^{(n)} X(n)依分布收敛到分布为 π \pi π的随机变量。

(二)MCMC原理

1、核心思想

MCMC 的核心是构建转移矩阵,使得我们的目标分布 (π) 满足细致平衡,也即目标分布是构建出的马尔科夫链的平稳分布。又因为遍历性,这条链的极限分布就是唯一的平稳分布。所以我们只要迭代次数足够大,就能假设达到了平稳分布,Xn 即为目标分布的样本。

2、连续状态

在连续分布情况下,对任一可测集 B B B,一步转移概率定义为:
P ( x → B ) = ∫ B p ( x , x ′ ) d x ′ P(x\rightarrow B)=\int_Bp(x,x')dx' P(xB)=Bp(x,x)dx
转移概率 p ( ⋅ , ⋅ ) p(·, ·) p(⋅,⋅)称为Markov链转移核。通常假定 p ( ⋅ , ⋅ ) p(·, ·) p(⋅,⋅) t t t无关,即基于该转移核的Markov链是时间齐次的。

例:根据转移核 P ( X ( t + 1 ) ∣ X ( t ) ) ∼ N ( 0.5 X ( t ) , 1 ) P(X^{(t+1)}|X^{(t)})\sim N(0.5X^{(t)},1) P(X(t+1)X(t))N(0.5X(t),1),产生平稳分布是$N(0,4/3)的Markov链。
在这里插入图片描述

3、MCMC估计期望的步骤

MCMC方法估计 E π f ( X ) = ∫ f ( x ) π ( x ) d x E_\pi f(X)=\int f(x)\pi(x)dx Eπf(X)=f(x)π(x)dx的基本流程:

  • 选择转移核 p ( ⋅ , ⋅ ) p(·, ·) p(⋅,⋅)(参数更新公式),使其平稳分布是 π ( x ) \pi(x) π(x)
  • 从某一点 X ( 0 ) X^{(0)} X(0)出发,用上述转移核 p ( ⋅ , ⋅ ) p(·, ·) p(⋅,⋅)产生Markov链序列 X ( 0 ) , X ( 1 ) , . . . X ( n ) X^{(0)},X^{(1)},...X^{(n)} X(0),X(1),...X(n)
  • 对较大的 n n n,选择合适的 m m m E π f ( X ) = ∫ f ( x ) π ( x ) d x E_\pi f(X)=\int f(x)\pi(x)dx Eπf(X)=f(x)π(x)dx的估计为:
    E ^ π f = 1 n − m ∑ t = m + 1 n f ( X ( t ) ) \hat{E}_\pi f=\frac{1}{n-m}\sum_{t=m+1}^nf(X^{(t)}) E^πf=nm1t=m+1nf(X(t))

上述步骤中构造的转移核 p ( ⋅ , ⋅ ) p(·, ·) p(⋅,⋅) 使得概率分布 π ( x ) π(x) π(x) 为其平稳分布最为重要。如何构造合适的转移核是 MCMC 方法主要研究的问题,不同的 MCMC 方法主要区别就是转移核的构造方法不同。

二、满条件分布

1、定义

MCMC方法中转移核 p ( x , x ′ ) p(x,x') p(x,x)的构造大多建立在形如 π ( x T ∣ x − T ) \pi(x_T|x_{-T}) π(xTxT)的条件分布上:

  • x T = { x i , i ∈ T } x_T=\{x_i,i\in T\} xT={xi,iT} x − T = { x i , i ∉ T } x_{-T}=\{x_i,i \notin T\} xT={xi,i/T} T ⊂ I = { 1 , ⋅ ⋅ ⋅ , p } T⊂ I = \{1, · · · , p\} TI={1,⋅⋅⋅,p}
  • 在条件分布 π ( x T ∣ x − T ) \pi(x_T|x_{-T}) π(xTxT)中,p个所有的变量或者出现在条件中,或者出现在变元中,这种条件分布称为满条件分布
    π ( x T ∣ x − T ) = π ( x ) ∫ π ( x ) d x T ∝ π ( x ) \pi(x_T|x_{-T})=\frac{\pi(x)}{\int \pi(x)dx_T}∝ π(x) π(xTxT)=π(x)dxTπ(x)π(x)
    进而对任意 x , x ′ ∈ X x,x'\in X x,xX x − T = x − T ′ x_{-T}=x'_{-T} xT=xT,有:
    π ( x T ′ ∣ x − T ′ ) π ( x T ∣ x − T ) = π ( x ′ ) π ( x ) \frac{\pi(x'_T|x'_{-T})}{\pi(x_T|x_{-T})}=\frac{\pi(x')}{π(x)} π(xTxT)π(xTxT)=π(x)π(x)

2、考虑MCMC中的应用

一般情况下 y = ( x , z , θ ) y=(x,z,\theta) y=(x,z,θ),其中 x x x表示观测数据, z z z表示缺失数据, θ \theta θ表示参数。令 p ( x , z ∣ θ ) p(x,z|\theta) p(x,zθ)表示完全数据的密度函数, π ( θ ) \pi(\theta) π(θ)表示 θ \theta θ的先验分布, f ( y ) = p ( x , z ∣ θ ) π ( θ ) f(y)=p(x,z|\theta)\pi(\theta) f(y)=p(x,zθ)π(θ),则 y y y的满条件分布为:
f ( z i ∣ z − i , x , θ ) ∝ p ( x , z ∣ θ ) π ( θ i ∣ θ − i , x , z ) ∝ f ( y ) ∝ p ( x , z ∣ θ ) π ( θ i ∣ θ − i ) f ( x i ∣ x − i , z , θ ) ∝ p ( x , z ∣ θ ) f(z_i|z_{-i},x,\theta) ∝ p(x, z|θ)\\ π(θ_i|θ_{−i}, x, z) ∝ f(y) ∝ p(x, z|θ)π(θ_i|θ_{−i})\\ f(x_i|x_{−i}, z, θ) ∝ p(x, z|θ) f(zizi,x,θ)p(x,zθ)π(θiθi,x,z)f(y)p(x,zθ)π(θiθi)f(xixi,z,θ)p(x,zθ)
其中 θ − i = { θ j : j ≠ i } , z − i = { z j : j ≠ i } , x − i = { x j : j ≠ i } \theta_{-i}=\{\theta_j:j\neq i\},z_{-i}=\{z_j:j\neq i\},x_{-i}=\{x_j:j\neq i\} θi={θj:j=i},zi={zj:j=i},xi={xj:j=i}

例:设 x = ( x 1 , x 2 ) x=(x_1,x_2) x=(x1,x2)的密度函数为:
π ( x 1 , x 2 ) ∝ e x p { − 1 2 ( x 1 − 1 ) 2 ( x 2 − 1 ) 2 } \pi(x_1,x_2) ∝ exp\{-\frac{1}{2}(x_1-1)^2(x_2-1)^2\} π(x1,x2)exp{21(x11)2(x21)2}
则满条件分布为:

  • π ( x 1 ∣ x 2 ) ∝ π ( x 1 , x 2 ) ∝ e x p { − 1 2 ( x 1 − 1 ) 2 ( x 2 − 1 ) 2 } = N ( 1 , ( x 2 − 1 ) − 2 ) π (x_1 | x_2) ∝ π (x_1, x_2)∝ exp\{−\frac{1}{2}(x_1 − 1)^2(x_2-1)^2\} = N(1,(x_2 − 1)^{−2}) π(x1x2)π(x1,x2)exp{21(x11)2(x21)2}=N(1,(x21)2)
  • π ( x 2 ∣ x 1 ) ∝ π ( x 1 , x 2 ) ∝ e x p { − 1 2 ( x 2 − 1 ) 2 ( x 1 − 1 ) 2 } = N ( 1 , ( x 1 − 1 ) − 2 ) π (x_2 | x_1) ∝ π (x_1, x_2)∝ exp\{−\frac{1}{2}(x_2 − 1)^2(x_1-1)^2\} = N(1,(x_1 − 1)^{−2}) π(x2x1)π(x1,x2)exp{21(x21)2(x11)2}=N(1,(x11)2)

3、伽玛分布

伽玛分布 X ∼ Γ ( α , λ ) X\sim \Gamma(\alpha,\lambda) XΓ(α,λ)的密度函数为:
f ( x ) = x ( α − 1 ) λ α e ( − λ x ) Γ ( α ) , x > 0 f(x)=\frac{x^{(\alpha-1)}\lambda^\alpha e^{(-\lambda x)}}{\Gamma(\alpha)},x>0 f(x)=Γ(α)x(α1)λαe(λx),x>0
其中Gamma函数之特征为:
{ Γ ( α ) = ( α − 1 ) ! , if  α is  Z + Γ ( α ) = ( α − 1 ) Γ ( α − 1 ) , if  α is  R + Γ ( 1 2 ) = π \begin{cases} \Gamma(\alpha)=(\alpha-1)!, & \text{if $\alpha$ is $Z^+$} \\ \Gamma(\alpha)=(\alpha-1)\Gamma(\alpha-1), & \text{if $\alpha$ is $R^+$} \\ \Gamma(\frac{1}{2})=\sqrt{\pi}\\ \end{cases} Γ(α)=(α1)!,Γ(α)=(α1)Γ(α1),Γ(21)=π if α is Z+if α is R+
指数分布为 α = 1 \alpha=1 α=1的伽马分布

例:设 y 1 , . . . , y n y_1,..., y_n y1,...,yn 独立同分布,来自正态分布 N ( µ , τ − 1 ) N(µ, τ^{−1}) N(µ,τ1)。 参数 µ , τ − 1 µ, τ^{−1} µ,τ1的先验分布分别为正态分布 µ ∼ N ( 0 , 1 ) µ ∼ N(0, 1) µN(0,1), 伽马分布 τ ∼ Γ ( 2 , 1 ) τ ∼ Γ(2, 1) τΓ(2,1), 且 µ µ µ τ τ τ 独立,计算满条件分布。
在这里插入图片描述

满条件分布并不都能表示为显示形式。

例:样本 y i , i = 1 , . . . , n y_i,i=1,...,n yi,i=1,...,n独立且 y i ∼ B ( 1 , p i ) y_i\sim B(1,p_i) yiB(1,pi)(伯努利分布),其中 p i = ( 1 + e x p ( − ( α + β x i ) ) ) − 1 p_i=(1+exp(-(\alpha+\beta x_i)))^{-1} pi=(1+exp((α+βxi)))1 x i x_i xi假设是固定的,已知的。参数 α , β \alpha,\beta α,β的先验分布分别为 α ∼ N ( 0 , 1 ) \alpha \sim N(0,1) αN(0,1),且 α , β \alpha,\beta α,β独立。
在这里插入图片描述

三、Metropolis–Hastings 算法

(一)基本概念

1、思想

Metropolis-Hastings 算法是马尔可夫链蒙特卡洛 (MCMC) 方法的一种,用于从难直接抽样的概率分布中获取一系列随机样本。它通过构建一个Markov 链来实现,使其平衡分布等于所需分布。

Metropolis-Hastings 方法转移核的构造:
p ( x , x ′ ) = q ( x ′ ∣ x ) α ( x , x ′ ) p(x,x')=q(x'|x)\alpha(x,x') p(x,x)=q(xx)α(x,x)
潜在的转移核 q ( x ′ ∣ x ) q(x'|x) q(xx)作为 x ′ x' x的函数是一个概率密度或概率分布,被称为提案分布。提案分布可以取各种形式,常把它取为易于产生随机数的分布。

2、具体实施

  • 如果链在时刻 t t t处于状态 x x x,即 X t = x X_t=x Xt=x
  • q ( ⋅ ∣ x ) q(· | x) q(x)产生一个潜在的转移 x → x ′ x\rightarrow x' xx
  • 根据概率 α ( x , x ′ ) \alpha(x,x') α(x,x)决定是否转移

α ( x , x ′ ) \alpha(x,x') α(x,x)是接受概率,满足 0 < α ( x , x ′ ) ≤ 1 0<\alpha(x,x')\leq 1 0<α(x,x)1。也就是说,在潜在转移点 x ′ x' x找到后,以概率 α ( x , x ′ ) \alpha(x,x') α(x,x)接受 x ′ x' x作为链在下一时刻 t + 1 t+1 t+1的状态值,而以概率 1 − α ( x , x ′ ) 1-\alpha(x,x') 1α(x,x)拒绝转移到 x ′ x' x,从而链在下一时刻 t + 1 t+1 t+1仍处于状态 x x x。在实际计算中,产生区间 [ 0 , 1 ] [0,1] [0,1]上均匀分布的随机数 u u u,令:
X t + 1 = { x ′ u ≤ α ( x , x ′ ) x u > α ( x , x ′ ) X_{t+1} = \begin{cases} x' & \text{$u\leq \alpha(x,x')$} \\ x & \text{$u>\alpha(x,x')$} \end{cases} Xt+1={xxuα(x,x)u>α(x,x)

3、算法过程

假设 π ( x ) \pi(x) π(x)为目标概率分布,MH算法的过程为:

  • 初始化:选定初始状态 x 0 x_0 x0,令 t = 0 t=0 t=0
  • 迭代过程:
    • 生成:从某一容易抽样的分布 q ( x ′ ∣ x t ) q(x'|x_t) q(xxt)中随机生成候选状态 x ′ x' x
    • 计算:计算是否采纳候选状态的概率
      α ( x t , x ′ ) = m i n ( 1 , π ( x ′ ) π ( x t ) q ( x t ∣ x ′ ) q ( x ′ ∣ x t ) ) \alpha(x_t,x')=min(1,\frac{\pi(x')}{\pi(x_t)}\frac{q(x_t|x')}{q(x'|x_t)}) α(xt,x)=min(1,π(xt)π(x)q(xxt)q(xtx))
    • 接受或拒绝
      • [ 0 , 1 ] [0,1] [0,1]的均匀分布中生成随机数 u u u
      • 如果 u ≤ α ( x t , x ′ ) u\leq \alpha(x_t,x') uα(xt,x),则接受该状态,并令 x t + 1 = x ′ x_{t+1}=x' xt+1=x
      • 如果 u > α ( x t , x ′ ) u>\alpha(x_t,x') u>α(xt,x),则拒绝该状态,并令 x t + 1 = x t x_{t+1}=x_t xt+1=xt
    • 增量:令 t = t + 1 t=t+1 t=t+1

MH算法的推导:
在这里插入图片描述

(二)Metropolis 选择(提案分布 q ( x , x ′ ) q(x, x′) q(x,x) 的选取)

  • Metropolis 建议 q ( x , x ′ ) q(x, x′) q(x,x) 为对称分布,即:
    q ( x , x ′ ) = q ( x ′ , x ) , ∀ x , x ′ q(x,x')=q(x',x),∀x, x′ q(x,x)=q(x,x),x,x
    • α ( x , x ′ ) \alpha(x,x') α(x,x)简化为:
      α ( x , x ′ ) = min ⁡ { 1 , π ( x ′ ) π ( x t ) q ( x ′ , x ) q ( x , x ′ ) } = min ⁡ { 1 , π ( x ′ ) π ( x ) } \alpha(x,x')=\min\{1,\frac{\pi(x')}{\pi(x_t)}\frac{q(x',x)}{q(x,x')}\}=\min\{1,\frac{\pi(x')}{\pi(x)}\} α(x,x)=min{1,π(xt)π(x)q(x,x)q(x,x)}=min{1,π(x)π(x)}
    • 常用的对称分布形式为: q ( x , x ′ ) = f ( ∣ x − x ′ ∣ ) q(x,x')=f(|x-x'|) q(x,x)=f(xx),被称为随机移动Metropolis算法,具体例子为 q ( x , x ′ ) ∝ e x p { − ( x ′ − x ) 2 / e } q(x,x')∝exp\{-(x'-x)^2/e\} q(x,x)exp{(xx)2/e}

例:生成一个Markov链,使得其平稳分布为柯西分布,
f ( x ) = 1 π 1 1 + x 2 f(x)=\frac{1}{\pi}\frac{1}{1+x^2} f(x)=π11+x21
选定的提案分布 q ( x , x ′ ) q(x,x') q(x,x) N ( x , b 2 ) N(x,b^2) N(x,b2),其中 b b b为任意常数,如0.1,1,10,此时:
α ( x , x ′ ) = m i n { 1 , π ( x ′ ) π ( x ) } = m i n { 1 , 1 + x 2 1 + ( x ′ ) 2 } \alpha(x,x')=min\{1,\frac{\pi(x')}{\pi(x)}\}=min\{1,\frac{1+x^2}{1+(x')^2}\} α(x,x)=min{1,π(x)π(x)}=min{1,1+(x)21+x2}

  • 独立抽样:如果 q ( x , x ′ ) q(x,x') q(x,x)与当前状态 x x x无关,即 q ( x , x ′ ) = q ( x ′ ) q(x,x')=q(x') q(x,x)=q(x),则由此提案分布导出的MH算法称为独立抽样。此时的 α ( x , x ′ ) \alpha(x,x') α(x,x)为:
    α ( x , x ′ ) = min ⁡ { 1 , π ( x ′ ) π ( x t ) q ( x ′ ) q ( x ) } \alpha(x,x')=\min\{1,\frac{\pi(x')}{\pi(x_t)}\frac{q(x')}{q(x)}\} α(x,x)=min{1,π(xt)π(x)q(x)q(x)}
    如果 q ( x ) q(x) q(x)接近 π ( x ) \pi(x) π(x),基于独立抽样获取的Markov链的收敛效果更好。

给定数据 Y 1 , . . . , Y n ∼ i . i . d N ( θ , 1 ) Y_1,... , Y_n\stackrel{i.i.d}{\sim}N(θ, 1) Y1,...,Yni.i.dN(θ,1), 先验分布 π ( θ ) = 1 / π ( 1 + θ 2 ) π(θ) = 1/{π(1 + θ^2)} π(θ)=1/π(1+θ2),此时后验分布为
π ( θ ∣ Y 1 , . . . , Y n ) ∝ p ( Y 1 , . . . , Y n ∣ θ ) π ( θ ) ∝ e x p { − n ( θ − y ˉ ) 2 / 2 } × 1 1 + θ 2 \begin{aligned} \pi(\theta|Y_1,...,Y_n)&∝ p(Y_1, . . . , Y_n | θ)π(θ)\\ &∝ exp\{-n(\theta-\bar{y})^2/2\}\times\frac{1}{1+\theta^2} \end{aligned} π(θY1,...,Yn)p(Y1,...,Ynθ)π(θ)exp{n(θyˉ)2/2}×1+θ21
假设已有数据给定 n = 40 , y ˉ = 0.14 n=40,\bar{y}=0.14 n=40,yˉ=0.14,此时使用 x x x x ′ x' x的记号, θ \theta θ的后验分布为:
π ( θ ) ∝ e x p { − 40 ( x − 0.14 ) 2 / 2 } × 1 1 + x 2 \pi(\theta)∝ exp\{-40(x-0.14)^2/2\}\times\frac{1}{1+x^2} π(θ)exp{40(x0.14)2/2}×1+x21
求后验期望 E ( x ) E(x) E(x)
在这里插入图片描述

(三)单元素 Metropolis-Hastings 算法

在 X 是 p 维的情况,同时产生整个 X 有时是困难的 (接受率特别低),而将 X 根据其分量逐个进行抽样则简单得多,这就要用到条件分布,特别是满条件分布性质。

单元素 Metropolis-Hastings 算法的核心思想:对于 p p p 维变量 X X X, 基于 p − 1 p − 1 p1维变量 X − i X_{−i } Xi的条件分布 X i ∣ X − i , i = 1 , . . . , p X_i| X_{−i}, i = 1, ... , p XiXi,i=1,...,p,选择转移核 q i ( x i → x i ′ ∣ X − i = x − i ) q_i (x_i → x^′_i|X_{-i}=x_{-i}) qi(xixiXi=xi)。由转移核 q i ( x i → x i ′ ∣ X − i = x − i ) q_i (x_i → x^′_i|X_{-i}=x_{-i}) qi(xixiXi=xi) 产生可能的 x i ′ x^′_i xi, 以概率
α i ( x i → x i ′ ∣ x − i ) = m i n ( 1 , π ( x ′ ) q i ( x i ′ → x i ∣ x − i ) π ( x ) 1 i ( x i → x i ′ ∣ x − i ) ) \alpha_i(x_i\rightarrow x^′_i|x_{-i})=min(1,\frac{\pi(x')q_i( x^′_i\rightarrow x_i|x_{-i})}{\pi(x)1_i(x_i\rightarrow x^′_i|x_{-i})}) αi(xixixi)=min(1,π(x)1i(xixixi)π(x)qi(xixixi))
决定是否接受 x ′ x' x为链的下一状态。

即每次只更新一个元素,其他保持不变:
在这里插入图片描述

四、Gibbs 算法

(一)基本概念

Gibbs 抽样是一种单元素 Metropolis-Hastings 算法的特殊情况,单元素Metropolis-Hastings 算法中取 q i ( x i → x i ′ ∣ X − i = x − i ) q_i (x_i → x^′_i|X_{-i}=x_{-i}) qi(xixiXi=xi) π ( x i ∣ x − i ) \pi(x_i|x_{-i}) π(xixi)。(也就是说直接定义为其满条件分布)
在这里插入图片描述此时的接受率为1,这意味着,我们不需要舍弃样本,每个更新后的值都为样本。

(二)采样过程

Gibbs采样的过程:

  • 确定初始值 X ( 1 ) X^{(1)} X(1)
  • 假设已得到样本 X ( i ) X^{(i)} X(i),记下一个样本为 X ( i + 1 ) = ( x 1 ( i + 1 ) , x 2 ( i + 1 ) , . . . , x n ( i + 1 ) ) X^{(i+1)}=(x_1^{(i+1)},x_2^{(i+1)},...,x_n^{(i+1)}) X(i+1)=(x1(i+1),x2(i+1),...,xn(i+1)),对其中某一分量 x j ( i + 1 ) x_j^{(i+1)} xj(i+1)可通过在其他分量已知的条件下该分量的概率分布来抽取该分量。对于此条件概率,我们使用样本 X ( i + 1 ) X^{(i+1)} X(i+1)中已得到的分量 x 1 ( i + 1 ) x_1^{(i+1)} x1(i+1) x j − 1 ( i + 1 ) x_{j-1}^{(i+1)} xj1(i+1)以及上一样本 X ( i ) X^{(i)} X(i)中的分量 x j + 1 ( i ) x_{j+1}^{(i)} xj+1(i) x n ( i ) x_n^{(i)} xn(i),即:
    f ( x j ( i + 1 ) ∣ x 1 ( i + 1 ) , . . . , x j − 1 ( i + 1 ) , x j + 1 ( i ) , . . . , x n ( i ) ) f(x_j^{(i+1)}|x_1^{(i+1)},...,x_{j-1}^{(i+1)},x_{j+1}^{(i)},...,x_n^{(i)}) f(xj(i+1)x1(i+1),...,xj1(i+1),xj+1(i),...,xn(i))
  • 重复上述过程
    在这里插入图片描述

(三)延展

  • 更新顺序:
    X 元素的更新顺序对于不同的循环是可以变化的. 有时候对每个循环而言, 使用随机顺序是比较合理的. 这被称作为随机扫描 Gibbs 抽样. 事实上, 甚至没有必要对每个循环中的每个元素都进行更新, 而只要每个元素的更新足够地频繁就可以了.
  • 区组化:
    当 X 的元素相关时, 区组化特别有用, 用其构造的算法能够使更相关的元素在同一个区组中被一起抽样出来.
    在这里插入图片描述
  • 混合Gibbs:在适当的时候使用不同的MH采样
    • 用某Gibbs迭代更新 X 1 ( t + 1 ) ∣ ( x 2 ( t ) , x 3 ( t ) , x 4 ( t ) , x 5 ( t ) , x 6 ( t ) ) X_1^{(t+1)}|(x_2^{(t)},x_3^{(t)},x_4^{(t)},x_5^{(t)},x_6^{(t)}) X1(t+1)(x2(t),x3(t),x4(t),x5(t),x6(t))
    • 用某Metropolis迭代更新 ( X 2 ( t + 1 ) , X 3 ( t + 1 ) ) ∣ ( x 1 ( t + 1 ) , x 4 ( t ) , x 5 ( t ) , x 6 ( t ) ) (X_2^{(t+1)},X_3^{(t+1)})|(x_1^{(t+1)},x_4^{(t)},x_5^{(t)},x_6^{(t)}) (X2(t+1),X3(t+1))(x1(t+1),x4(t),x5(t),x6(t))
    • 用某Metropolis迭代更新 X 4 ( t + 1 ) ∣ x 1 ( t + 1 ) , x 2 ( t + 1 ) ∣ ( x 3 ( t + 1 ) , x 5 ( t ) , x 6 ( t ) ) X_4^{(t+1)}|x_1^{(t+1)},x_2^{(t+1)}|(x_3^{(t+1)},x_5^{(t)},x_6^{(t)}) X4(t+1)x1(t+1),x2(t+1)(x3(t+1),x5(t),x6(t))
    • 用某Gibbs迭代更新 ( X 5 ( t + 1 ) , X 6 ( t + 1 ) ) ∣ ( x 1 ( t + 1 ) , x 2 ( t + 1 ) , x 3 ( t + 1 ) , x 4 ( t + 1 ) ) (X_5^{(t+1)},X_6^{(t+1)})|(x_1^{(t+1)},x_2^{(t+1)},x_3^{(t+1)},x_4^{(t+1)}) (X5(t+1),X6(t+1))(x1(t+1),x2(t+1),x3(t+1),x4(t+1))

当 X 的一个或者多个元素的一元边际密度没有显示表达的时候,Gibbs算法中的 Metropolis-Hastings 迭代特别有用. 有时也是 Gibbs 跳出局部最优的好方法

Gibbs虽然是寻找参数的后验分布,但是其本质为一个采样算法

(四)示例

1、简单正态例子:

(已知参数,但假设我们只能单独产生单变量正态分布的随机数,如何从二维正态采样)给定一个目标分布:
X = [ X 1 X 2 ] ∼ N ( [ u 1 μ 2 ] , [ σ 1 2 σ 12 σ 12 σ 2 2 ] ) X=\begin{bmatrix} X_1 \\ X_2 \\ \end{bmatrix}\sim N(\begin{bmatrix} u_1 \\ \mu_2 \\ \end{bmatrix},\begin{bmatrix} \sigma_1^2 & \sigma_{12} \\ \sigma_{12} & \sigma_2^2 \\ \end{bmatrix}) X=[X1X2]N([u1μ2],[σ12σ12σ12σ22])
在这里插入图片描述

2、多项分布示例:

p { X 1 = m 1 , X 2 = m 2 , . . . , X n = m n } = N ! m 1 ! m 2 ! . . . m n ! p 1 m 1 p x m 2 . . . p n m n p\{ X_1=m_1,X_2=m_2,...,X_n=m_n\}=\frac{N!}{m_1!m_2!...m_n!}p_1^{m1}p_x^{m2}...p_n^{m_n} p{X1=m1,X2=m2,...,Xn=mn}=m1!m2!...mn!N!p1m1pxm2...pnmn
其中 N = ∑ i = 1 n m i N=\sum_{i=1}^nm_i N=i=1nmi,某实验服从上述多项分布, N = 22 , n = 7 N=22,n=7 N=22,n=7,7个结果出现的概率分别为:
p : = ( p 1 , p 2 , . . . , p 7 ) = ( θ 4 , 1 8 , θ 4 , 3 8 , 1 2 ( 1 − θ − η ) ) p:=(p_1,p_2,...,p_7)=(\frac{\theta}{4},\frac{1}{8},\frac{\theta}{4},\frac{3}{8},\frac{1}{2}(1-\theta-\eta)) p:=(p1,p2,...,p7)=(4θ,81,4θ,83,21(1θη))
现有观测数据为 y = ( y 1 , y 2 , y 3 , y 4 , y 5 ) = ( 14 , 1 , 1 , 1 , 5 ) y=(y_1,y_2,y_3,y_4,y_5)=(14,1,1,1,5) y=(y1,y2,y3,y4,y5)=(14,1,1,1,5),且
( z 1 , y 1 − z 1 , y 2 , y 3 , z 2 , y 4 − z 2 , y 5 ) ∼ M ( 22 ; p ) (z_1,y_1-z_1,y_2,y_3,z_2,y_4-z_2,y_5)\sim M(22;p) (z1,y1z1,y2,y3,z2,y4z2,y5)M(22;p)
其中 M M M表示多项分布
在这里插入图片描述

3、分类消费者例子

数据: X g j X_{gj} Xgj 代表第 j j j 个消费者在上个月,一共从第 g g g 个类别中买了 X g j X_{gj} Xgj个物品。
问题:怎么通过消费偏好,将消费者归为 K K K 类?
在这里插入图片描述 c j c_j cj 代表第 j j j 个消费者属于的类别,取值范围为 1 , 2 , . . . , K 1, 2, ..., K 1,2,...,K。简单来看,其实就是有 n n n 个样本, G G G 个特征,目标是将这些样本分成 K组。

  • 缺失数据 c j c_j cj的边际分布: p ( c j = k ) = π k , ∑ k = 1 K π k = 1 p(c_j=k)=\pi_k,\sum_{k=1}^K\pi_k=1 p(cj=k)=πk,k=1Kπk=1
  • λ g k \lambda_{gk} λgk表示第k类消费者买第g类物品的均值: X g j ∣ c j = k ∼ P o i s ( λ g k ) X_{gj}|c_j=k\sim Pois(\lambda_{gk}) Xgjcj=kPois(λgk)

在贝叶斯统计中,共轭先验分布是指这样一种先验分布:当它与某个特定的似然函数结合后,后验分布仍然属于与先验分布相同的分布族。这种性质使得贝叶斯推断中的计算变得更加简便,因为后验分布的形式与先验分布一致。常见的共轭先验分布对:

  • 二项分布和Beta分布
  • 泊松分布和Gamma分布
  • 正态分布和正态分布
  • 多项分布和狄利克雷分布

考虑完整数据的似然函数:
在这里插入图片描述

设定参数的先验分布:
在这里插入图片描述

推出联合后验分布:
在这里插入图片描述

推导各参数的全条件后验分布:
在这里插入图片描述

实际迭代:
在这里插入图片描述

五、实施

(一)混合和收敛

在马尔可夫链蒙特卡罗(Markov Chain Monte Carlo,MCMC)算法中,混合(mixing)和收敛(convergence)是相关但不同的概念。

  • 混合指的是马尔可夫链在状态空间中有效地探索和转移的能力。
    • 混合良好的链可以自由、快速地在状态空间中移动,访问不同的区域并从目标分布的不同部分进行采样。这表明马尔可夫链能够高效地探索分布并生成代表性的样本。
    • 混合不佳意味着链在某些区域陷入困境,无法充分探索分布的整个范围,并可能产生偏误或不准确的样本。(样本自相关性太强)
  • 收敛是指马尔可夫链随着迭代次数增加,逐渐接近并稳定在目标分布周围的特性。
    • 收敛意味着马尔可夫链生成的样本随着算法的进行越来越能够代表目标分布。
    • 它表明算法已经达到一种状态,在进一步迭代中估计到的分布不会显著改变。
  • 混合是收敛的前提条件。如果马尔可夫链混合不好,它将无法收敛到目标分布。然而,即使链混合良好,也不能保证收敛。收敛需要良好的混合以及足够的迭代次数,以确保链充分探索状态空间并稳定在目标分布周围。

(二)相关术语

  • 预烧 (burn-in): 我们通常假定 MCMC 要经过一段时间的迭代才能收敛到平稳分布,这段过程我们称为 burn-in.
  • 轨迹图 (trace plot): 画出每次迭代时参数的值.
  • 对数似然函数或后验分布函数图:随着迭代,对数似然函数的变化.
  • 多链: 使用不同初值的多条短链画出变量的轨迹图,观测 f 的主要特征 (比如多峰,高度集中的支撑域). 之后选取一个好的初始值,运作一个相当长的单链计算并公布结果. (一般由最终 likelihood 大小筛选)
  • 自相关性图: 描述样本序列在不同迭代延迟下的相关性.

在这里插入图片描述
为确定预烧期和运行长度,Gelman 和 Rubin 提出一种统计量判断MCMC 是否已经收敛到平稳分布:

假设感兴趣的变量是 X X X, 其中 x 1 ( j ) , x 2 ( j ) , . . . x^{(j)}_1, x^{(j)}_2, . . . x1(j),x2(j),...是第 j j j 个马尔可夫链的样本,并假设有 J J J 个链并行运行:

  • 对于每个链,首先丢弃 D D D 值作为“预烧”并保留剩余的 L L L 值, x D ( j ) x^{(j)}_D xD(j) , x D + 1 ( j ) , . . . , x D + L − 1 ( j ) x^{(j)}_{D+1},...,x^{(j)}_{D+L−1} xD+1(j),...,xD+L1(j)
  • 计算:
    在这里插入图片描述
  • Gelman-Rubin 统计量是:
    R = L − 1 L W + 1 L B W R=\frac{\frac{L-1}{L}W+\frac{1}{L}B}{W} R=WLL1W+L1B

L → ∞ L → ∞ L 并且 B → 0 B → 0 B0 时, R R R 是趋近于 1 的。 实际应用中,某些学者建议可以接受 R < 1.2 \sqrt R < 1.2 R <1.2。但使用这种方法有一些潜在的困难:当 f 是多峰分布的情况下, 如何选择合适的初始值也许较为困难, 如果选择不恰当, 则会导致大部分的链都长期停留在同样的子域或者峰的附近。

稳妥方法:结合轨迹图和对数似然函数图,多条链进行肉眼观测分析。

.

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/pingmian/19002.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

echarts性能优化

echarts数据量多的时候优化方案&#xff1a; 渲染的数据太多时&#xff0c;渲染的速度会变慢。 let data [];for (let i 0; i < 100000; i) {let style {};if (i % 2 0) {style.color "red";}data.push({value: i,itemStyle: style,}); } myEcharts init(c…

STM32-13-MPU

STM32-01-认识单片机 STM32-02-基础知识 STM32-03-HAL库 STM32-04-时钟树 STM32-05-SYSTEM文件夹 STM32-06-GPIO STM32-07-外部中断 STM32-08-串口 STM32-09-IWDG和WWDG STM32-10-定时器 STM32-11-电容触摸按键 STM32-12-OLED模块 文章目录 STM32-12-MPU1. 内存保护单元MPU1. M…

交换机的三层交换技术

现有pc1与pc2不在同一个网段之下&#xff0c;通过交换机相连接。 进人交换机1&#xff0c;创建两个vlan 10和vlan 20 &#xff0c;进入串口2设置串口模式为access&#xff0c;并且设置默认vlan为10.进入串口3设置串口模式为access&#xff0c;并且设置默认vlan为20. 进入串口1…

深度解析搜索引擎广告(SEM)与社交媒体广告(SMM):NetFarmer助力企业数字化出海

在当今数字化时代&#xff0c;企业出海已经成为了一个必然趋势。然而&#xff0c;如何有效地在海外市场中推广品牌、吸引潜在客户&#xff0c;成为了众多企业面临的重要挑战。搜索引擎广告&#xff08;SEM&#xff09;和社交媒体广告&#xff08;SMM&#xff09;作为两种主要的…

如何下载b站(哔哩哔哩bilibili)的学习视频教程

方法1&#xff1a; 打开粘贴视频链接下载即可哔哩哔哩(bilibili)视频解析下载 - 保存B站视频到手机、电脑哔哩哔哩高清视频解析下载工具是一个免费的B站视频在线解析提取工具,支持提取B站APP和bilibili网站上的任何视频,提取出来的视频无水印.我们可以借助此下载器方便地将视频…

最大公约数和最小公倍数(函数)(C语言)

一、运行结果&#xff1b; 二、源代码&#xff1b; # define _CRT_SECURE_NO_WARNINGS # include <stdio.h>//声明函数&#xff1b; //最大公约数&#xff1b; int greatdivisor(int x, int y);//最小公倍数&#xff1b; int leastmultiple(int x, int y);int main() {/…

618精选编程书单:学好代码是用好大模型的基础

大家好&#xff0c;我是爱编程的喵喵。双985硕士毕业&#xff0c;现担任全栈工程师一职&#xff0c;热衷于将数据思维应用到工作与生活中。从事机器学习以及相关的前后端开发工作。曾在阿里云、科大讯飞、CCF等比赛获得多次Top名次。现为CSDN博客专家、人工智能领域优质创作者。…

如何选择适合自己需求的云服务器

最近明月接了一个跨境电商的代维业务&#xff0c;发现他们的云服务器很有代表性&#xff0c;今天就以此为例给大家分享一下应该如何选择适合自己需求的云服务器。像明月这样专做代维业务的可以说什么云服务器都体验过了&#xff0c;也发现大家在选择自己的云服务器的时候有很大…

加密资产私钥安全完整手册(一) ,bitget钱包为例

比特币和以太坊等加密货币的兴起开创了数字金融的新时代&#xff0c;但也带来了独特的安全挑战。这些代表现实世界价值的数字资产已成为黑客和窃贼的主要目标。为了安全地应对这种情况&#xff0c;了解私钥的基本概念至关重要。 私钥是加密货币所有权和安全性的基石。它们相当于…

VSCode小技巧,忽略不想格式化的代码行

零&#xff0e;格式化工具文档 1 . Black Ignoring sections功能 2 . autopep8 disabling-line-by-line功能&#xff1b;&#xff1b;–line-range选项 3 . Prettier prettier-ignore功能(例&#xff1a;适用于JS的// prettier-ignore&#xff0c;适用于CSS的/* prettier-igno…

HTML新春烟花盛宴

目录 写在前面 烟花盛宴 完整代码 修改文字

三步问题 | 动态规划

1.三步问题 题目连接&#xff1a;面试题 08.01. 三步问题 三步问题。有个小孩正在上楼梯&#xff0c;楼梯有n阶台阶&#xff0c;小孩一次可以上1阶、2阶或3阶。实现一种方法&#xff0c;计算小孩有多少种上楼梯的方式。结果可能很大&#xff0c;你需要对结果模1000000007。 2…

C语言 指针——指针变量的定义、初始化及解引用

目录 指针 内存如何编址&#xff1f; 如何对变量进行寻址&#xff1f; 用什么类型的变量来存放变量的地址? 如何显示变量的地址?​编辑 使用未初始化的指针会怎样&#xff1f; NULL是什么&#xff1f; 如何访问指针变量指向的存储单元中的数据&#xff1f; 指针变量的…

APP原生开发与框架开发的优劣势

电话管家APP商用也有几年时间了&#xff0c;但是客户一直都有遇到一些问题。 为什么我们的APP老是要升级&#xff1f; 为什么有些手机使用体验不好&#xff1f; 为什么有些公司的APP几天就开发出来上线了&#xff1f; 为什么有些公司的APP那么便宜&#xff1f; 今天就来从…

家政预约小程序08服务详情

目录 1 创建页面2 创建URL参数3 配置数据详情组件4 从分类页跳转到详情页5 搭建详情页总结 现在我们的小程序已经在首页和分类页展示了服务的列表信息&#xff0c;当用户点击具体的内容的时候需要打开详情页&#xff0c;本篇介绍一下详情页的开发。 1 创建页面 打开应用编辑器…

中学生学人工智能系列:如何用AI学英语

经常有读者朋友给公众号《人工智能怎么学》留言咨询如何使用人工智能学习语文、数学、英语等科目。这些都是中学教师、中学生朋友及其家长们普遍关注的问题。仅仅使用留言回复的方式&#xff0c;不可能对这些问题做出具体和透彻的解答&#xff0c;因此本公众号近期将推出中学生…

如何在phpMy管理对Joomla后台的登录密码进行重置

本周有一个客户&#xff0c;购买Hostease的虚拟主机&#xff0c;询问我们的在线客服&#xff0c;如何在phpMy管理对Joomla后台的登录密码进行重置&#xff1f;我们为用户提供相关教程&#xff0c;用户很快解决了遇到的问题。在此&#xff0c;我们分享这个操作教程&#xff0c;希…

LeetCode题练习与总结:平衡二叉树--110

一、题目描述 给定一个二叉树&#xff0c;判断它是否是平衡二叉树。 示例 1&#xff1a; 输入&#xff1a;root [3,9,20,null,null,15,7] 输出&#xff1a;true示例 2&#xff1a; 输入&#xff1a;root [1,2,2,3,3,null,null,4,4] 输出&#xff1a;false示例 3&#xff1a…

【Java用法】java中计算两个时间差

java中计算两个时间差 不多说&#xff0c;直接上代码&#xff0c;可自行查看示例 package org.example.calc;import java.time.LocalDateTime; import java.time.format.DateTimeFormatter; import java.time.temporal.ChronoUnit;public class MinusTest {public static void…

97.网络游戏逆向分析与漏洞攻防-ui界面的设计-通过逆向分析确认角色信息

免责声明&#xff1a;内容仅供学习参考&#xff0c;请合法利用知识&#xff0c;禁止进行违法犯罪活动&#xff01; 如果看不懂、不知道现在做的什么&#xff0c;那就跟着做完看效果&#xff0c;代码看不懂是正常的&#xff0c;只要会抄就行&#xff0c;抄着抄着就能懂了 内容…