docs = [[0, 1], [0, 0], [1, 0, 1]]
D = len(docs)
z_d_n = [[(0) for _ in range(len(d))] for d in docs]
l_d_n = [[(0) for _ in range(len(d))] for d in docs]
V = 2
T = 2
S = 2
n_m_j_k = numpy.zeros((V, T, S))
n_j_k_d = numpy.zeros((T, S, D))
n_j_k = numpy.zeros((T, S))
n_k_d = numpy.zeros((S, D))
n_d = numpy.zeros(D)
beta = 0.1
alpha = 0.1
gamma = 0.1
for d, doc in enumerate(docs):
    for n, m in enumerate(doc):
        j = z_d_n[d][n]
k = l_d_n[d][n]
n_m_j_k[m][j][k] += 1
n_j_k_d[j][k][d] += 1
n_j_k[j][k] += 1
n_k_d[k][d] += 1
n_d[d] += 1
for d, doc in enumerate(docs):
    for n, m in enumerate(doc):
        j = z_d_n[d][n]
k = l_d_n[d][n]
n_m_j_k[m][j][k] -= 1
n_j_k_d[j][k][d] -= 1
n_j_k[j][k] -= 1
n_k_d[k][d] -= 1
n_d[d] -= 1
p_left = (n_m_j_k[m] + beta) / (n_j_k + V * beta)
p_mid = (n_j_k_d[:, :, (d)] + alpha) / numpy.tile(n_k_d[:, (d)] + T * alpha,
    (T, 1))
p_right = numpy.tile(n_k_d[:, (d)] + gamma, (T, 1)) / numpy.tile(n_d[d] + S *
    gamma, (T, S))
p = p_left * p_mid * p_right
p /= numpy.sum(p)
new_jk = numpy.random.multinomial(1, numpy.reshape(p, T * S)).argmax()
j = new_jk / T
k = new_jk % T
z_d_n[d][n] = j
l_d_n[d][n] = k
n_m_j_k[m][j][k] += 1
n_j_k[j][k] += 1
n_k_d[k][d] += 1
n_d[d] += 1
Implementing Topic Model with Python (numpy)