<jupyter_start><jupyter_text>3-1. Latent Dirichlet Allocation1&nbsp;&nbsp;Data: Reuters News Titles2&nbsp;&nbsp;Model: Basic LDA2.1&nbsp;&nbsp;Variational EM2.1.1&nbsp;&nbsp;E-step2.1.2&nbsp;&nbsp;M-step2.1.3&nbsp;&nbsp;Variational lower bound2.1.4&nbsp;&nbsp;Training3&nbsp;&nbsp;Model: Smoothed LDA3.1&nbsp;&nbsp;Collapsed Gibbs Sampling3.1.1&nbsp;&nbsp;Run Gibbs sampler3.1.2&nbsp;&nbsp;Recover $\beta$ and $\theta$ from the sample3.2&nbsp;&nbsp;Variational EM3.2.1&nbsp;&nbsp;E-step3.2.2&nbsp;&nbsp;M-step  Data: Reuters News Titles Reuters is a multi-class, multi-label dataset.* 90 classes* 10788 documents    * 7769 training documents    * 3019 testing documents<jupyter_code>from nltk.corpus import reuters
from nltk.corpus import stopwords<jupyter_output><empty_output><jupyter_text>* train-test split: The data is already splitted. Just sort it out.<jupyter_code>stops = stopwords.words("english")
stops += [
    "a", "about", "above", "across", "after", "afterwards", "again", "against",
    "all", "almost", "alone", "along", "already", "also", "although", "always",
    "am", "among", "amongst", "amoungst", "amount", "an", "and", "another",
    "any", "anyhow", "anyone", "anything", "anyway", "anywhere", "are",
    "around", "as", "at", "back", "be", "became", "because", "become",
    "becomes", "becoming", "been", "before", "beforehand", "behind", "being",
    "below", "beside", "besides", "between", "beyond", "bill", "both",
    "bottom", "but", "by", "call", "can", "cannot", "cant", "co", "con",
    "could", "couldnt", "cry", "de", "describe", "detail", "do", "done",
    "down", "due", "during", "each", "eg", "eight", "either", "eleven", "else",
    "elsewhere", "empty", "enough", "etc", "even", "ever", "every", "everyone",
    "everything", "everywhere", "except", "few", "fifteen", "fifty", "fill",
    "find", "fire", "first", "five", "for", "former", "formerly", "forty",
    "found", "four", "from", "front", "full", "further", "get", "give", "go",
    "had", "has", "hasnt", "have", "he", "hence", "her", "here", "hereafter",
    "hereby", "herein", "hereupon", "hers", "herself", "him", "himself", "his",
    "how", "however", "hundred", "i", "ie", "if", "in", "inc", "indeed",
    "interest", "into", "is", "it", "its", "itself", "keep", "last", "latter",
    "latterly", "least", "less", "ltd", "made", "many", "may", "me",
    "meanwhile", "might", "mill", "mine", "more", "moreover", "most", "mostly",
    "move", "much", "must", "my", "myself", "name", "namely", "neither",
    "never", "nevertheless", "next", "nine", "no", "nobody", "none", "noone",
    "nor", "not", "nothing", "now", "nowhere", "of", "off", "often", "on",
    "once", "one", "only", "onto", "or", "other", "others", "otherwise", "our",
    "ours", "ourselves", "out", "over", "own", "part", "per", "perhaps",
    "please", "put", "rather", "re", "same", "see", "seem", "seemed",
    "seeming", "seems", "serious", "several", "she", "should", "show", "side",
    "since", "sincere", "six", "sixty", "so", "some", "somehow", "someone",
    "something", "sometime", "sometimes", "somewhere", "still", "such",
    "system", "take", "ten", "than", "that", "the", "their", "them",
    "themselves", "then", "thence", "there", "thereafter", "thereby",
    "therefore", "therein", "thereupon", "these", "they", "thick", "thin",
    "third", "this", "those", "though", "three", "through", "throughout",
    "thru", "thus", "to", "together", "too", "top", "toward", "towards",
    "twelve", "twenty", "two", "un", "under", "until", "up", "upon", "us",
    "very", "via", "was", "we", "well", "were", "what", "whatever", "when",
    "whence", "whenever", "where", "whereafter", "whereas", "whereby",
    "wherein", "whereupon", "wherever", "whether", "which", "while", "whither",
    "who", "whoever", "whole", "whom", "whose", "why", "will", "with",
    "within", "without", "would", "yet", "you", "your", "yours", "yourself",
    "yourselves", ".", "!", "?", ",", ";", ":", "[", "]", "{", "}", "-", "+", 
    "_", "/", "@", "#", "$", "%", "^", "&", "*", "(", ")", "<", ">", "|", "=",
    ".-", ".,", "'", '"', ',"'
]
reuters.words()[:10]
# fetch titles only
# fetch 2000 docs only
trainset, testset = [], []
vocab = []

i = 0
for file_id in reuters.fileids():
    if file_id.startswith("train"):
        doc = [w.lower() for w in reuters.words(file_id) \
                 if (w.isupper()) \
                 if (w.lower() not in stops) \
                 and (not w.isnumeric())]
        if doc:
            trainset.append(doc)
            vocab += doc
            i += 1
    else:
        testset.append([w.lower() for w in reuters.words(file_id) \
                         if (w.isupper()) \
                         and (w.lower() not in stops) \
                         and (not w.isnumeric())])
    if i >= 2000:
        break
vocab = list(set(vocab))
word_to_ix = {w: i for i, w in enumerate(vocab)}
def seq_to_ix(seq, vocab=vocab):
    # len(vocab), which is the last index, is for the <unk> (unknown) token
    unk_idx = len(vocab)
    return np.array(list(map(lambda w: word_to_ix.get(w, unk_idx), seq)))

data = {
    "train": list(map(seq_to_ix, trainset)),
    "test": list(map(seq_to_ix, testset))
}
data["train"][0][:5]  # word indices<jupyter_output><empty_output><jupyter_text>Model: Basic LDA  Variational EM For each document $\mathbf{w}$ in a corpus $D$, generate$$N \sim \mathcal{P}(\xi) \\\theta \sim \text{Dir}(\alpha)$$and for $n = 1, \cdots, N$, generate$$z_n \sim \text{Multi}(\theta) \\w_n \sim P(w_n | z_n, \beta)$$where $\beta \in \mathbb{R}^{k \times V}$, $\beta_{ij} = P(w^j = 1| z^i = 1)$. * $\alpha, \beta$: hyperparameters (Dirichlet, Multinomial).* $N$: The number of words in the document. (ancillary variable)* $\theta$: A topic mixture.* (For $i$ in $1\cdots N$)  * $z_n$: A topic variable.  * $w_n$: A generated word.<jupyter_code>def init_lda(docs, vocab, n_topic, gibbs=False, random_state=0):
    if gibbs:
        global V, k, N, M, alpha, eta, n_iw, n_di
    else:
        global V, k, N, M, alpha, beta, gamma, phi
        
    np.random.seed(random_state)

    V = len(vocab)
    k = n_topic  # number of topics
    N = np.array([doc.shape[0] for doc in docs])
    M = len(docs)

    print(f"V: {V}\nk: {k}\nN: {N[:10]}...\nM: {M}")

    # initialize α, β
    
    if gibbs:
        alpha = np.random.gamma(shape=100, scale=0.01, size=1)  # one for all k
        eta = np.random.gamma(shape=100, scale=0.01, size=1)  # one for all V
        print(f"α: {alpha}\nη: {eta}")
        
        n_iw = np.zeros((k, V), dtype=int)
        n_di = np.zeros((M, k), dtype=int)
        print(f"n_iw: dim {n_iw.shape}\nn_di: dim {n_di.shape}")
    else:
        alpha = np.random.gamma(shape=100, scale=0.01, size=k) #np.random.rand(k)
        beta = np.random.dirichlet(np.ones(V), k)
        print(f"α: dim {alpha.shape}\nβ: dim {beta.shape}")

        # initialize ϕ, γ
        ## ϕ: (M x max(N) x k) arrays with zero paddings on the right
        gamma = alpha + np.ones((M, k)) * N.reshape(-1, 1) / k

        phi = np.ones((M, max(N), k)) / k
        for m, N_d in enumerate(N):
            phi[m, N_d:, :] = 0  # zero padding for vectorized operations

        print(f"γ: dim {gamma.shape}\nϕ: dim ({len(phi)}, N_d, {phi[0].shape[1]})")<jupyter_output><empty_output><jupyter_text>E-step Let $\phi_d \in \mathbb{R}^{N \times k}, \gamma_d \in \mathbb{R}^k$ be variational parameters for $\alpha, \beta$. For a document $\mathbf{w}_d$, $d = 1,\cdots,M$,1. initialize $\phi_{dni}^0 := 1/k$ for all $i,n$.2. initialize $\gamma_{di}^0 := \alpha_i + N/k$ for all $i$.3. **repeat until** convergence    1. for $n=1$ to $N$        1. for $i=1$ to $k$            1. $\phi_{dni}^{t+1} := \beta_{i\mathbf{w}_{dn}}\exp\left(\Psi(\gamma_{di}^t) - \Psi(\sum_{j=1}^k \gamma_{dj}^t)\right)$        2. normalize $\phi_{dn}^{t+1}$ to sum to 1    2. $\gamma_d^{t+1} := \alpha + \sum_{n=1}^N \phi_{dn}^{t+1}$    where $\Psi$ is the first derivative of the $\log\Gamma$ function.<jupyter_code>def E_step(docs, phi, gamma, alpha, beta):
    """
    Minorize the joint likelihood function via variational inference.
    This is the E-step of variational EM algorithm for LDA.
    """
    # optimize phi
    for m in range(M):
        phi[m, :N[m], :] = (beta[:, docs[m]] * np.exp(psi(gamma[m, :]) - psi(gamma[m, :].sum())).reshape(-1, 1)).T

        # Normalize phi
        phi[m, :N[m]] /= phi[m, :N[m]].sum(axis=1).reshape(-1, 1)
        if np.any(np.isnan(phi)):
            raise ValueError("phi nan")

    # optimize gamma
    gamma = alpha + phi.sum(axis=1)

    return phi, gamma<jupyter_output><empty_output><jupyter_text>M-step $$\beta_{ij} \propto \sum_{d=1}^M \sum_{n=1}^N \phi_{dni} \mathbf{w}_{dn}^j$$ $\alpha$ is updated via Newton-Raphson method: $$\frac{\partial L}{\partial \alpha_i}   = M\left( \Psi\left(\sum_{j=1}^k \alpha_j\right) - \Psi(\alpha_i) \right)    - \sum_{d=1}^M \left( \Psi(\gamma_{di}) - \Psi\left(\sum_{j=1}^k \gamma_{dj}\right) \right) \\\frac{\partial^2 L}{\partial \alpha_i \alpha_j} = M \left( \Psi'\left(\sum_{j=1}^k \alpha_j\right) - \delta(i,j) \Psi'(\alpha_i) \right)$$where $\delta(i,j) = 1$ if $i=j$, $0$ otherwise.<jupyter_code>def M_step(docs, phi, gamma, alpha, beta, M):
    """
    maximize the lower bound of the likelihood.
    This is the M-step of variational EM algorithm for (smoothed) LDA.
    
    update of alpha follows from appendix A.2 of Blei et al., 2003.
    """
    # update alpha
    alpha = _update(alpha, gamma, M)
    
    # update beta
    for j in range(V):
        beta[:, j] = np.array([_phi_dot_w(docs, phi, m, j) for m in range(M)]).sum(axis=0)
    beta /= beta.sum(axis=1).reshape(-1, 1)

    return alpha, beta
import warnings

def _update(var, vi_var, const, max_iter=10000, tol=1e-6):
    """
    From appendix A.2 of Blei et al., 2003.
    For hessian with shape `H = diag(h) + 1z1'`
    
    To update alpha, input var=alpha and vi_var=gamma, const=M.
    To update eta, input var=eta and vi_var=lambda, const=k.
    """
    for _ in range(max_iter):
        # store old value
        var0 = var.copy()
        
        # g: gradient 
        psi_sum = psi(vi_var.sum(axis=1)).reshape(-1, 1)
        g = const * (psi(var.sum()) - psi(var)) \
            + (psi(vi_var) - psi_sum).sum(axis=0)

        # H = diag(h) + 1z1'
        z = const * polygamma(1, var.sum())  # z: Hessian constant component
        h = -const * polygamma(1, var)       # h: Hessian diagonal component
        c = (g / h).sum() / (1./z + (1./h).sum())

        # update var
        var -= (g - c) / h
        
        # check convergence
        err = np.sqrt(np.mean((var - var0) ** 2))
        crit = err < tol
        if crit:
            break
    else:
        warnings.warn(f"max_iter={max_iter} reached: values might not be optimal.")
    
    #print(err)
    return var

def _phi_dot_w(docs, phi, d, j):
    """
    \sum_{n=1}^{N_d} ϕ_{dni} w_{dn}^j
    """
    # doc = np.zeros(docs[m].shape[0] * V, dtype=int)
    # doc[np.arange(0, docs[m].shape[0] * V, V) + docs[m]] = 1
    # doc = doc.reshape(-1, V)
    # lam += phi[m, :N[m], :].T @ doc
    return (docs[d] == j) @ phi[m, :N[d], :]<jupyter_output><empty_output><jupyter_text>Variational lower bound $$\begin{align}L(\gamma, \phi; \alpha, \beta)  &= \log\Gamma(\sum_{j=1}^k \alpha_j) - \sum_{i=1}^k \log\Gamma(\alpha_i) + \sum_{i=1}^k (\alpha_i - 1) \left(\Psi(\gamma_i) - \Psi(\sum_{i=1}^k \gamma_i)\right) \\  &+ \sum_{n=1}^N \sum_{i=1}^k \phi_{ni} \left(\Psi(\gamma_i) - \Psi(\sum_{i=1}^k \gamma_i)\right) \\  &+ \sum_{n=1}^N \sum_{i=1}^k \sum_{j=1}^V \phi_{ni} \mathbf{w}_{n}^j \log\beta_{ij} \\  &- \log\Gamma(\sum_{i=1}^k \gamma_i) + \sum_{i=1}^k \log\Gamma(\gamma_i) - \sum_{i=1}^k (\gamma_i - 1) \left(\Psi(\gamma_i) - \Psi(\sum_{i=1}^k \gamma_i)\right) \\  &- \sum_{n=1}^N \sum_{i=1}^k \phi_{ni} \log\phi_{ni}\end{align}$$<jupyter_code>def dg(gamma, d, i):
    """
    E[log θ_t] where θ_t ~ Dir(gamma)
    """
    return psi(gamma[d, i]) - psi(np.sum(gamma[d, :]))


def dl(lam, i, w_n):
    """
    E[log β_t] where β_t ~ Dir(lam)
    """
    return psi(lam[i, w_n]) - psi(np.sum(lam[i, :]))

def vlb(docs, phi, gamma, alpha, beta, M, N, k):
    """
    Average variational lower bound for joint log likelihood.
    """
    lb = 0
    for d in range(M):
        lb += (
            gammaln(np.sum(alpha))
            - np.sum(gammaln(alpha))
            + np.sum([(alpha[i] - 1) * dg(gamma, d, i) for i in range(k)])
        )

        lb -= (
            gammaln(np.sum(gamma[d, :]))
            - np.sum(gammaln(gamma[d, :]))
            + np.sum([(gamma[d, i] - 1) * dg(gamma, d, i) for i in range(k)])
        )

        for n in range(N[d]):
            w_n = int(docs[d][n])

            lb += np.sum([phi[d][n, i] * dg(gamma, d, i) for i in range(k)])
            lb += np.sum([phi[d][n, i] * np.log(beta[i, w_n]) for i in range(k)])
            lb -= np.sum([phi[d][n, i] * np.log(phi[d][n, i]) for i in range(k)])

    return lb / M<jupyter_output><empty_output><jupyter_text>Training * Only on 2,000 documents<jupyter_code>import os
# for beutiful plot later, reorder training set
# (don't need to do this)

if "lda_trainset.idx" in os.listdir():
    with open("lda_trainset.idx") as r:
        idx = eval(r.read())
    docs = np.array(data["train"])[idx].tolist()
else:
    docs = data["train"]
init_lda(docs, vocab, n_topic=10)
%%time

N_EPOCH = 1000
TOL = 0.1

verbose = True
lb = -np.inf

for epoch in range(N_EPOCH): 
    # store old value
    lb_old = lb
    
    # Variational EM
    phi, gamma = E_step(docs, phi, gamma, alpha, beta)
    alpha, beta = M_step(docs, phi, gamma, alpha, beta, M)
    
    # check anomaly
    if np.any(np.isnan(alpha)):
        print("NaN detected: alpha")
        break
    
    # check convergence
    lb = vlb(docs, phi, gamma, alpha, beta, M, N, k)
    err = abs(lb - lb_old)
    
    # check anomaly
    if np.isnan(lb):
        print("NaN detected: lb")
        break
        
    if verbose:
        print(f"{epoch: 04}:  variational_lb: {lb: .3f},  error: {err: .3f}")
    
    if err < TOL:
        break
else:
    warnings.warn(f"max_iter reached: values might not be optimal.")

print(" ========== TRAINING FINISHED ==========")<jupyter_output>000:  variational_lb: -63.016,  error:  inf
 001:  variational_lb: -62.239,  error:  0.777
 002:  variational_lb: -61.789,  error:  0.450
 003:  variational_lb: -61.427,  error:  0.362
 004:  variational_lb: -61.108,  error:  0.319
 005:  variational_lb: -60.816,  error:  0.292
 006:  variational_lb: -60.545,  error:  0.270
 007:  variational_lb: -60.285,  error:  0.261
 008:  variational_lb: -60.018,  error:  0.267
 009:  variational_lb: -59.734,  error:  0.284
 010:  variational_lb: -59.426,  error:  0.308
 011:  variational_lb: -59.092,  error:  0.334
 012:  variational_lb: -58.737,  error:  0.355
 013:  variational_lb: -58.366,  error:  0.371
 014:  variational_lb: -57.988,  error:  0.379
 015:  variational_lb: -57.614,  error:  0.374
 016:  variational_lb: -57.259,  error:  0.355
 017:  variational_lb: -56.931,  error:  0.328
 018:  variational_lb: -56.632,  error:  0.299
 019:  variational_lb: -56.362,  error:  0.270
 020:  variational_lb: -56.120,  error:  0.242
 021:  variatio[...]<jupyter_text>* Training result 1. Topic extraction<jupyter_code>def n_most_important(beta_i, n=30):
    """
    find the index of the largest `n` values in a list
    """
    
    max_values = beta_i.argsort()[-n:][::-1]
    return np.array(vocab)[max_values]
for i in range(k):
    print(f"TOPIC {i:02}: {n_most_important(beta[i], 9)}")<jupyter_output>TOPIC 00: ['ec' 'csr' 'loss' 'bank' 'icco' 'unit' 'cocoa' 'fob' 'petroleum']
TOPIC 01: ['raises' 'acquisition' 'prices' 'prime' 'rate' 'w' 'completes' 'mar'
 'imports']
TOPIC 02: ['year' 'sets' 'sees' 'net' 'stock' 'dividend' 'l' 'industries' 'corp']
TOPIC 03: ['pct' 'cts' 'gdp' 'opec' 'shr' 'february' 'plc' 'ups' 'rose']
TOPIC 04: ['u' 'japan' 'fed' 'trade' 'gaf' 'ems' 'says' 'dlr' 'gnp']
TOPIC 05: ['usda' 'f' 'ag' 'international' 'sell' 'report' 'ge' 'corn' 'wheat']
TOPIC 06: ['k' 'market' 'money' 'rate' 'eep' 'treasury' 'prime' 'says' 'mln']
TOPIC 07: ['qtr' 'note' '4th' 'net' 'loss' 'ico' 'corp' '1st' 'group']
TOPIC 08: ['dlrs' 'mln' 'bp' 'march' 'corp' 'week' 'canada' 'e' 'bank']
TOPIC 09: ['unit' 'buy' 'says' 'sale' 'c' 'sells' 'completes' 'american' 'grain']<jupyter_text>2. Topic-word & document-document distribution<jupyter_code>n_sample = 10000
theta_hat = np.array([np.random.dirichlet(gamma[d], n_sample).mean(0) for d in range(M)])
theta_hat /= theta_hat.sum(1).reshape(-1, 1)
plt.figure(figsize=(8,8))
plt.subplot(121)
n_plot_words = 150
sns.heatmap(beta.T[:n_plot_words], xticklabels=[], yticklabels=[])
plt.xlabel("Topics", fontsize=14)
plt.ylabel(f"Words[:{n_plot_words}]", fontsize=14)
plt.title("topic-word distribution", fontsize=16)

plt.subplot(122)
sns.heatmap(theta_hat, xticklabels=[], yticklabels=[])
plt.xlabel("Topics", fontsize=14)
plt.ylabel("Documents", fontsize=14)
plt.title("document-topic distribution", fontsize=16)

plt.tight_layout();<jupyter_output><empty_output><jupyter_text>---  Model: Smoothed LDA For each document $\mathbf{w}$ in a corpus $D$, generate$$N \sim \mathcal{P}(\xi) \\\beta \sim \text{Dir}(\lambda) \\\theta \sim \text{Dir}(\alpha)$$and for $n = 1, \cdots, N$, generate$$z_n \sim \text{Multi}(\theta) \\w_n \sim P(w_n | z_n, \beta)$$where $\beta \in \mathbb{R}^{k \times V}$, $\beta_{ij} = P(w^j = 1| z^i = 1)$. * $\alpha, \eta$: Dirichlet hyperparameters.* $\beta$: Unsmoothed multinomial hyperparameter.* $N$: The number of words in the document. (ancillary variable)* $\theta$: A topic mixture.* (For $i$ in $1\cdots N$)  * $z_n$: A topic variable.  * $w_n$: A generated word.  Collapsed Gibbs Sampling $$P(z_{dn}^i=1 | \mathbf{z}_{(-dn)},\mathbf{w}) \propto \frac{n_{(-dn),iw_{dn}}+\eta}{n_{(-dn),i\cdot}+V\eta} \frac{n_{(-dn),dj}+\alpha}{n_{(-dn),d\cdot}+k\alpha}.$$  Run Gibbs sampler<jupyter_code>def _init_gibbs(docs, vocab, n_topic, n_gibbs=2000):
    """
    Initialize t=0 state for Gibbs sampling.
    Replace initial word-topic assignment ndarray (M, N, N_GIBBS) in-place.
    """
    # initialize variables
    init_lda(docs, vocab, n_topic=n_topic, gibbs=True)
    
    # word-topic assignment
    global assign
    N_max = max(N)
    assign = np.zeros((M, N_max, n_gibbs+1), dtype=int)
    print(f"assign: dim {assign.shape}")
    
    # initial assignment
    for d in range(M):
        for n in range(N[d]):
            # randomly assign topic to word w_{dn}
            w_dn = docs[d][n]
            assign[d, n, 0] = np.random.randint(k)

            # increment counters
            i = assign[d, n, 0]
            n_iw[i, w_dn] += 1
            n_di[d, i] += 1
def _conditional_prob(w_dn, d):
    """
    P(z_{dn}^i=1 | z_{(-dn)}, w)
    """
    prob = np.empty(k)
    
    for i in range(k):
        # P(w_dn | z_i)
        _1 = (n_iw[i, w_dn] + eta) / (n_iw[i, :].sum() + V*eta)
        # P(z_i | d)
        _2 = (n_di[d, i] + alpha) / (n_di[d, :].sum() + k*alpha)
        
        prob[i] = _1 * _2
    
    return prob / prob.sum()
def run_gibbs(docs, vocab, n_topic, n_gibbs=2000, verbose=True):
    """
    Run collapsed Gibbs sampling
    """
    # initialize required variables
    _init_gibbs(docs, vocab, n_topic, n_gibbs)
    
    if verbose:
        print("\n", "="*10, "START SAMPLER", "="*10)
    
    # run the sampler
    for t in range(n_gibbs):
        for d in range(M):
            for n in range(N[d]):
                w_dn = docs[d][n]
                
                # decrement counters
                i_t = assign[d, n, t]  # previous assignment
                n_iw[i_t, w_dn] -= 1
                n_di[d, i_t] -= 1

                # assign new topics
                prob = _conditional_prob(w_dn, d)
                i_tp1 = np.argmax(np.random.multinomial(1, prob))

                # increment counter according to new assignment
                n_iw[i_tp1, w_dn] += 1
                n_di[d, i_tp1] += 1
                assign[d, n, t+1] = i_tp1
        
        # print out status
        if verbose & ((t+1) % 50 == 0):
            print(f"Sampled {t+1}/{n_gibbs}")
run_gibbs(docs, vocab, n_topic=10, n_gibbs=2000)<jupyter_output>V: 4048
k: 10
N: [16  5  6  7  9  7  8  8  6  5]...
M: 2000
α: [1.18335424]
η: [1.03715177]
n_iw: dim (10, 4048)
n_di: dim (2000, 10)
assign: dim (2000, 88, 2001)

 ========== START SAMPLER ==========
Sampled 50/2000
Sampled 100/2000
Sampled 150/2000
Sampled 200/2000
Sampled 250/2000
Sampled 300/2000
Sampled 350/2000
Sampled 400/2000
Sampled 450/2000
Sampled 500/2000
Sampled 550/2000
Sampled 600/2000
Sampled 650/2000
Sampled 700/2000
Sampled 750/2000
Sampled 800/2000
Sampled 850/2000
Sampled 900/2000
Sampled 950/2000
Sampled 1000/2000
Sampled 1050/2000
Sampled 1100/2000
Sampled 1150/2000
Sampled 1200/2000
Sampled 1250/2000
Sampled 1300/2000
Sampled 1350/2000
Sampled 1400/2000
Sampled 1450/2000
Sampled 1500/2000
Sampled 1550/2000
Sampled 1600/2000
Sampled 1650/2000
Sampled 1700/2000
Sampled 1750/2000
Sampled 1800/2000
Sampled 1850/2000
Sampled 1900/2000
Sampled 1950/2000
Sampled 2000/2000<jupyter_text>Recover $\beta$ and $\theta$ from the sample<jupyter_code>β̂ = np.empty((k, V))
θ̂ = np.empty((M, k))

for j in range(V):
    for i in range(k):
        β̂[i, j] = (n_iw[i, j] + eta) / (n_iw[i, :].sum() + V*eta)

for d in range(M):
    for i in range(k):
        θ̂[d, i] = (n_di[d, i] + alpha) / (n_di[d, :].sum() + k*alpha)
plt.figure(figsize=(8,8))
plt.subplot(121)
n_plot_words = 150
sns.heatmap(β̂.T[:n_plot_words], xticklabels=[], yticklabels=[])
plt.xlabel("Topics", fontsize=14)
plt.ylabel(f"Words[:{n_plot_words}]", fontsize=14)
plt.title("topic-word distribution", fontsize=16)

plt.subplot(122)
sns.heatmap(θ̂, xticklabels=[], yticklabels=[])
plt.xlabel("Topics", fontsize=14)
plt.ylabel("Documents", fontsize=14)
plt.title("document-topic distribution", fontsize=16)

plt.tight_layout();<jupyter_output><empty_output><jupyter_text>Variational EM  E-step Let $\phi_d \in \mathbb{R}^{N \times k}, \gamma_d \in \mathbb{R}^k, \lambda \in \mathbb{R}^{k \times V}$ be variational parameters for $\alpha, \beta, \eta$.  Suppose further that for $\beta \in \mathbb{R}^{k \times V}$, $\beta_i^0 \sim \text{Dir}(\lambda^0)$ where $\lambda_i^0 = \eta$ for all $i$. For a document $\mathbf{w}_d$, $d = 1,\cdots,M$,1. initialize $\phi_{dni}^0 := 1/k$ for all $i,n$.2. initialize $\gamma_{di} := \alpha_i + N/k$ for all $i$.3. **repeat until** convergence    1. for $n=1$ to $N$        1. for $i=1$ to $k$            1. $\phi_{dni}^{t+1} := \exp\left(\Psi(\lambda_{iw_{dn}}^t) - \Psi(\sum_{j=1}^V \lambda_{ij}^t) + \Psi(\gamma_{di}^t) - \Psi(\sum_{j=1}^k \gamma_{dj}^t)\right)$            1. for $j=1$ to $V$                1. $\lambda_{ij} = \eta + \sum_{d=1}^M \sum_{n=1}^{N_d} \phi_{dni} w_{dn}^j$        2. normalize $\phi_{dn}^{t+1}$ to sum to 1    2. $\gamma_d^{t+1} := \alpha + \sum_{n=1}^N \phi_{dn}^{t+1}$    where $\Psi$ is the first derivative of the $\log\Gamma$ function.  M-step $$\beta_{ij} \propto \sum_{d=1}^M \sum_{n=1}^N \phi_{dni} \mathbf{w}_{dn}^j$$ $\alpha$ is updated via Newton-Raphson method: $$\frac{\partial L}{\partial \alpha_i}   = M\left( \Psi\left(\sum_{j=1}^k \alpha_j\right) - \Psi(\alpha_i) \right)    - \sum_{d=1}^M \left( \Psi(\gamma_{di}) - \Psi\left(\sum_{j=1}^k \gamma_{dj}\right) \right) \\\frac{\partial^2 L}{\partial \alpha_i \alpha_j} = M \left( \Psi'\left(\sum_{j=1}^k \alpha_j\right) - \delta(i,j) \Psi'(\alpha_i) \right)$$where $\delta(i,j) = 1$ if $i=j$, $0$ otherwise.<jupyter_code>class SmoothedLDA:
    
    def __init__(self, docs, vocab, k):
        self.docs = docs
        
        self.V = len(vocab)
        self.k = k  # number of topics
        self.N = np.array([doc.shape[0] for doc in docs])
        self.M = len(docs)
        
        V = self.V
        N = self.N
        M = self.M
        
        # initialize model parameters
        ##self.beta = np.ones((k, V)) / V
        self.alpha = np.random.gamma(100, 0.01, k)
        self.eta = np.ones(V)

        # initialize variational parameters
        # ϕ: (M x max(N) x k) arrays with zero paddings on the right
        self.phi = [np.ones((N[d], k)) / k for d in range(M)]
        self.gamma = alpha + (N / k).reshape(-1, 1)
        self.lam = np.random.gamma(shape=100, scale=0.01, size=(k, V))
        
        
    def _update_phi(self):
        """
        Update variational parameter phi
        ϕ_{n, j} ∝ e^[ (Ψ(λ_j) - Ψ(Σλ_j)) + ( Ψ(γ_j) - Ψ(Σγ_j) ) ]
        """
        M = self.M
        N = self.N
        k = self.k

        phi = self.phi
        gamma = self.gamma
        docs = self.docs

        for d in range(M):
            for n in range(N[d]):
                for i in range(k):
                    w_n = int(docs[d][n])
                    phi[d][n, i] = np.exp(dl(lam, i, w_n) + dg(gamma, d, i))

                # Normalize over topics
                phi[d][n, :] = phi[d][n, :] / np.sum(phi[d][n, :])
                
        return phi
    
    def _update_gamma(self):
        """
        Update variational parameter gamma
        γ_t = α_t + Σ_{n=1}^{N_d} ϕ_{t, n}
        """
        M = self.M
        phi = self.phi
        alpha = self.alpha

        gamma = alpha + np.array(
            list(map(lambda x: x.sum(axis=0), phi))
        )
        
        return gamma
    
    
    def _update_lam(self):
        V = self.V
        N = self.N
        M = self.M
        phi = self.phi
        lam = self.lam
        eta = self.eta
        docs = self.docs
        
        lam[:] = eta
        for d in range(M):  #, desc="MINORIZE lam"):
            doc = np.zeros(N[d] * V, dtype=int)
            doc[np.arange(0, N[d] * V, V) + docs[d]] = 1
            doc = doc.reshape(-1, V)
            
            lam += phi[d].T @ doc
        
        return lam
        
    
    def _update_alpha(self, max_iter=1000, tol=0.1):
        """
        Update alpha with linear time Newton-Raphson.
        """
        M = self.M
        k = self.k

        alpha = self.alpha
        gamma = self.gamma

        for _ in range(max_iter):
            alpha_old = alpha

            #  Calculate gradient
            g = M * (psi(np.sum(alpha)) - psi(alpha)) +\
                (psi(gamma) - psi(np.sum(gamma, axis=1)).reshape(-1, 1)).sum(axis=0)

            #  Calculate Hessian diagonal component
            h = -M * polygamma(1, alpha)

            #  Calculate Hessian constant component
            z = M * polygamma(1, np.sum(alpha))

            #  Calculate constant
            c = np.sum(g / h) / (z ** (-1.0) + np.sum(h ** (-1.0)))

            #  Update alpha
            alpha = alpha - (g - c) / h
            
            #  Check convergence
            if np.sqrt(np.mean(np.square(alpha - alpha_old))) < tol:
                break
        else:
            warnings.warn("_update_alpha(): max_iter reached.")

        return alpha
    
    
    def _update_eta(self, max_iter=1000, tol=0.1):
        """
        Update eta with linear time Newton-Raphson.
        """
        M = self.M
        k = self.k

        eta = self.eta
        lam = self.lam

        for _ in range(max_iter):
            eta_old = eta

            #  Calculate gradient
            g = k * (psi(np.sum(eta)) - psi(eta)) +\
                (psi(lam) - psi(np.sum(lam, axis=1)).reshape(-1, 1)).sum(axis=0)

            #  Calculate Hessian diagonal component
            h = -k * polygamma(1, eta)

            #  Calculate Hessian constant component
            z = k * polygamma(1, np.sum(eta))

            #  Calculate constant
            c = np.sum(g / h) / (z ** (-1.0) + np.sum(h ** (-1.0)))

            #  Update alpha
            eta = eta - (g - c) / h

            #  Check convergence
            if np.sqrt(np.mean(np.square(eta - eta_old))) < tol:
                break
        else:
            warnings.warn("_update_eta(): max_iter reached.")

        return eta
    
    
    def _E_step(self):
        """
        E-step of the variational EM algorithm.
        Update ϕ, γ, λ.
        """
        self.phi = self._update_phi()
        self.gamma = self._update_gamma()
        self.lam = self._update_lam()
        
        
    def _M_step(self):
        """
        M-step of the variational EM algorithm.
        Update α, η.
        """
        self.alpha = self._update_alpha()
        self.eta = self._update_eta()
    
    
    def vlb(self):
        """
        lower bound from variational inference
        """
        phi = self.phi
        gamma = self.gamma
        lam = self.lam
        alpha = self.alpha
        eta = self.eta
        docs = self.docs
        
        M = self.M
        k = self.k
        N = self.N
        
        a0, a1, a2, a3_1, a3_2, a4, a5 = 0., 0., 0., 0., 0., 0., 0.
        for d in range(M):
            a0 += (
                k * (
                    gammaln(np.sum(eta)) 
                    - np.sum(gammaln(eta))
                )
                + np.sum([(eta[j] - 1) * dl(lam, i, j) for j in range(V) for i in range(k)])
            )
            a1 += (
                gammaln(np.sum(alpha))
                - np.sum(gammaln(alpha))
                + np.sum([(alpha[i] - 1) * dg(gamma, d, i) for i in range(k)])
            )

            a4 += (
                gammaln(np.sum(gamma[d, :]))
                - np.sum(gammaln(gamma[d, :]))
                + np.sum([(gamma[d, i] - 1) * dg(gamma, d, i) for i in range(k)])
            )
            
            for i in range(k):
                for j in range(V):
                    a3_2 += (
                        gammaln(np.sum(lam[i, j]))
                        - np.sum(gammaln(lam[i, :]))
                        + np.sum((lam[i, j] - 1) * dl(lam, i, j))
                    )

            for n in range(N[d]):
                w_n = int(docs[d][n])
                a2 += np.sum([phi[d][n, i] * dg(gamma, d, i) for i in range(k)])
                a3_1 += np.sum([phi[d][n, i] * dl(lam, i, w_n) for i in range(k)])
                a5 += np.sum([phi[d][n, i] * np.log(phi[d][n, i]) for i in range(k)])

        return a0 + a1 + a2 + a3_1 - a3_2 - a4 - a5
    
    
    def train(self, max_iter=1000, tol=5, verbose=True):
        vlb = -np.inf
        
        for it in range(max_iter):
            old_vlb = vlb
            self._E_step()
            self._M_step()
            
            vlb = self.vlb()
            err = vlb - old_vlb
            
            if verbose:
                print(f"Iteration {it+1}: {vlb: .3f} (delta: {err: .2f})") 
            
            if err < tol:
                break
        else:
            warnings.warn("max_iter reached.")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.19261522459858285, 'fraction_numerical': 0.03555873366703504, 'mean_word_length': 2.40154781647319, 'pattern_counts': {'":': 3, '<': 53, '<?xml version=': 0, '>': 49, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '634138', 'license': 'MIT', 'n_tokens_mistral': 12754, 'n_tokens_neox': 11705, 'n_words': 2858, 'path': 'NLP_with_PyTorch/3_document-embedding/3-1. latent dirichlet allocation.ipynb', 'repo_name': 'naturale0/NLP-Do-It-Yourself', 'size': '30766'}