最近，我使用numpy在Python上实现了LDA主题模型的Gibbs采样，参考了一些网站上的代码。在Gibbs采样的每次迭代中，我们移除一个（当前的）词，根据从LDA模型中推断出的后验条件概率分布为该词采样一个新主题，并更新词-主题计数，如下所示：

for m, doc in enumerate(docs): #m: 文档ID
    for n, t in enumerate(doc): #n: 文档内词的ID, t: 全局词的ID
        # 减少与主题z关联的词t的计数
        z = z_m_n[m][n]
        n_m_z[m][z] -= 1
        n_z_t[z, t] -= 1 
        n_z[z] -= 1
        n_m[m] -= 1
        # 为多项式采样新主题
        p_z_left = (n_z_t[:, t] + beta) / (n_z + V * beta)
        p_z_right = (n_m_z[m] + alpha) / ( n_m[m] + alpha * K)
        p_z = p_z_left * p_z_right
        p_z /= numpy.sum(p_z)
        new_z = numpy.random.multinomial(1, p_z).argmax() 
        # 将z设置为新主题并增加计数
        z_m_n[m][n] = new_z
        n_m_z[m][new_z] += 1
        n_z_t[new_z, t] += 1
        n_z[new_z] += 1
        n_m[m] += 1

在上面的代码中，我们使用多项式scipy函数采样一个新的（单一的）z。

现在，我想实现这篇论文中的联合情感主题模型。现在，我需要以下结构来跟踪所需的计数：

包含每个主题、每个情感的词出现次数的3D矩阵
包含每个情感、每个文档的每个主题出现次数的3D矩阵
包含每个情感的每个主题出现次数的2D矩阵
包含每个文档的每个情感出现次数的2D矩阵

现在问题来了：在这个Gibbs采样器中，对于在文档中看到的每个词，现在都需要从条件后验分布中同时采样一个新主题和一个情感标签（论文第4页的方程5）。我现在如何在Python中“采样这两个值”？

提前感谢…

回答：

试试这个。从主题和情感标签的联合分布中采样只意味着整个T x S矩阵的总和应为1。

docs=[[0,1],[0,0],[1,0,1]]
D=len(docs)
z_d_n=[[0 for _ in xrange(len(d))] for d in docs]
l_d_n=[[0 for _ in xrange(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=.1
alpha=.1
gamma=.1
for d, doc in enumerate(docs): #d: 文档ID
    for n, m in enumerate(doc): #i: 文档内词的索引, m: 词汇表中词的ID
        # j是主题
        j = z_d_n[d][n]
        # k是情感
        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): #d: 文档ID
    for n, m in enumerate(doc): #i: 文档内词的索引, m: 词汇表中词的ID
        # j是主题
        j = z_d_n[d][n]
        # k是情感
        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 
        
        # 联合采样一个新主题和情感标签
        # T是主题数量
        # S是情感数量
        p_left = (n_m_j_k[m] + beta) / (n_j_k + V * beta) # T x S数组
        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