<jupyter_start><jupyter_text>Based on https://github.com/Savvysherpa/slda .  Modified to the prediction of multivariate normal (diagonal covariance) responses.<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
V = 25   # number of vocabulary
K = 10   # number of topics
N = 100  # number of words in each document
D = 1000 # total number of documents<jupyter_output><empty_output><jupyter_text>Generate topics<jupyter_code>topics = []
topic_base = np.concatenate((np.ones((1, 5)) * 0.2, np.zeros((4, 5))), axis=0).ravel()
for i in range(5):
    topics.append(np.roll(topic_base, i * 5))
topic_base = np.concatenate((np.ones((5, 1)) * 0.2, np.zeros((5, 4))), axis=1).ravel()
for i in range(5):
    topics.append(np.roll(topic_base, i))
topics = np.array(topics)

fig, axes = plt.subplots(figsize=(12,6), nrows=2, ncols=5)
for k in range(K):
    row_ind = int(k / 5)
    col_ind = k % 5
    axes[row_ind, col_ind].imshow(topics[k, :].reshape((5, 5)))
plt.show()<jupyter_output><empty_output><jupyter_text>Generate documents from topics<jupyter_code>alpha = np.ones(K)
thetas = np.random.dirichlet(alpha, size=D)
topic_assignments = np.array([np.random.choice(range(K), size=N, p=theta)
                              for theta in thetas])
word_assignments = np.array([[np.random.choice(range(V), size=1, p=topics[topic_assignments[d, n]])[0]
                              for n in range(N)] for d in range(D)])
doc_term_matrix = np.array([np.histogram(word_assignments[d], bins=V, range=(0, V - 1))[0] for d in range(D)])

fig, ax = plt.subplots(figsize=(4,3))
ax.imshow(doc_term_matrix, aspect='auto')
plt.show()<jupyter_output><empty_output><jupyter_text>Generate responses<jupyter_code>nu2 = 10
sigma2 = 1

eta1 = np.random.normal(scale=nu2, size=K)
eta2 = np.random.normal(scale=nu2, size=K)
fig, axes = plt.subplots(figsize=(9,3), ncols=2)
axes[0].plot(range(K), eta1, color='b')
axes[1].plot(range(K), eta2, color='r')
axes[0].set_title('eta1'); axes[1].set_title('eta2')
plt.show()
y1 = [np.dot(eta1, thetas[i]) for i in range(D)] + np.random.normal(scale=sigma2, size=D)
y2 = [np.dot(eta2, thetas[i]) for i in range(D)] + np.random.normal(scale=sigma2, size=D)
y = np.hstack((y1[:, np.newaxis], y2[:, np.newaxis]))

fig, ax = plt.subplots()
ax.hist(y[:, 0], bins=20, alpha=.5, label='y1')
ax.hist(y[:, 1], bins=20, alpha=.5, label='y2')
ax.legend()
plt.show()
np.savetxt('y.txt', y)
with open('train.dat', 'w') as ofp:
    for i in range(D):
        ofp.write(' '.join([str(ID)+':'+str(Cnt) for ID, Cnt in zip(range(1,V+1), doc_term_matrix[i,:])])+'\n')<jupyter_output><empty_output><jupyter_text>Estimate parameters<jupyter_code>!../src/mvslda -I 300 -K 10 -Y 2 ./train.dat ./y.txt ./model
# likelihood
lik = np.loadtxt('./model.lik')
fig, ax = plt.subplots()
ax.plot(range(len(lik)), lik)
plt.show()
# phi
phi = np.loadtxt('./model.phi')
fig, axes = plt.subplots(figsize=(12,6), nrows=2, ncols=5)
for k in range(K):
    row_ind = int(k / 5)
    col_ind = k % 5
    axes[row_ind, col_ind].imshow(phi[:, k].reshape((5, 5)))
plt.show()
# topic reordering
from scipy.spatial.distance import cdist
topic_reorder = np.argsort([np.argmin(cdist([phi[:, k]], topics)[0]) for k in range(K)])

phi = phi[:, topic_reorder]
fig, axes = plt.subplots(figsize=(12,6), nrows=2, ncols=5)
for k in range(K):
    row_ind = int(k / 5)
    col_ind = k % 5
    axes[row_ind, col_ind].imshow(phi[:, k].reshape((5, 5)))
plt.show()
# eta
eta = np.loadtxt('./model.eta')
eta = eta[topic_reorder, :]

fig, axes = plt.subplots(figsize=(9,3), ncols=2)
axes[0].plot(range(K), eta1, color='b', label='truth')
axes[0].plot(range(K), eta[:, 0], color='g', linestyle=':', label='predicted')
axes[1].plot(range(K), eta2, color='r', label='truth')
axes[1].plot(range(K), eta[:, 1], color='g', linestyle=':', label='predicted')
axes[0].set_title('eta1'); axes[1].set_title('eta2')
axes[0].legend(); axes[1].legend()
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.20297904569553143, 'fraction_numerical': 0.031810148952284775, 'mean_word_length': 5.0488549618320615, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 21, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '581567', 'license': 'MIT', 'n_tokens_mistral': 1720, 'n_tokens_neox': 1633, 'n_words': 335, 'path': 'example/test_mvslda.ipynb', 'repo_name': 'khigashi1987/mvsLDA', 'size': '3961'}