Wenjie Fan     04/08/2023

Markov chain Monte Carlo (MCMC) is a sampling method widely used for computation of hierarchical models in Bayesian statistics. One of the examples, Gibbs sampler, is simple compared to others like Metropolis-Hastings algorithm when the full conditionals (i.e. one parameter conditioned on others) probabilities can be explicitly evaluated. In Latent Dirichlet Allocation (LDA), a collapsed Gibbs sampler (CGS) is used where one marginalize out the topic distribution parameter \(\theta\) and the topic-word distribution parameter \(\phi\) in [1] (or \(\theta, \beta\) in [2]), leaving the topic assignments \(\mathbf{z}\) to be sampled (and to be fully conditioned on).

Load AP corpora

We use the same dataset (AP corpus) tested in LDA with VEM, and maintain the notations:

Considering the algorithmic structure of CGS, the way we store the documents is different from the form word id :occurrence in the case of VEM. Instead, we directly store the words with repetitions in each document, i.e. length(corpus$docs[[d]]$words) will be equal to corpus$docs[[d]]$N. Stopwords provided by stopwords-iso are also removed from each document in the corpus.

# Read file
file <- file('data/ap.dat', open='r')
M <- 0 # Number of documents
V <- 0 # Size of vocabulary, set to be the max representing number of the words
docs <- list()

# List of stopwords
SW <- stopwords("en", source = "stopwords-iso")
addSW <- c('percent','monday','tuesday','wednesday','thursday','friday','saturday','sunday')
vocab = read.table('data/vocab.txt', stringsAsFactors=FALSE)
SW <- which(vocab$V1 %in% SW) - 1

# Iterate through the documents
while (length(doc <- readLines(file, n=1, warn=FALSE)) != 0) {
    M <- M + 1
    doc <- unlist(strsplit(doc, ' '))
    nwords <- as.integer(doc[1]) # Number of distinct words
    words <- c()
    occur <- c()
    sw.count <- 0
    sw.total <- 0
    N <- 0 # Word count of document
    # Iterate through the words
    for (n in 1:nwords) {
        word.occur <- unlist(strsplit(doc[n+1], ':'))
        if (word.occur[1] %in% SW) { # If the word is a stop word
            sw.count <- sw.count + 1
            sw.total <- sw.total + as.integer(word.occur[2])
        }
        else { # Or, the word is not a stop word
            words <- c(words, rep(as.integer(word.occur[1]), as.integer(word.occur[2])))
            occur <- c(occur, as.integer(word.occur[2]))
            N <- N + occur[length(occur)]
            if (words[length(words)] > V) {V <- words[length(words)]}
        }
    }
    nwords <- nwords - sw.count
    docs[[M]] <- list(words=words, occur=occur, nwords=nwords, N=N)
}
close(file)

# Generate corpus: list of documents as sublists
V <- V + 1 # Representing number starts with 0
corpus <- list(docs=docs, M=M, V=V)

The runtime for the gibbs sampler could be long (each sample approximately \(2\) minutes). For experimental purposes, we could eliminate the corpus to half by randomly removing the documents. The following chunk will NOT be compiled otherwise (eval=FALSE).

removing <- sample(1:M, as.integer(M * 0.75), replace=FALSE) # Halving the corpus
corpus$docs <- corpus$docs[-removing]
corpus$M <- M - as.integer(M * 0.75)
M <- corpus$M

Collapsed Gibbs sampler

The structure of LDA model has been described in “Latent Dirichlet Allocation (1)” so we only present some essential specifications here. The Dirichlet prior for each column of the topic-word matrix \(\phi\) is just the “smoothing” technique in [1]. Notice that here we are not rigorously considering \(w_n\) as a vector in \(\mathbb{R}^V\). For readability, for example, by \(\phi_{iw_n}\) we are referring to the \((i,j)\)-th of the matrix if the word \(w_n\) has index \(j\) in the vocabulary (i.e. \(w_n = \text{the j-th vector in canonical basis}\)). \[ \text{For} ~~ \phi_i \in \mathbb{R}^V, ~i=1,...,K: \quad \quad p(\phi_{i}) = \frac{\Gamma(V\beta)}{\Gamma(\beta)^V} \prod_{v=1}^V \phi_{iv}^{\beta-1}, \quad \quad p(w_n \mid z_n, \phi_i) = \phi_{iw_n}, ~\forall n \] \[ \text{For} ~~ \theta_d \in \mathbb{R}^K, ~d=1,...,M: \quad \quad p(\theta_d) = \frac{\Gamma(K\alpha)}{\Gamma(\alpha)^K} \prod_{i=1}^K \theta_{di}^{\alpha-1}, \quad \quad p(z_n \mid \theta_d) = \theta_{dz_n},~ \forall n \]

Thus following these definitions, we can obtain \(p(\mathbf{z}\mid\mathbf{w})\) by \(p(\mathbf{w}\mid\mathbf{z}) p(\mathbf{z})\). Denote \(n_i^{(v)}\) the number of times the word \(v \in \{1,...,V\}\) was assigned to topic \(i\), then by argumenting the Dirichlet distribution we have: \[ p(\mathbf{w}\mid\mathbf{z}) = \left(\frac{\Gamma(V\beta)}{\Gamma(\beta)^V}\right)^K \prod_{i=1}^K \int \prod_{v=1}^V \phi_{iv}^{\beta -1 + n_i^{(v)}} d\phi_i = \left(\frac{\Gamma(V\beta)}{\Gamma(\beta)^V}\right)^K \prod_{i=1}^K \frac{\prod_{v=1}^V \Gamma(n_i^{(v)} + \beta)}{\Gamma(n_i^{(\cdot)} + V\beta)} \] Also denote \(n_i^{(d)}\) the number of times some word from the document \(d\in \{1,...,M\}\) was assigned to topic \(i\), similarly we have: \[ p(\mathbf{z}) = \left(\frac{\Gamma(K\alpha)}{\Gamma(\alpha)^K}\right)^M \prod_{d=1}^M \int \prod_{i=1}^K \theta_{di}^{\alpha -1 + n_i^{(d)}} d\theta_d = \left(\frac{\Gamma(K\alpha)}{\Gamma(\alpha)^K}\right)^M \prod_{d=1}^M \frac{\prod_{i=1}^K \Gamma(n_i^{(d)} + \alpha)}{\Gamma(n_{~\cdot}^{(d)} + K\alpha)} \]

As mentioned by [2], the posterior \(p(\mathbf{z}\mid\mathbf{w})\) can be computed by \(\frac{p(\mathbf{w}, \mathbf{z})}{\sum_\mathbf{z} p(\mathbf{w}, \mathbf{z})}\) but the denominator is intractable. Gibbs sampling solves this problem.

Full conditional probabilities

The core of Gibbs sampling is the derivation of full conditional probabilities \(P(z_n = i \mid \mathbf{z}_{-n}, \mathbf{w})\) for each topic \(i\in \{1,...,K\}\). Through some manipulation, we can represent this in the form of two probabilities we derived above: \[ P(z_n = i \mid \mathbf{z}_{-n}, \mathbf{w}) = \frac{p(\mathbf{z}, \mathbf{w})}{p(\mathbf{z}_{-n}, \mathbf{w})} = \frac{p(\mathbf{z}) p(\mathbf{w}\mid\mathbf{z})}{p(\mathbf{z}_{-n}) p(\mathbf{w}_{-n}\mid \mathbf{z}_{-n}) p(w_n)} \propto \frac{p(\mathbf{z}) p(\mathbf{w}\mid\mathbf{z})}{p(\mathbf{z}_{-n}) p(\mathbf{w}_{-n}\mid \mathbf{z}_{-n})} \]

It suffices to calculate the product of two terms: \(\frac{p(\mathbf{z})}{p(\mathbf{z}_{-n})}\) and \(\frac{p(\mathbf{w}\mid\mathbf{z})}{p(\mathbf{w}_{-n}\mid \mathbf{z}_{-n})}\) for each topic \(i\). Since the product will then be considered as the weights for sampling the next state (topic) of the current word in MCMC. To simplify and derive an explicit expression for these two terms, we need the identity \(\Gamma(x+1) = x\Gamma(x)\).

Start with the first term. There should be \(\prod_{d=1}^M\) and \(\prod_{i=1}^K\) in the ratio below. However, if we assume the current word \(w_n\) is from document \(\tilde{d}\), and has been assigned to topic \(\tilde{i}\), then terms involving documents other than \(\tilde{d}\) and topics other than \(\tilde{i}\) will be eliminated. Hence the ratio can be simplified as follows: \[ \frac{p(\mathbf{z})}{p(\mathbf{z}_{-n})} = \frac{\Gamma(n_\tilde{i}^{(\tilde{d})} + \alpha)}{\Gamma(n_{~\cdot}^{(\tilde{d})} + K\alpha)} \cdot \frac{\Gamma(n_{~\cdot, -n}^{(\tilde{d})} + K\alpha)}{ \Gamma(n_{\tilde{i},-n}^{(\tilde{d})} + \alpha)} = \frac{n_{\tilde{i},-n}^{(\tilde{d})} + \alpha}{n_{~\cdot,-n}^{(\tilde{d})} + K\alpha} \]

The second term is similar. Assuming the current word is \(\tilde{v}\) in the vocabulary \(\{1,...,V\}\): \[ \frac{p(\mathbf{w}\mid\mathbf{z})}{p(\mathbf{w}_{-n}\mid \mathbf{z}_{-n})} = \frac{\Gamma(n_\tilde{i}^{(\tilde{v})} + \beta)}{\Gamma(n_{\tilde{i}}^{(\cdot)} + V\beta)} \cdot \frac{\Gamma(n_{\tilde{i}, -n}^{(\cdot)} + V\beta)}{ \Gamma(n_{\tilde{i},-n}^{(\tilde{v})} + \beta)} = \frac{n_{\tilde{i},-n}^{(\tilde{v})} + \beta}{n_{\tilde{i},-n}^{(\cdot)} + V\beta} \]

These derivation makes the definition of the full conditional probability in the Gibbs sampler concise in two lines:

# This is Eq [5] in Griffiths & Steyvers, 2004
# 'sample' function in R receive unnormalized weights
first.term <- (n.di[d, ] + alpha) / (sum(n.di[d,]) + K * alpha) 
second.term <- (n.wi[w_dn, ] + beta) / (colSums(n.wi) + V * beta)

In order the check the convergence of the Gibbs sampler, we check the ‘perplexity’ index used in [1]. The following chunk is the function for computing perplexity given the word-topic and document-topic counts of the current iteration.

perplexity <- function(n.wi, n.di) {
    ll <- 0
    for (i in 1:K) {
        ll <- ll + lgamma(V*beta) - V * lgamma(beta)
        ll <- ll + sum(lgamma(n.wi[,i] + beta))
        ll <- ll - lgamma(sum(n.wi[,i]) + V*beta)
    }
    for (d in 1:M) {
        ll <- ll + lgamma(K*alpha) - K * lgamma(alpha)
        ll <- ll + sum(lgamma(n.di[d,] + alpha))
        ll <- ll - lgamma(sum(n.di[d,]) + K*alpha)
    }
    n.tokens <- sum(sapply(corpus$docs, function(x) x$N))
    return(exp(- ll / n.tokens))
}

We follow the popular setting of \(\alpha=50/K\) and \(\beta=0.1\). The initial states are assigned randomly. Note that rigorously one should keep a track on the states of each iteration, so that afterwards the mean or mode will be taken for the parameter estimation. However, this involves computation using tensors or large arrays in R. For the sampling process to be carried out in a reasonable amount of time, we will use a single state after the burn-in period of MCMC. And actually this still produce great and accurate results in topic modelling, as we shall see in the model interpretation output.

# MCMC setup
set.seed(123)
K <- 30
N.iter <- 100 # Number of MCMC samples
beta <- 0.1 
alpha <- 50/K

states <- lapply(corpus$docs, function(x) rep(0, x$N)) # Zeros according to doc lengths
n.wi <- matrix(0, V, K) # word-topic matrix
n.di <- matrix(0, M, K) # document-topic matrix

# Initial state 
# Iterate through documents
for (d in 1:M) { 
    doc <- corpus$docs[[d]]
    N <- doc$N
    # Iterate through words
    for (n in 1:N) {
        w_dn <- doc$words[n] + 1 # Account for word number starting from 0
        # States initialization
        states[[d]][n] <- sample(1:K, 1)
        # Counts initialization
        i.alloc <- states[[d]][n]
        # Count of word w being assigned to topic i throughout
        n.wi[w_dn, i.alloc] <- n.wi[w_dn, i.alloc] + 1
    }
    # Count of some words from document d being assigned to topic i 
    for (i in 1:K) { n.di[d, i] <- sum(states[[d]] == i) }    
}
# Initialise perplexity
perp <- c(perplexity(n.wi,n.di))

# Sampler main loop
for (t in 1:N.iter) {
    # print(sprintf('Iteration: %d', t))
    # Iterate through documents
    for (d in 1:M) {
        doc <- corpus$docs[[d]]
        N <- doc$N
        # Iterate through words
        for (n in 1:N) {
            w_dn <- doc$words[n] + 1
            
            # Leave-current-word-out
            i.old <- states[[d]][n] # From previous iteration
            n.wi[w_dn, i.old] <- n.wi[w_dn, i.old] - 1
            n.di[d, i.old] <- n.di[d, i.old] - 1
            
            # MCMC walk
            term1 <- (n.wi[w_dn, ] + beta) / (colSums(n.wi) + V * beta)
            term2 <- (n.di[d, ] + alpha) / (sum(n.di[d,]) + K * alpha) 
            cond.prob.weight <- term1 * term2
            i.new <- sample(1:K, 1, prob=cond.prob.weight)
            
            # Update state
            n.wi[w_dn, i.new] <- n.wi[w_dn, i.new] + 1
            n.di[d, i.new] <- n.di[d, i.new] + 1
            states[[d]][n] <- i.new
        }
    }
    # Record perplexity of current iteration
    perp <- c(perp, perplexity(n.wi,n.di))
}

We can examine the topics by looking at the top words (the words with the highest probability given the topic).

# Parameter estimation
phi <- t(t(n.wi + beta) / (colSums(n.wi) + V * beta)) # Size V*K, the transposes are important
theta <- (n.di + alpha) / (rowSums(n.di) + K * alpha) # Size M*K

# Topic with highest probability for each document
topic <- apply(theta, 1, which.max)

# Top words in each topic
top.n <- 20 # Top 20 words each topic being shown
word.id = matrix(0, top.n, K)
probs = matrix(0, top.n, K)
for (i in 1:K) {
  topic.top = sort(phi[,i], decreasing=TRUE, index.return=TRUE)
  word.id[,i] = topic.top$ix[1:top.n]
  probs[,i] = topic.top$x[1:top.n]
}
topic.top.words = matrix(vocab[word.id], top.n, K)
# First row
kable(topic.top.words[,1:10], booktabs = TRUE) %>%
    column_spec(1:10, border_left = T, border_right = T) %>%
    add_header_above(paste0("Topic ", 1:10)) %>% 
    add_header_above(c("Topic top words, 1~10" = 10)) %>%
    kable_styling()
Topic top words, 1~10
Topic 1
Topic 2
Topic 3
Topic 4
Topic 5
Topic 6
Topic 7
Topic 8
Topic 9
Topic 10
market john air health miles company news government people oil
stock film officials medical water corp time military program cents
york record plane hospital southern business television president tax prices
exchange york flight people city offer tuesday defense programs farmers
trading book force aids north announced monday human taxes food
issues magazine aircraft drug coast share week rights cut futures
stocks died eastern study northern chairman network north social trade
wall art spokesman dr san companies wednesday security income cent
prices movie airport care reported executive expected force service farm
shares war navy disease damage president calls national country production
rose editor airlines treatment river vice angeles rebels director tons
jones vietnam base center county chief series report housing market
investors martin accident doctors california deal broadcast tuesday private agriculture
bond newspaper ship heart park companys cbs reagan cities price
average history planes blood central firm los dec security corn
week william service patients texas sale spokesman country poor supplies
volume director pilots cancer homes management nbc sept land july
board wrote airline researchers fair sell scheduled june homeless gasoline
securities published travel children officials agreement thursday secret public products
street george report virus people industry sunday political store drought 
# Second row
kable(topic.top.words[,11:20], booktabs = TRUE) %>%
    column_spec(1:10, border_left = T, border_right = T) %>%
    add_header_above(paste0("Topic ", 11:20)) %>% 
    add_header_above(c("Topic top words, 11~20" = 10)) %>%
    kable_styling()
Topic top words, 11~20
Topic 11
Topic 12
Topic 13
Topic 14
Topic 15
Topic 16
Topic 17
Topic 18
Topic 19
Topic 20
united court people bush money soviet family school south space
bush decision government dukakis federal party children students government project
soviet rights israel campaign bank government yearold university united building
east law israeli democratic plan political wife black africa mission
president ruling police republican financial gorbachev day city african launch
west aid killed president pay union time women national shuttle
agreement supreme capital jackson claims communist life college president scientists
germany judge peace election debt leader son white north day
german ruled china presidential fund opposition mother education nations time
reagan appeals protest vice banks minister days student black earth
europe civil jewish george loan elections parents schools panama nasa
summit lawyers army convention funds power woman board korea feet
administration legal official candidates savings prime child public noriega center
aid request minister support estate moscow husband community korean planned
trade commission violence vote assets country daughter mayor country field
economic nicaragua security voters trust parties left job political director
union contra demonstrators candidate insurance leaders living hotel fighting venus
nations abortion reported democrats loans central father professor mandela built
foreign courts palestinian primary cash economic months teachers remains start
baker appeal west political week gorbachevs apartment citys people ground 
# Third row
kable(topic.top.words[,21:30], booktabs = TRUE) %>%
    column_spec(1:10, border_left = T, border_right = T) %>%
    add_header_above(paste0("Topic ", 21:30)) %>% 
    add_header_above(c("Topic top words, 21~30" = 10)) %>%
    kable_styling()
Topic top words, 21~30
Topic 21
Topic 22
Topic 23
Topic 24
Topic 25
Topic 26
Topic 27
Topic 28
Topic 29
Topic 30
late house meeting prison workers percent military department people police
dollar committee officials trial union sales war office american people
japan congress told drug plant rate iraq federal america killed
japanese senate official charges contract increase troops attorney play city
yen rep news court environmental report forces investigation lot shot
london budget conference convicted strike economy kuwait district coming authorities
thursday chairman talks death labor months iraqi filed british died
gold sen president charged officials month united justice radio night
trade legislation friday judge jobs economic gulf government tour church
major administration reporters attorney employees reported iran law brown death
bid congressional week guilty plants department saudi documents played car
friday vote agency charge safety rose american hearing music told
close measure met accused department growth soldiers judge audience officers
west policy press arrested waste rates army washington life victims
tokyo approved time cocaine pay june force suit looked saturday
british president foreign jury industry total attack service fall hospital
rates action day crime company quarter held fraud stage injured
french lawmakers washington murder local annual arabia indictment religious shooting
dealers proposal visit sentence energy survey hostages lawsuit hands yearold
dollars defense efforts arrest cars increased crisis assistant playing attack

References

[1]
D. M. Blei, A. Y. Ng, and M. I. Jordan, “Latent dirichlet allocation,” Journal of machine Learning research, vol. 3, no. Jan, pp. 993–1022, 2003.
[2]
T. L. Griffiths and M. Steyvers, “Finding scientific topics,” Proceedings of the National academy of Sciences, vol. 101, no. suppl_1, pp. 5228–5235, 2004.