# Latent Dirichlet allocation¶

Latent Dirichlet allocation is a widely used topic model. The data is a collection of documents which contain words. The goal of the analysis is to find topics (distribution of words in topics) and document topics (distribution of topics in documents).

## Data¶

The data consists of two vectors of equal length. The elements in these vectors correspond to the words in all documents combined. If there were documents and each document had words, the vectors contain elements. Let be the number of documents in total. The first vector gives each word a document index defining to which document the word belongs. Let be the size of the whole available vocabulary. The second vector gives each word a vocabulary index defining which word it is from the vocabulary.

For this demo, we will just generate an artificial dataset for simplicity. We use the LDA model itself to generate the dataset. First, import relevant packages:

>>> import numpy as np
>>> from bayespy import nodes


Let us decide the number of documents and the number of words in those documents:

>>> n_documents = 10
>>> n_words = 10000


Randomly choose into which document each word belongs to:

>>> word_documents = nodes.Categorical(np.ones(n_documents)/n_documents,
...                                    plates=(n_words,)).random()


Let us also decide the size of our vocabulary:

>>> n_vocabulary = 100


Also, let us decide the true number of topics:

>>> n_topics = 5


Generate some random distributions for the topics in each document:

>>> p_topic = nodes.Dirichlet(1e-1*np.ones(n_topics),
...                           plates=(n_documents,)).random()


Generate some random distributions for the words in each topic:

>>> p_word = nodes.Dirichlet(1e-1*np.ones(n_vocabulary),
...                          plates=(n_topics,)).random()


Sample topic assignments for each word in each document:

>>> topic = nodes.Categorical(p_topic[word_documents],
...                           plates=(n_words,)).random()


And finally, draw vocabulary indices for each word in all the documents:

>>> corpus = nodes.Categorical(p_word[topic],
...                            plates=(n_words,)).random()


Now, our dataset consists of word_documents and corpus, which define the document and vocabulary indices for each word in our dataset.

Todo

Use some large real-world dataset, for instance, Wikipedia.

## Model¶

Variable for learning the topic distribution for each document:

>>> p_topic = nodes.Dirichlet(np.ones(n_topics),
...                           plates=(n_documents,),
...                           name='p_topic')


Variable for learning the word distribution for each topic:

>>> p_word = nodes.Dirichlet(np.ones(n_vocabulary),
...                          plates=(n_topics,),
...                          name='p_word')


The document indices for each word in the corpus:

>>> from bayespy.inference.vmp.nodes.categorical import CategoricalMoments
>>> document_indices = nodes.Constant(CategoricalMoments(n_documents), word_documents,
...                                   name='document_indices')


Variable for learning the topic assignments of each word in the corpus:

>>> topics = nodes.Categorical(nodes.Gate(document_indices, p_topic),
...                            plates=(len(corpus),),
...                            name='topics')


The vocabulary indices for each word in the corpus:

>>> words = nodes.Categorical(nodes.Gate(topics, p_word),
...                           name='words')


## Inference¶

Observe the corpus:

>>> words.observe(corpus)


Break symmetry by random initialization:

>>> p_topic.initialize_from_random()
>>> p_word.initialize_from_random()


Construct inference engine:

>>> from bayespy.inference import VB
>>> Q = VB(words, topics, p_word, p_topic, document_indices)


Run the VB learning algorithm:

>>> Q.update(repeat=1000)
Iteration ...


## Results¶

Use bayespy.plot to plot the results:

>>> import bayespy.plot as bpplt


Plot the topic distributions for each document:

>>> bpplt.pyplot.figure()
<matplotlib.figure.Figure object at 0x...>
>>> bpplt.hinton(Q['p_topic'])
>>> bpplt.pyplot.title("Posterior topic distribution for each document")
<matplotlib.text.Text object at 0x...>
>>> bpplt.pyplot.xlabel("Topics")
<matplotlib.text.Text object at 0x...>
>>> bpplt.pyplot.ylabel("Documents")
<matplotlib.text.Text object at 0x...>


Plot the word distributions for each topic:

>>> bpplt.pyplot.figure()
<matplotlib.figure.Figure object at 0x...>
>>> bpplt.hinton(Q['p_word'])
>>> bpplt.pyplot.title("Posterior word distributions for each topic")
<matplotlib.text.Text object at 0x...>
>>> bpplt.pyplot.xlabel("Words")
<matplotlib.text.Text object at 0x...>
>>> bpplt.pyplot.ylabel("Topics")
<matplotlib.text.Text object at 0x...>


Todo

Create more illustrative plots.

## Stochastic variational inference¶

LDA is a popular example for stochastic variational inference (SVI). Using SVI for LDA is quite simple in BayesPy. In SVI, only a subset of the dataset is used at each iteration step but this subset is “repeated” to get the same size as the original dataset. Let us define a size for the subset:

>>> subset_size = 1000


Thus, our subset will be repeat this many times:

>>> plates_multiplier = n_words / subset_size


Note that this multiplier doesn’t need to be an integer.

Now, let us repeat the model construction with only one minor addition. The following variables are identical to previous:

>>> p_topic = nodes.Dirichlet(np.ones(n_topics),
...                           plates=(n_documents,),
...                           name='p_topic')
>>> p_word = nodes.Dirichlet(np.ones(n_vocabulary),
...                          plates=(n_topics,),
...                          name='p_word')


The document indices vector is now a bit shorter, using only a subset:

>>> document_indices = nodes.Constant(CategoricalMoments(n_documents),
...                                   word_documents[:subset_size],
...                                   name='document_indices')


Note that at this point, it doesn’t matter which elements we chose for the subset. For the topic assignments of each word in the corpus we need to use plates_multiplier because these topic assignments for the subset are “repeated” to recover the full dataset:

>>> topics = nodes.Categorical(nodes.Gate(document_indices, p_topic),
...                            plates=(subset_size,),
...                            plates_multiplier=(plates_multiplier,),
...                            name='topics')


Finally, the vocabulary indices for each word in the corpus are constructed as before:

>>> words = nodes.Categorical(nodes.Gate(topics, p_word),
...                           name='words')


This node inherits the plates and multipliers from its parent topics, so there is no need to define them here. Again, break symmetry by random initialization:

>>> p_topic.initialize_from_random()
>>> p_word.initialize_from_random()


Construct inference engine:

>>> from bayespy.inference import VB
>>> Q = VB(words, topics, p_word, p_topic, document_indices)


In order to use SVI, we need to disable some lower bound checks, because the lower bound doesn’t anymore necessarily increase at each iteration step:

>>> Q.ignore_bound_checks = True


For the stochastic gradient ascent, we’ll define some learning parameters:

>>> delay = 1
>>> forgetting_rate = 0.7


Run the inference:

>>> for n in range(1000):
...     # Observe a random mini-batch
...     subset = np.random.choice(n_words, subset_size)
...     Q['words'].observe(corpus[subset])
...     Q['document_indices'].set_value(word_documents[subset])
...     # Learn intermediate variables
...     Q.update('topics')
...     # Set step length
...     step = (n + delay) ** (-forgetting_rate)
...     # Stochastic gradient for the global variables

>>> bpplt.pyplot.figure()