gibbs samplingのコード供養

ポアソン混合モデルの gibbs sampling を書いたので供養のため以下に載せておきます。

    for i in range(MAXITER):
        for n in range(N):
            # (4.37), (4.38)
            for k in range(K):
                etas[n][k] = exp(Xs[n] * ln(lambdas[k]) - lambdas[k] + ln(pis[k]))
            # normalization
            subtotal = sum(etas[n])
            for k in range(K):
                etas[n][k] /= subtotal
            #
            Ss[n] = np.random.multinomial(n=1, pvals=etas[n])
        for k in range(K):
            # (4.42)
            hat_a[k] = sum(Ss[n][k] * Xs[n] for n in range(N)) + a
            hat_b[k] = sum(Ss[n][k] for n in range(N)) + b
            # (4.41)
            lambdas[k] = np.random.gamma(shape=hat_a[k], scale=1 / hat_b[k])
        # (4.44)
        for k in range(K):
            hat_alphas[k] = sum(Ss[n][k] for n in range(N)) + alphas[k]
        pis = np.random.dirichlet(alpha=hat_alphas, size=1).reshape(K)