<jupyter_start><jupyter_text>Generate synthetic data and compare Gibbs to RJMCMC<jupyter_code>true_model = NeymanScottModel(mu, alpha, beta, lmbda0, 
                              bkgd_class(**bkgd_hypers), 
                              DefaultClusterFactory(cluster_class, cluster_hypers))

# Sample from the model
sample_kwargs = dict(time_kwargs=dict(), mark_kwargs=dict())
data, parents, clusters, weights = true_model.generate(**sample_kwargs)
nodes, times, marks = data[:, 0], data[:, 1], data[:, 2:]

# Only keep the valid times
# valid = times < T
valid = np.ones_like(times, dtype=bool)
nodes = nodes[valid]
times = times[valid]
marks = marks[valid]
parents = parents[valid]

# Sort the data by time
perm = np.argsort(times)
nodes = nodes[perm]
times = times[perm]
marks = marks[perm]
parents = parents[perm]
data = np.column_stack((nodes, times, marks))
%debug
plt.figure(figsize=(12, 15))
ylim = 1.1 * abs(marks[:, 0]).max()
for m in range(M):
    plt.subplot(M, 1, m+1)
    for z in np.unique(parents):
        inds = np.where((parents == z) & (nodes == m))[0]
        plt.plot(times[inds], marks[inds, 0], 'o' if z >= 0 else 'ko')
    
    plt.ylabel("Mark Dim 1")
    plt.ylim(-ylim, ylim)
    if m == M-1:
        plt.xlabel("time")
    else:
        plt.xticks([])
        
    plt.title("Node {}".format(m+1))

plt.tight_layout()
# # Fit the correct model
# num_samples = 100
# # nsm = NeymanScottModel(mu, alpha, beta, lmbda0, 
# #                        bkgd_class(**bkgd_hypers), 
# #                        DefaultClusterFactory(cluster_class, cluster_hypers))

# # Fit the model
# print("true num clusters: ", len(clusters))

# gibbs_samples = true_model.gibbs_sample_posterior(data, num_samples=num_samples, init_method="prior")
# Fit the correct model with RJMCMC
num_samples = 10000
# nsm = NeymanScottModel(mu, alpha, beta, lmbda0, 
#                        bkgd_class(**bkgd_hypers), 
#                        DefaultClusterFactory(cluster_class, cluster_hypers))
rjmcmc_samples = true_model.rjmcmc_sample_posterior(data, num_samples=num_samples, init_method="background")
# Plot the number of clusters over samples
plt.figure()
# plt.plot([s["num_clusters"] for s in gibbs_samples], label="nsm")
plt.plot([s["num_clusters"] for s in rjmcmc_samples], label="nsm")
plt.plot([0, len(rjmcmc_samples)], [len(clusters), len(clusters)], '--k', label="True")
plt.xlabel("Iteration")
plt.ylabel("Num clusters")
plt.legend(loc="lower right")
plt.show()
k = int(mu)
k_smpls = [k]
n_smpls = 20000
for i in range(n_smpls):
    move = npr.rand() < 0.5
    if move == 1:
        p_acc = mu / (k + 1)
        if npr.rand() < p_acc:
            k += 1
            
    else:
        # remove
        p_acc = k / mu
        if npr.rand() < p_acc:
            k -= 1
            
    k_smpls.append(k)
    
plt.figure()
plt.plot(k_smpls)
plt.plot(mu * np.ones(n_smpls), '--')

from scipy.stats import poisson
plt.figure()
plt.hist(k_smpls, density=True)
plt.plot(np.arange(2 * mu), poisson.pmf(np.arange(2 * mu), mu))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1481844946025515, 'fraction_numerical': 0.014066077854105331, 'mean_word_length': 3.206327372764787, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '175548', 'license': 'MIT', 'n_tokens_mistral': 1175, 'n_tokens_neox': 1096, 'n_words': 249, 'path': 'analyses/synthetic/2019_05_14_synthetic_comparison_3.ipynb', 'repo_name': 'slinderman/neymanscott', 'size': '3057'}