I am writing a pymc3-based implementation of Latent Dirichlet Allocation, and am referencing this CrossValidated answer (modified for pymc3) as well as pymc3's own tutorial on LDA, in addition to the Wikipedia article on LDA.
But I am hitting a problem conceptually.
In the linked CrossValidated answer above, the final distribution for all of the word-position-specific Categorical (multinomial) outputs looks like this:
[
pm.Categorical(
"w_%i_%i" % (d, i),
p=pm.Lambda(
'phi_z_%i_%i' % (d, i),
lambda z=z[d][i], phi=phi: phi[z]
),
value=data[d][i],
observed=True
)
for d in range(D) for i in range(Wd[d])
]
where z
and phi
are defined just as in the Wikipedia formulation of LDA, data[d]
is a vector of word counts (bag of words vector) for document number d
, D
is the number of documents, and Wd[d]
is the document length of document d
.
You don't need to dig around for what every variable name means. The problem I have is that the value
parameter is set to data[d][i]
, for this observed=True
random variable -- meaning whatever is in data[d][i]
will define how the likelihood is penalized for matching the observed value, in tun affecting the gradient for choosing proposals for samples to draw, etc.
There are two approaches you could take to feeding the data into the model. First, which seems standard in both the sklearn LDA and the pymc3 'official' tutorial linked above, is to use a standard bag-of-words representation, so that data[d][i]
would be the count of word i
appearing in document d
.
But the nested for
loop of the model above shows that there is a Categorical parameter for each (document, position) pair, and not for each (document, word index) pair.
In fact, if a given document happened to have more word positions in it than there are total words in the vocabulary, then the indexing data[d][i]
might even cause an indexing error, because the data are counts for each word, and are not indexed by the position within the document.
This would be more in line with the second way of feeding in data: as a sequence of word index values, where the sequence has an entry for each position in a document, so its length is the length of the document.
In more mathematical terms, this is the set of distributional relationships defining the graphical model for LDA:
$$\boldsymbol\phi_{k=1 \dots K} \sim \operatorname{Dirichlet}_V(\boldsymbol\beta)\\ \boldsymbol\theta_{d=1 \dots M} \sim \operatorname{Dirichlet}_K(\boldsymbol\alpha)\\ z_{d=1 \dots M,w=1 \dots N_d} \sim \operatorname{Categorical}_K(\boldsymbol\theta_d) \\ w_{d=1 \dots M,w=1 \dots N_d} \sim \operatorname{Categorical}_V(\boldsymbol\phi_{z_{dw}}) $$
and the parameter $N_{d}$ is the length of document $d$. So in this setting, the parameters clearly should correspond to all of the word positions, and not to the word counts.
Over in the linked pymc3 tutorial, a custom log-likelihood function is written that accounts for this, by multiplying the count by the log likelihood on the other conditional terms.
But then this makes it seem that pymc3's tutorial conflicts with the graphical model from Wikipedia and from the linked CrossValidated answer: those both assume that the output prediction is a Categorical selection of a word index for each position in each document.
I hope it's clear what the conflict is. It looks to me like a problem where you "can't have it both ways."
You cannot both input the data in a bag-of-words format (where count values are used as the observed values in the likelihood function), and also have a generative model for each (document, position) word location.
You could either use count data as the input and then generate predictions about the occurrence of each word in the vocabulary as output (like unscaled bag-of-words vectors)...
Or, you could provide the observed word sequences (not counts), and then have a generative model that matches the Wikipedia and CrossValidated links, predicting a word for each (document, index) position.
If my either-or statement is correct, then does this mean the code from linked CrossValidated answer is erroneous when it uses value=data[d][i]
(the word count of vocbulary position i
in document d
) as the observational data?
Under that CrossValidated code, shouldn't the observed value
parameter be the word index for the word that actually occurs in document d
at position i
? Feeding in bag-of-words data to that particular observational model is wrong?
data[d][i]
has to be wrong. I think instead, you have to usedata[d][vocab_index[original_document[d][i]]]
, wherevocab_index
is a hash table mapping a word to an index in the vocabulary, andoriginal_document[d]
is the actual sequence of words (not the counts as are stored indata
). $\endgroup$