系列文章目录

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

---

一、基本概念

（一）马尔科夫链

1、定义

考虑在每个时间段有一个值的随机过程，令$X_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​=j∣Xn​=i,Xn−1​=in−1​,...,X1​=i1​}P{Xn+1​=j∣Xn​=i}Pij​​
这样的随机过程称为马尔科夫链。$P$为一步转移概率$P_{ij}$的矩阵，n次转移后的矩阵为$P^n$。

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

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

2、性质

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

• 常返：称能以概率1回来的状态是常返的：
  $$f_{ii}=\sum _{i=1}^{\infty }{f}_{ii}^{(n)}=1$$
  其中$f_{ij}^{(n)}=p(X^{(n)}=j,X_k\ne j,k=1,2,...,n-1|X_0=i)$为马氏链在0时从状态 $i$ 出发，经过n步转移后，首次到达状态 $j$ 的概率。若返回概率小于1，即$f_{ii}<1$，则为非常返态。

• 正常返和零常返：一个状态的平均返回时间 ${\mu }_i$ 为
  $${\mu }_i=\sum _{n=1}^{\infty }n{f}_{ii}^{(n)}$$
  若${\mu }_i<\infty$，则为正常返，反之${\mu }_i=\infty$为零常返。如果状态空间有限，其常返状态都是非零常返

• 周期性：

  • 如果马尔可夫链只能以一定的规则间隔访问状态空间的某些部分，则它是周期性的。如果由状态 $j$ 经非 $d$ 整数倍步到达$j$ 的概率为 0，则称状态 $j$ 具有周期 $d$，周期可定义为集合 { n : n ≥ 0 , p i i ( n ) > 0 } \{n : n ≥ 0, p^{(n)}_{ii} > 0\} {n:n≥0,pii(n)​>0} 的最大公约数，即所有返回可能的次数的最大公约数。
  • 如果一条马氏链的每一个状态的周期都为 1, 则称此链为非周期的。

• 遍历的：如果一条马氏链是不可约、非周期，且其所有状态都是非零常返的，则称之为遍历的
  

3、常返

设 $i$ 是常返态，则：

• $i$ 是零常返态的充要条件为：$li{m}_{n\to \infty }{p}_{ii}^{(n)}=0$
• $i$ 是正常返态的充要条件为：$li{m}_{n\to \infty }{p}_{ii}^{(n)}>0$
  

4、平稳分布

平稳分布定义：若一离散分布 $\pi$ 满足 ${\pi }^{T}P={\pi }^T$，则称之为转移概率矩阵为 $P$ 的马氏链的平稳分布。

• （充分条件）如果一条时间齐性的马氏链满足：
  $${\pi }_i{p}_{ij}={\pi }_j{p}_{ji},\forall i,j\in S,$$
  则$\pi$是此链的平稳分布，且称此链为可逆的，上述方程也称为细致平衡。

如果一个转移概率阵为$P$，平稳分布为$\pi$的马氏链是不可约的，正常返的，且非周期的（遍历），则$\pi$唯一，且满足：
$$\underset{n\to \infty }{\lim }P[X_n=j]={\pi }_j$$
也就是极限分布即为平稳分布。其中${\pi }_j$是$\pi$的第$j$个元素，且满足如下方程组：
$${\pi }_j\ge 0,\sum _{i\in S}{\pi }_i=1,且{\pi }_j=\sum _{i\in S}{\pi }_i{p}_{ij},\forall j\in S$$

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

（二）MCMC原理

1、核心思想

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

2、连续状态

在连续分布情况下，对任一可测集$B$，一步转移概率定义为：
$$P(x\to B)={\int }_{B}p(x,x^{\prime })d{x}^{\prime }$$
转移概率$p(\cdot ,\cdot )$称为Markov链转移核。通常假定$p(\cdot ,\cdot )$与$t$无关，即基于该转移核的Markov链是时间齐次的。

例：根据转移核$P(X^{(t+1)}|X^{(t)})\sim N(0.5X^{(t)},1)$，产生平稳分布是$N(0,4/3)的Markov链。

3、MCMC估计期望的步骤

MCMC方法估计$E_{\pi }f(X)=\int f(x)\pi (x)dx$的基本流程：

• 选择转移核$p(\cdot ,\cdot )$（参数更新公式），使其平稳分布是$\pi (x)$
• 从某一点$X^{(0)}$出发，用上述转移核$p(\cdot ,\cdot )$产生Markov链序列$X^{(0)},X^{(1)},...X^{(n)}$
• 对较大的$n$，选择合适的$m$，$E_{\pi }f(X)=\int f(x)\pi (x)dx$的估计为：
  $${\hat{E}}_{\pi }f=\frac{1}{n-m}\sum _{t=m+1}^{n}f(X^{(t)})$$

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

二、满条件分布

1、定义

MCMC方法中转移核$p(x,x^{\prime })$的构造大多建立在形如$\pi (x_T|x_{-T})$的条件分布上：

• x T = { x i , i ∈ T } x_T=\{x_i,i\in T\} xT​={xi​,i∈T}， x − T = { x i , i ∉ T } x_{-T}=\{x_i,i \notin T\} x−T​={xi​,i∈/T}， T ⊂ I = { 1 , ⋅ ⋅ ⋅ , p } T⊂ I = \{1, · · · , p\} T⊂I={1,⋅⋅⋅,p}
• 在条件分布$\pi (x_T|x_{-T})$中，p个所有的变量或者出现在条件中，或者出现在变元中，这种条件分布称为满条件分布
  $$\pi (x_T|x_{-T})=\frac{\pi (x)}{\int \pi (x)d{x}_T}\propto \pi (x)$$
  进而对任意$x,x^{\prime }\in X$且$x_{-T}=x_{-T}^{\prime }$，有：
  $$\frac{\pi (x_T^{\prime }|x_{-T}^{\prime })}{\pi (x_T|x_{-T})}=\frac{\pi (x^{\prime })}{\pi (x)}$$

2、考虑MCMC中的应用

一般情况下$y=(x,z,\theta )$，其中$x$表示观测数据，$z$表示缺失数据，$\theta$表示参数。令$p(x,z|\theta )$表示完全数据的密度函数，$\pi (\theta )$表示$\theta$的先验分布，$f(y)=p(x,z|\theta )\pi (\theta )$，则$y$的满条件分布为：
$$\begin{gathered} f(z_i|z_{-i},x,\theta )\propto p(x,z|\theta ) \\ \pi ({\theta }_i|{\theta }_{-i},x,z)\propto f(y)\propto p(x,z|\theta )\pi ({\theta }_i|{\theta }_{-i}) \\ f(x_i|x_{-i},z,\theta )\propto p(x,z|\theta ) \end{gathered}$$
其中 θ − 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},z−i​={zj​:j=i},x−i​={xj​:j=i}。

例：设$x=(x_1,x_2)$的密度函数为：
π ( 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​(x1​−1)2(x2​−1)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}) π(x1​∣x2​)∝π(x1​,x2​)∝exp{−21​(x1​−1)2(x2​−1)2}=N(1,(x2​−1)−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}) π(x2​∣x1​)∝π(x1​,x2​)∝exp{−21​(x2​−1)2(x1​−1)2}=N(1,(x1​−1)−2)

3、伽玛分布

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

例：设 $y_1,...,y_n$ 独立同分布，来自正态分布 $N(\mu ,{\tau }^{-1})$。 参数$\mu ,{\tau }^{-1}$的先验分布分别为正态分布 $\mu \sim N(0,1)$, 伽马分布 $\tau \sim \Gamma (2,1)$, 且 $\mu$ 与 $\tau$ 独立，计算满条件分布。

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

例：样本$y_i,i=1,...,n$独立且$y_i\sim B(1,p_i)$（伯努利分布），其中$p_i=(1+exp(-(\alpha +\beta x_i)){)}^{-1}$，$x_i$假设是固定的，已知的。参数$\alpha ,\beta$的先验分布分别为$\alpha \sim N(0,1)$，且$\alpha ,\beta$独立。

三、Metropolis–Hastings 算法

（一）基本概念

1、思想

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

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

2、具体实施

• 如果链在时刻$t$处于状态$x$，即$X_t=x$
• 由$q(\cdot |x)$产生一个潜在的转移$x\to x^{\prime }$
• 根据概率$\alpha (x,x^{\prime })$决定是否转移

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

3、算法过程

假设$\pi (x)$为目标概率分布，MH算法的过程为：

• 初始化：选定初始状态$x_0$，令$t=0$
• 迭代过程：
  • 生成：从某一容易抽样的分布$q(x^{\prime }|x_t)$中随机生成候选状态$x^{\prime }$
  • 计算：计算是否采纳候选状态的概率
    $$\alpha (x_t,x^{\prime })=min(1,\frac{\pi (x^{\prime })}{\pi (x_t)}\frac{q(x_t|x^{\prime })}{q(x^{\prime }|x_t)})$$
  • 接受或拒绝
    • 从$[0,1]$的均匀分布中生成随机数$u$
    • 如果$u\le \alpha (x_t,x^{\prime })$，则接受该状态，并令$x_{t+1}=x^{\prime }$
    • 如果$u>\alpha (x_t,x^{\prime })$，则拒绝该状态，并令$x_{t+1}=x_t$
  • 增量：令$t=t+1$

MH算法的推导：

（二）Metropolis 选择（提案分布 $q(x,x\prime )$ 的选取）

• Metropolis 建议 $q(x,x\prime )$ 为对称分布，即：
  $$q(x,x^{\prime })=q(x^{\prime },x),\forall x,x\prime$$
  • $\alpha (x,x^{\prime })$简化为：
    α ( 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^{\prime })=f(|x-x^{\prime }|)$，被称为随机移动Metropolis算法，具体例子为 q ( x , x ′ ) ∝ e x p { − ( x ′ − x ) 2 / e } q(x,x')∝exp\{-(x'-x)^2/e\} q(x,x′)∝exp{−(x′−x)2/e}

例：生成一个Markov链，使得其平稳分布为柯西分布，
$$f(x)=\frac{1}{\pi }\frac{1}{1+x^2}$$
选定的提案分布$q(x,x^{\prime })$是$N(x,b^2)$，其中$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^{\prime })$与当前状态$x$无关，即$q(x,x^{\prime })=q(x^{\prime })$，则由此提案分布导出的MH算法称为独立抽样。此时的$\alpha (x,x^{\prime })$为：
  α ( 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)$接近$\pi (x)$，基于独立抽样获取的Markov链的收敛效果更好。

给定数据 $Y_1,...,Y_n\stackrel{i.i.d}{\sim }N(\theta ,1)$， 先验分布 $\pi (\theta )=1/\pi (1+{\theta }^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,\bar{y}=0.14$，此时使用$x$和$x^{\prime }$的记号，$\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(x−0.14)2/2}×1+x21​
求后验期望$E(x)$

（三）单元素 Metropolis-Hastings 算法

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

单元素 Metropolis-Hastings 算法的核心思想：对于 $p$ 维变量 $X$, 基于 $p-1$维变量 $X_{-i}$的条件分布 $X_i|X_{-i},i=1,...,p$，选择转移核$q_i(x_i\to x_i^{\prime }|X_{-i}=x_{-i})$。由转移核$q_i(x_i\to x_i^{\prime }|X_{-i}=x_{-i})$ 产生可能的 $x_i^{\prime }$, 以概率
$${\alpha }_i(x_i\to x_i^{\prime }|x_{-i})=min(1,\frac{\pi (x^{\prime })q_i(x_i^{\prime }\to x_i|x_{-i})}{\pi (x)1_i(x_i\to x_i^{\prime }|x_{-i})})$$
决定是否接受$x^{\prime }$为链的下一状态。

即每次只更新一个元素，其他保持不变：

四、Gibbs 算法

（一）基本概念

Gibbs 抽样是一种单元素 Metropolis-Hastings 算法的特殊情况，单元素Metropolis-Hastings 算法中取$q_i(x_i\to x_i^{\prime }|X_{-i}=x_{-i})$为$\pi (x_i|x_{-i})$。（也就是说直接定义为其满条件分布）
此时的接受率为1，这意味着，我们不需要舍弃样本，每个更新后的值都为样本。

（二）采样过程

Gibbs采样的过程：

• 确定初始值$X^{(1)}$
• 假设已得到样本$X^{(i)}$，记下一个样本为$X^{(i+1)}=(x_1^{(i+1)},x_2^{(i+1)},...,x_n^{(i+1)})$，对其中某一分量$x_j^{(i+1)}$可通过在其他分量已知的条件下该分量的概率分布来抽取该分量。对于此条件概率，我们使用样本$X^{(i+1)}$中已得到的分量$x_1^{(i+1)}$到$x_{j-1}^{(i+1)}$以及上一样本$X^{(i)}$中的分量$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)})$$
• 重复上述过程
  

（三）延展

• 更新顺序：
  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)})$
  • 用某Metropolis迭代更新 $(X_2^{(t+1)},X_3^{(t+1)})|(x_1^{(t+1)},x_4^{(t)},x_5^{(t)},x_6^{(t)})$
  • 用某Metropolis迭代更新 $X_4^{(t+1)}|x_1^{(t+1)},x_2^{(t+1)}|(x_3^{(t+1)},x_5^{(t)},x_6^{(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 的一个或者多个元素的一元边际密度没有显示表达的时候，Gibbs算法中的 Metropolis-Hastings 迭代特别有用. 有时也是 Gibbs 跳出局部最优的好方法

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

（四）示例

1、简单正态例子：

（已知参数，但假设我们只能单独产生单变量正态分布的随机数，如何从二维正态采样）给定一个目标分布：
$$X=\left[\begin{array}{c}{X}_1\\ X_2\end{array}\right]\sim N(\left[\begin{array}{c}{u}_1\\ {\mu }_2\end{array}\right],\left[\begin{array}{cc}{\sigma }_1^2& {\sigma }_{12}\\ {\sigma }_{12}& {\sigma }_2^2\end{array}\right])$$

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!​p1m1​pxm2​...pnmn​​
其中$N={\sum }_{i=1}^n{m}_i$，某实验服从上述多项分布，$N=22,n=7$，7个结果出现的概率分别为：
$$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 ))$$
现有观测数据为$y=(y_1,y_2,y_3,y_4,y_5)=(14,1,1,1,5)$，且
$$(z_1,y_1-z_1,y_2,y_3,z_2,y_4-z_2,y_5)\sim M(22;p)$$
其中$M$表示多项分布

3、分类消费者例子

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

• 缺失数据$c_j$的边际分布：$p(c_j=k)={\pi }_k,{\sum }_{k=1}^K{\pi }_k=1$
• ${\lambda }_{gk}$表示第k类消费者买第g类物品的均值：$X_{gj}|c_j=k\sim Pois({\lambda }_{gk})$

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

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

考虑完整数据的似然函数：

设定参数的先验分布：

推出联合后验分布：

推导各参数的全条件后验分布：

实际迭代：

五、实施

（一）混合和收敛

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

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

（二）相关术语

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

为确定预烧期和运行长度，Gelman 和 Rubin 提出一种统计量判断MCMC 是否已经收敛到平稳分布：

假设感兴趣的变量是 $X$， 其中 $x_1^{(j)},x_2^{(j)},...$是第 $j$ 个马尔可夫链的样本，并假设有 $J$ 个链并行运行：

• 对于每个链，首先丢弃 $D$ 值作为“预烧”并保留剩余的 $L$ 值，$x_D^{(j)}$ , $x_{D+1}^{(j)},...,x_{D+L-1}^{(j)}$
• 计算：
  
• Gelman-Rubin 统计量是：
  $$R=\frac{\frac{L-1}{L}W+\frac{1}{L}B}{W}$$

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

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

.