已知样本服从多项分布运用贝叶斯循环迭代更新后验分布的代码

假设我们有 $n$ 个样本，每个样本有 $k$ 种可能的离散取值。我们用多项分布来描述这些样本，其中 $\boldsymbol{\theta}=(\theta_1,\theta_2,\dots,\theta_k)$ 表示每个类别的概率，也就是多项分布的参数。假设我们有先验知识 $p(\boldsymbol{\theta})$，那么我们可以通过贝叶斯定理得到后验分布：

$$
p(\boldsymbol{\theta}|\mathbf{x})=\frac{p(\mathbf{x}|\boldsymbol{\theta})p(\boldsymbol{\theta})}{p(\mathbf{x})}
$$

其中，$\mathbf{x}=(\mathbf{x}_1,\mathbf{x}_2,\dots,\mathbf{x}_n)$ 表示样本集合，$\mathbf{x}_i$ 是一个 $k$ 维向量，表示第 $i$ 个样本的类别。

假设我们采用狄利克雷分布作为先验分布，即

$$
p(\boldsymbol{\theta})=\frac{1}{\mathrm{B}(\alpha)}\prod_{i=1}^k\theta_i^{\alpha_i-1}
$$

其中，$\mathrm{B}$ 表示贝塔函数，$\alpha=(\alpha_1,\alpha_2,\dots,\alpha_k)$ 是一个 $k$ 维向量，表示先验分布的超参数。可以证明，狄利克雷分布是多项分布的共轭先验分布，因此，后验分布也是狄利克雷分布，即

$$
p(\boldsymbol{\theta}|\mathbf{x})=\frac{1}{\mathrm{B}(\boldsymbol{\alpha}+\mathbf{n})}\prod_{i=1}^k\theta_i^{\alpha_i+n_i-1}
$$

其中，$\boldsymbol{\alpha}$ 是先验分布的超参数向量，$\mathbf{n}=(n_1,n_2,\dots,n_k)$ 是样本中每个类别出现的次数，即

$$
n_i=\sum_{j=1}^n\mathbf{1}_{\{\mathbf{x}_j=i\}}
$$

这个后验分布的参数可以通过贝叶斯循环迭代更新来求解。具体来说，假设我们已经有一个先验分布 $p(\boldsymbol{\theta}^{(t-1)})$，其中 $\boldsymbol{\theta}^{(t-1)}$ 表示第 $t-1$ 次迭代后的参数值。我们从样本 $\mathbf{x}$ 中抽取一个样本 $\mathbf{x}_i$，然后计算其对应的似然函数 $p(\mathbf{x}_i|\boldsymbol{\theta}^{(t-1)})$ 和先验分布 $p(\boldsymbol{\theta}^{(t-1)})$，然后通过贝叶斯定理计算出后验分布 $p(\boldsymbol{\theta}^{(t)}|\mathbf{x}_i)$。由于 $\mathbf{x}_i$ 是从样本 $\mathbf{x}$ 中随机抽取的，因此我们可以多次重复这个过程，直到后验分布收敛。具体的迭代公式如下：

$$
\begin{aligned}
\alpha_i^{(t)}&=\alpha_i^{(t-1)}+n_i\\
\boldsymbol{\theta}^{(t)}&\sim\mathrm{Dirichlet}(\boldsymbol{\alpha}^{(t)})
\end{aligned}
$$

其中，$\mathrm{Dirichlet}$ 表示狄利克雷分布。这个迭代过程可以用代码实现，如下所示：

import numpy as np
from scipy.stats import dirichlet

def posterior_update(alpha, x, num_iters=1000, burn_in=100):
    k = alpha.shape[0]
    n = x.shape[0]
    n_k = np.zeros(k)
    for i in range(n):
        n_k[x[i]] += 1
    alpha_t = alpha + n_k
    theta_t = dirichlet.rvs(alpha_t)
    for t in range(num_iters):
        i = np.random.randint(n)
        x_i = x[i]
        p = theta_t[x_i] * (alpha_t[x_i] + n_k[x_i]) / (alpha_t.sum() + n)
        if np.random.rand() < p:
            n_k[x_i] += 1
        alpha_t = alpha + n_k
        theta_t = dirichlet.rvs(alpha_t)
        if t > burn_in and t % 10 == 0:
            print('iter={}, theta={}, alpha={}'.format(t, theta_t, alpha_t))
    return alpha_t, theta_t

其中，`alpha` 是先验分布的超参数向量，`x` 是样本集合，`num_iters` 表示总共迭代的次数，`burn_in` 表示前面几次迭代不考虑。我们首先计算 $\mathbf{n}$，然后通过循环迭代更新 $\alpha$ 和 $\boldsymbol{\theta}$，最后返回后验分布的参数。在每次迭代时，我们随机抽取一个样本 $\mathbf{x}_i$，然后计算其对应的似然函数和先验分布，并根据贝叶斯定理计算出后验分布。如果后验概率比当前参数更高，则接受这个样本，否则保持原来的参数不变。最后，我们输出一些中间结果，以便观察迭代过程。