3月 05

Reading Note : Parameter estimation for text analysis 暨LDA学习小结

伟大的Parameter estimation for text analysis!当把这篇看的差不多的时候,也就到了LDA基础知识终结的时刻了,意味着LDA基础模型的基本了解完成了。所以对该模型的学习告一段落,下一阶段就是了解LDA无穷无尽的变种,不过那些不是很有用了,因为LDA已经被人水遍了各大“论坛”……

抛开LDA背后复杂深入的数学背景不说,光就LDA的内容,确实不多,虽然变分法还是不懂,不过现在终于还是理解了“LDA is just a simple model”这句话。

总结一下学习过程:
1.概率的基本概念:CDF、PDF、Bayes’rule、各种简单的分布Bernoulli,binomial,multinomial、包括对prior、likelihood、postprior的理解(PRML1.2)
​3.概率图模型 Probabilistic Graphical Models: PRML Chapter 8 基本概念即可
4.采样算法:Basic Sampling,Sampling Methods(PRML Chapter 11),马尔科夫蒙特卡洛 MCMC,Gibbs Sampling
​6.进阶资料:《Gibbs Sampling for the Uninitiated》、本文
 
——————————————– 伟大的分割线 !PETA! ​——————————————–

一、前面无关部分

关于ML、MAP、Bayesian inference

二、模型进一步记忆

从本图来看,需要记住:

1.$\theta_m$是每一个document单独一个$\theta$,所以M个doc共有M个$\theta_m$,整个$\theta$是一个M*K的矩阵(M个doc,每个doc一个K维topic分布向量)。

2.$\varphi_k$总共只有K个,对于每一个topic,有一个$\varphi_k$,这些参数是独立于文档的,也就是对于整个corpus只sample一次。不像$\theta_m$那样每一个都对应一个文档,每个文档都不同,$\varphi_k$对于所有文档都相同,是一个K*V的矩阵(K个topic,每个topic一个V维从topic产生词的概率分布)。

就这些了。

三、推导

公式(39):$P(p|\alpha)=Dir(p|\alpha)$意思是从参数为$\alpha$的狄利克雷分布,采样一个多项分布参数p的概率是多少,概率是标准狄利克雷PDF。这里Dirichlet delta function为:

$\Delta(\vec{\alpha})=\frac{\Gamma(\alpha_1)*\Gamma(\alpha_2)*…*\Gamma(\alpha_k)}{\Gamma(\sum_1^K \alpha_k)}$

这个function要记住,下面一溜烟全是这个。

公式(43)是一元语言模型的likelihood,意思是如果提供了语料库W,知道了W里面每个词的个数,那么使用最大似然估计最大化L就可以估计出参数多项分布p。

公式(44)是考虑了先验的情形,假如已知语料库W和参数$\alpha$,那么他们产生多项分布参数p的概率是$Dir(p|\alpha+n)$,这个推导我记得在PRML2.1中有解释,抛开复杂的数学证明,只要参考标准狄利克雷分布的归一化项,很容易想出式(46)的归一化项就是$\Delta(\alpha+n)$。这时如果要通过W估计参数p,那么就要使用贝叶斯推断,用这个狄利克雷pdf输出一个p的期望即可。

最关键的推导(63)-(78):从63-73的目标是要求出整个LDA的联合概率表达式,这样(63)就可以被用在Gibbs Sampler的分子上。首先(63)把联合概率拆成相互独立的两部分$p(w|z,\beta)$和$p(z|\alpha)$,然后分别对这两部分布求表达式。式(64)、(65)首先不考虑超参数$\beta$,而是假设已知参数$\Phi$。这个$\Phi$就是那个K*V维矩阵,表示从每一个topic产生词的概率。然后(66)要把$\Phi$积分掉,这样就可以求出第一部分$p(w|z,\beta)$为表达式(68)。从66-68的积分过程一直在套用狄利克雷积分的结果,反正整篇文章套来套去始终就是这么一个狄利克雷积分。$\vec{n}_z$是一个V维的向量,对于topic z,代表每一个词在这个topic里面有几个。从69到72的道理其实和64-68一模一样了。$\vec{n}_m$是一个K维向量,对于文档m,代表每一个topic在这个文档里有几个词。

最后(78)求出了Gibbs Sampler所需要的条件概率表达式。这个表达式还是要贴出来的,为了和代码里面对应:

具体选择下一个新topic的方法是:通过计算每一个topic的新的产生概率$p(z_i=k|z_{\urcorner i},w)$也就是代码中的p[k]产生一个新topic。比如有三个topic,算出来产生新的p的概率值为{0.3,0.2,0.4},注意这个条件概率加起来并不一定是一。然后我为了按照这个概率产生一个新topic,我用random函数从uniform distribution产生一个0至0.9的随机数r。如果0<=r<0.3,则新topic赋值为1,如果0.3<=r<0.5,则新topic赋值为2,如果0.5<=r<0.9,那么新topic赋值为3。

四、代码

/*
 * (C) Copyright 2005, Gregor Heinrich (gregor :: arbylon : net) 
 * LdaGibbsSampler is free software; you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by the Free
 * Software Foundation; either version 2 of the License, or (at your option) any
 * later version.
 * LdaGibbsSampler is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
 * details.
 * You should have received a copy of the GNU General Public License along with
 * this program; if not, write to the Free Software Foundation, Inc., 59 Temple
 * Place, Suite 330, Boston, MA 02111-1307 USA
 */
import java.text.DecimalFormat;
import java.text.NumberFormat;

public class LdaGibbsSampler {
    /**
     * document data (term lists)
     */
    int[][] documents;
    /**
     * vocabulary size
     */
    int V;
    /**
     * number of topics
     */
    int K;
    /**
     * Dirichlet parameter (document--topic associations)
     */
    double alpha;
    /**
     * Dirichlet parameter (topic--term associations)
     */
    double beta;
    /**
     * topic assignments for each word.
     * N * M 维,第一维是文档,第二维是word
     */
    int z[][];
    /**
     * nw[i][j] number of instances of word i (term?) assigned to topic j.
     */
    int[][] nw;
    /**
     * nd[i][j] number of words in document i assigned to topic j.
     */
    int[][] nd;
    /**
     * nwsum[j] total number of words assigned to topic j.
     */
    int[] nwsum;
    /**
     * nasum[i] total number of words in document i.
     */
    int[] ndsum;
    /**
     * cumulative statistics of theta
     */
    double[][] thetasum;
    /**
     * cumulative statistics of phi
     */
    double[][] phisum;
    /**
     * size of statistics
     */
    int numstats;
    /**
     * sampling lag (?)
     */
    private static int THIN_INTERVAL = 20;

    /**
     * burn-in period
     */
    private static int BURN_IN = 100;

    /**
     * max iterations
     */
    private static int ITERATIONS = 1000;

    /**
     * sample lag (if -1 only one sample taken)
     */
    private static int SAMPLE_LAG;

    private static int dispcol = 0;

    /**
     * Initialise the Gibbs sampler with data.
     * 
     * @param V
     *            vocabulary size
     * @param data
     */
    public LdaGibbsSampler(int[][] documents, int V) {

        this.documents = documents;
        this.V = V;
    }

    /**
     * Initialisation: Must start with an assignment of observations to topics ?
     * Many alternatives are possible, I chose to perform random assignments
     * with equal probabilities
     * 
     * @param K
     *            number of topics
     * @return z assignment of topics to words
     */
    public void initialState(int K) {
        int i;

        int M = documents.length;

        // initialise count variables.
        nw = new int[V][K];
        nd = new int[M][K];
        nwsum = new int[K];
        ndsum = new int[M];

        // The z_i are are initialised to values in [1,K] to determine the
        // initial state of the Markov chain.
        // 为了方便,他没用从狄利克雷参数采样,而是随机初始化了!

        z = new int[M][];
        for (int m = 0; m < M; m++) {
            int N = documents[m].length;
            z[m] = new int[N];
            for (int n = 0; n < N; n++) {
            	//随机初始化!
                int topic = (int) (Math.random() * K);
                z[m][n] = topic;
                // number of instances of word i assigned to topic j
                // documents[m][n] 是第m个doc中的第n个词
                nw[documents[m][n]][topic]++;
                // number of words in document i assigned to topic j.
                nd[m][topic]++;
                // total number of words assigned to topic j.
                nwsum[topic]++;
            }
            // total number of words in document i
            ndsum[m] = N;
        }
    }

    /**
     * Main method: Select initial state ? Repeat a large number of times: 1.
     * Select an element 2. Update conditional on other elements. If
     * appropriate, output summary for each run.
     * 
     * @param K
     *            number of topics
     * @param alpha
     *            symmetric prior parameter on document--topic associations
     * @param beta
     *            symmetric prior parameter on topic--term associations
     */
    private void gibbs(int K, double alpha, double beta) {
        this.K = K;
        this.alpha = alpha;
        this.beta = beta;

        // init sampler statistics
        if (SAMPLE_LAG > 0) {
            thetasum = new double[documents.length][K];
            phisum = new double[K][V];
            numstats = 0;
        }

        // initial state of the Markov chain:
        //启动马尔科夫链需要一个起始状态
        initialState(K);

        //每一轮sample
        for (int i = 0; i < ITERATIONS; i++) {

            // for all z_i
            for (int m = 0; m < z.length; m++) {
                for (int n = 0; n < z[m].length; n++) {

                    // (z_i = z[m][n])
                    // sample from p(z_i|z_-i, w)
                	//核心步骤,通过论文中表达式(78)为文档m中的第n个词采样新的topic
                    int topic = sampleFullConditional(m, n);
                    z[m][n] = topic;
                }
            }

            // get statistics after burn-in
            //如果当前迭代轮数已经超过 burn-in的限制,并且正好达到 sample lag间隔
            //则当前的这个状态是要计入总的输出参数的,否则的话忽略当前状态,继续sample
            if ((i > BURN_IN) && (SAMPLE_LAG > 0) && (i % SAMPLE_LAG == 0)) {
                updateParams();
            }
        }
    }

    /**
     * Sample a topic z_i from the full conditional distribution: p(z_i = j |
     * z_-i, w) = (n_-i,j(w_i) + beta)/(n_-i,j(.) + W * beta) * (n_-i,j(d_i) +
     * alpha)/(n_-i,.(d_i) + K * alpha)
     * 
     * @param m
     *            document
     * @param n
     *            word
     */
    private int sampleFullConditional(int m, int n) {

        // remove z_i from the count variables
    	//这里首先要把原先的topic z(m,n)从当前状态中移除
        int topic = z[m][n];
        nw[documents[m][n]][topic]--;
        nd[m][topic]--;
        nwsum[topic]--;
        ndsum[m]--;

        // do multinomial sampling via cumulative method:
        double[] p = new double[K];
        for (int k = 0; k < K; k++) {
        	//nw 是第i个word被赋予第j个topic的个数
        	//在下式中,documents[m][n]是word id,k为第k个topic
        	//nd 为第m个文档中被赋予topic k的词的个数
            p[k] = (nw[documents[m][n]][k] + beta) / (nwsum[k] + V * beta)
                * (nd[m][k] + alpha) / (ndsum[m] + K * alpha);
        }
        // cumulate multinomial parameters
        for (int k = 1; k < p.length; k++) {
            p[k] += p[k - 1];
        }
        // scaled sample because of unnormalised p[]
        double u = Math.random() * p[K - 1];
        for (topic = 0; topic < p.length; topic++) {
            if (u < p[topic])
                break;
        }

        // add newly estimated z_i to count variables
        nw[documents[m][n]][topic]++;
        nd[m][topic]++;
        nwsum[topic]++;
        ndsum[m]++;

        return topic;
    }

    /**
     * Add to the statistics the values of theta and phi for the current state.
     */
    private void updateParams() {
        for (int m = 0; m < documents.length; m++) {
            for (int k = 0; k < K; k++) {
                thetasum[m][k] += (nd[m][k] + alpha) / (ndsum[m] + K * alpha);
            }
        }
        for (int k = 0; k < K; k++) {
            for (int w = 0; w < V; w++) {
                phisum[k][w] += (nw[w][k] + beta) / (nwsum[k] + V * beta);
            }
        }
        numstats++;
    }

    /**
     * Retrieve estimated document--topic associations. If sample lag > 0 then
     * the mean value of all sampled statistics for theta[][] is taken.
     * 
     * @return theta multinomial mixture of document topics (M x K)
     */
    public double[][] getTheta() {
        double[][] theta = new double[documents.length][K];

        if (SAMPLE_LAG > 0) {
            for (int m = 0; m < documents.length; m++) {
                for (int k = 0; k < K; k++) {
                    theta[m][k] = thetasum[m][k] / numstats;
                }
            }

        } else {
            for (int m = 0; m < documents.length; m++) {
                for (int k = 0; k < K; k++) {
                    theta[m][k] = (nd[m][k] + alpha) / (ndsum[m] + K * alpha);
                }
            }
        }

        return theta;
    }

    /**
     * Retrieve estimated topic--word associations. If sample lag > 0 then the
     * mean value of all sampled statistics for phi[][] is taken.
     * 
     * @return phi multinomial mixture of topic words (K x V)
     */
    public double[][] getPhi() {
        double[][] phi = new double[K][V];
        if (SAMPLE_LAG > 0) {
            for (int k = 0; k < K; k++) {
                for (int w = 0; w < V; w++) {
                    phi[k][w] = phisum[k][w] / numstats;
                }
            }
        } else {
            for (int k = 0; k < K; k++) {
                for (int w = 0; w < V; w++) {
                    phi[k][w] = (nw[w][k] + beta) / (nwsum[k] + V * beta);
                }
            }
        }
        return phi;
    }

    /**
     * Configure the gibbs sampler
     * 
     * @param iterations
     *            number of total iterations
     * @param burnIn
     *            number of burn-in iterations
     * @param thinInterval
     *            update statistics interval
     * @param sampleLag
     *            sample interval (-1 for just one sample at the end)
     */
    public void configure(int iterations, int burnIn, int thinInterval,
        int sampleLag) {
        ITERATIONS = iterations;
        BURN_IN = burnIn;
        THIN_INTERVAL = thinInterval;
        SAMPLE_LAG = sampleLag;
    }

    /**
     * Driver with example data.
     * 
     * @param args
     */
    public static void main(String[] args) {
        // words in documents
        int[][] documents = { {1, 4, 3, 2, 3, 1, 4, 3, 2, 3, 1, 4, 3, 2, 3, 6},
            {2, 2, 4, 2, 4, 2, 2, 2, 2, 4, 2, 2},
            {1, 6, 5, 6, 0, 1, 6, 5, 6, 0, 1, 6, 5, 6, 0, 0},
            {5, 6, 6, 2, 3, 3, 6, 5, 6, 2, 2, 6, 5, 6, 6, 6, 0},
            {2, 2, 4, 4, 4, 4, 1, 5, 5, 5, 5, 5, 5, 1, 1, 1, 1, 0},
            {5, 4, 2, 3, 4, 5, 6, 6, 5, 4, 3, 2}};
        // vocabulary
        int V = 7;
        int M = documents.length;
        // # topics
        int K = 2;
        // good values alpha = 2, beta = .5
        double alpha = 2;
        double beta = .5;

        LdaGibbsSampler lda = new LdaGibbsSampler(documents, V);
        
        //设定sample参数,采样运行10000轮,burn-in 2000轮,第三个参数没用,是为了显示
        //第四个参数是sample lag,这个很重要,因为马尔科夫链前后状态conditional dependent,所以要跳过几个采样
        lda.configure(10000, 2000, 100, 10);
        
        //跑一个!走起!
        lda.gibbs(K, alpha, beta);

        //输出模型参数,论文中式 (81)与(82)
        double[][] theta = lda.getTheta();
        double[][] phi = lda.getPhi();
    }
}
3月 04

Reading Note : Gibbs Sampling for the Uninitiated

 这几天较空总算把《Gibbs Sampling for the Uninitiated》看明白了点,看完这个对其他模型的 Gibbs Sampling 感觉应该是不会有大问题了,之后可以再仔细看看《Parameter estimation for text analysis》,然后对照他的代码写一份总结。

在此真是万分感谢这个作者,拯救了万千Uninitiated!要是没有这样的tutorial,也没人教,谁特么知道这东西怎么弄。

基础知识:随机过程、马尔科夫蒙特卡洛、Gibbs Sampling的理论基础(待写)

注:本文档主要是Learning Note,而不是独立教程,因此没有:1、基础概念  2、符号解释  3、详细公式推导  这些都可以从《Gibbs Sampling for the Uninitiated》原文找到。我的水平也极其有限,同时我认为如果只看别人博客而不愿意读原文那么肯定是不成的!本文主要目的是写一个学习记录给自己,使得以后忘记的时候能够迅速回顾。若能达到与人交流的目的,那就更好啦。

————————— 伟大的分割线:Gibbs Sampling for the Uninitiated!—————————

1.模型解析

首先应该看清楚这个model的真实面目,如下图:

对于文中Figure4,展开来应该如上图,之前我对于Graph Model其实具体接触比较少,对于Figure4那样画成Plate还是有误导性(主要是会误以为W只是一个随机变量),在这个图里Plate被展开了,所以我就随便用虚线替代实线,对于该图展开后的解释如下:

1.每一个文档有一个Label(j),是文档的class,同时$\theta_0$和$\theta_1$是和$L_j$对应的,如果$L_j=1$则对应的就是$\theta_1$

2.在Gibbs Sampling中,所有的圆圈都是要被sample的,也就是每一个$L_j$都会被算作一个$z_i$(在Gibbs Sampling表达式里面$z_i^{(t+1)}\sim P(Z_i|z_1^{(t+1)},…,z_{i-1}^{(t+1)},z_{i+1}^{(t)},…)$中的随机变量z),同理,每一个$w_{jk}$也是一个$z_i$,只不过在具体sampling的时候,$L_j$是会被单独sample的,而$w_{jk}$是整体sample的,这个会在后面讲到。

3.在这个model中,Gibbs Sampling所谓的P(Z),就是产生图中这整个数据集的联合概率,也就是产生这N个文档整体联合概率,还要算上包括超参$\gamma$产生具体$\pi$和 $\theta$的概率。所以最后得到了上图中表达式与对应色彩。

2.联合概率公式推导

分别对红色、绿色、蓝色、紫色部分求表达式,得到如下结果:

1)红色部分:这个是从beta分布sample出一个伯努利分布,伯努利分布只有一个参数就是$\pi$,不要normalization项(要求的是整个联合概率,所以在这里纠结normalization是没有用的),得到:

$P(\pi|\gamma_{\pi1},\gamma_{\pi0})\varpropto \pi^{\gamma_{\pi1}-1}(1-\pi)^{\gamma_{\pi0}-1}$

2)绿色部分:这里L是一整个向量,其中值为0的有C0个,值为1的有C1个,多次伯努利分布就是二项分布啦,因此:

$P(L|\pi)\varpropto \pi^{C_1}(1-\pi)^{C_0}$

3)蓝色部分:对于0类和1类的两个$\theta$都采样自参数为$\gamma_\theta$的狄利克雷分布,注意所有这些都是向量,有V维,每一维度对应一个Word。根据狄利克雷的PDF得到以下表达式,其实这个表达式有两个,分别为$\theta_0$和$\theta_1$用相同的式子采样:

$P(\theta|\gamma_{\theta})\varpropto \prod\limits_{i=1}^V \theta_i^{\gamma_{\theta i}-1}$

4)紫色部分:这部分,首先要求对于单独一个文档n,产生所有word也就是$W_n$的概率。假设对于某个文档,$\theta=(0.2,0.5,0.3)$,意思就是word1产生概率0.2,word2产生概率0.5,假如这个文档里word1有2个,word2有3个,word3有2个,则这个文档的产生概率就是(0.2*0.2)*(0.5*0.5*0.5)*(0.3*0.3)。所以按照这个道理,一个文档整个联合概率如下:

$P(W_n|L,\theta_{L_n})=\prod\limits_{i=1}^V \theta_i^{W_{ni}}$

上面这个概率是针对单个文档而言的,把所有文档的这些概率乘起来,就得到了紫色部分:

$P(\mathbb{C}_x|L,\theta_x)=\prod\limits_{n\in \mathbb{C}_x}\prod\limits_{i=1}^V \theta_{x,i}^{W_{ni}}=\prod\limits_{i=1}^V \theta_{x,i}^{N_{C_x}(i)}$

其中x的取值可以是0或1,所以$\mathbb{C}_x$可以是$\mathbb{C}_0$或$\mathbb{C}_1$,当x=0时,$n\in\mathbb{C}_x$的意思就是对于所有class 0中的文档,然后对于这些文档中每一个word i,word i在该文档中出现$W_{ni}次$,求$\theta_{0,i}$的$W_{ni}$次方,所有这些乘起来就是紫色部分。后面那个表达式是规约后得到的结果,$N_{C_0}(i)$的意思就是word i出现在calss为0的所有文档中的总数,同理对于$N_{C_1}(i)$。

最后把上述4部分乘起来,就得到整个模型的联合概率。不再赘述,只是说明文中从式(29)跳到(30),以及式(31)到式(32)的计算,都用到了Beta分布与二项分布共轭的性质,以及多项分布与狄利克雷分布共轭性质,这些可以从PRML2.2和2.3中读到【1】。

3.将隐含变量$\pi$积出

为了方便,可以对隐含变量$\pi$进行积分,最后达到消去这个变量的目的。

图中带圈的(1)作者是直接通过Beta分布的Normalize Constant得出的,也就是和【1】是同一个道理。

4.吉布斯采样框架

包含三个步骤,按照顺序(1),(2),(3)进行

(1)对所有L变量采样一轮。即整个数据集合有N个文档,也就是有N个L实例,依次对每一个文档j,采样$L_j$:

$P(L_j)^{(t+1)}=P(L_j|L_1^{(t+1)},…,L_{j-1}^{(t+1)},L_{j+1}^{(t)},…L_N^{(t)},\mathbb{C}^{(t)},\theta_1^{(t)},\theta_0^{(t)};\mu)$

注意表达式后面那部分的样子,在条件概率中,没有$L_j$项,同时$L_j$之前的项是(t+1),就是基于已经采样过的$L_1 \to L_{j-1}$的结果,而$L_j$之后的所有项都是基于(t)轮,也就是前一轮采样的结果。

(2)对$\theta_0$采样,可以发现下式中是没有$\theta_0$的,这就是Gibbs采样的要求:

$P(\theta_0)^{(t+1)}=P(\theta_0^{(t+1)}|L_1^{(t+1)},L_2^{(t+1)},…,L_N^{(t+1)},\mathbb{C},\theta_1^{(t)};\mu)$

(3)对$\theta_1$采样:

$P(\theta_1)^{(t+1)}=P(\theta_1^{(t+1)}|L_1^{(t+1)},L_2^{(t+1)},…,L_N^{(t+1)},\mathbb{C},\theta_0^{(t+1)};\mu)$

对于其他变量,$\gamma$为超参数,是人为设定的,W为观察变量,不需采样。

然后就要求出$P(L_j^{(t+1)}|L^{(-j)},\mathbb{C}^{(-j)},\theta_0,\theta_1;\mu)$的表达式,这里有一些关键的理解,首先求出这个表达式之后要怎样对$L_j$进行采样。事实上他的具体做法是这样的,有了$P(L_j)^{(t+1)}$的表达式,我们要把$L_j=0$和$L_j=1$两个情况带入表达式,分别求出L为0和1的概率,这样就可以为$L_j^{(t+1)}$挑选一个新的值。怎么挑这个值,举例来说,如果$P(L_j^{(t+1)}=0)=0.2$,$P(L_j^{(t+1)}=1)=0.3$,那么就可以产生一个伯努利分布,其参数:

$\mu=P(L_j)^{(t+1)}被赋值为0的概率=P(L_j=0)/(P(L_j=0)+P(L_j=1))=0.4$

然后我用random函数从uniform distribution产生一个随机数r,如果r<=0.4,则将$L_j^{(t+1)}$赋值为0,否则赋值为1。

在这里不累赘写L的概率表达式了,记录一些关键内容。首先,计算$P(L_j^{(t+1)})$的时候,是要查看$L_j^{(t)}$的值的,因为我们要把文档$L_j$从原来的集合中去掉,那么如果$L_j^{(t)}==1$,就要从C1中扣除1,否则要从C0中扣除1。那么式(45)在计算的时候,为什么分子分母会被越成这个样子?因为对于特定的$L_j^{(t)}$等于0或1代入时(也就是Lj上一个状态的值),C0和C1只有其中一个会被减去1,而另外一个不变,不变化的那个就被分子分母约掉了!所以(45)中只剩下$C_x$,也就是被减了1的那个C。同理,对于式(47),看起来好像是$\theta_0$和$\theta_1$乘在一起,其实只有其中一个会被减去L中的词的个数,另外一个不变的就全部约掉了,所以最后只会剩下式(48)。

另外,对于式(52),其实作者还是没有写清楚,怎么从伽马分布$Gamma(\alpha_i,1)$采样yi呢?呵呵,这个正好在PRML11.1.2 Reject Sampling里面有的哦~

最后的最后,文章给出了很关键的理解:

​$P(L_j=x|L^{(-j)},\mathbb{C}^{(-j)},\theta_0,\theta_1;\mu)=\frac{C_x+\gamma_{\pi x}}{N+\gamma_{\pi1}+\gamma_{\pi2}-1}\prod\limits_{i=1}^V \theta_{x,i}^{W_{ji}}$

这个式子中,前半部分其实只有Cx是变量,所以如果C0大,则$P(L_j^{(t+1)}=0)$的概率就会大一点,所以下一次Lj的值就会倾向于C0,反之就会倾向于C1。而后半部分呢,是在判断当前$\theta$参数的情况下,整个文档的likelihood更倾向于C0还是C1。