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)