<jupyter_start><jupyter_text>Setup and EDA  Generate Data Sets<jupyter_code>independent_data = generate_data()
correlated_data = generate_data(0.9)<jupyter_output><empty_output><jupyter_text>Plot Mean and Variance Over Time<jupyter_code>fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 3))

ax1.plot(independent_data.days, independent_data.mu, label="mean")
ax1.set_xlabel("Date")
ax1.set_ylabel("Log Event Count")
ax1.set_title("E(Z)")

ax2.plot(independent_data.days, independent_data.sigma, label="standard deviation")
ax2.set_xlabel("Date")
plt.ylabel("Log Event Count")
ax2.set_title("sd(Z)")

plt.tight_layout()
plt.savefig("./figures/mean_var_functions.png")
plt.show()
# The data set is underdispersed
print(np.exp(5))
print(np.sqrt(np.exp(5)))
print(np.exp(0.15))<jupyter_output>148.4131591025766
12.182493960703473
1.161834242728283<jupyter_text>Compare Data Sets<jupyter_code>fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 3))

ax1.plot(independent_data.days, independent_data.y)
ax1.set_ylim((0, 1.1*np.max(independent_data.y)))
ax1.set_xlabel("Date")
ax1.set_ylabel("Event Count")
ax1.set_title("Independent Event Counts")

ax2.plot(correlated_data.days, correlated_data.y)
ax2.set_ylim((0, 1.1*np.max(correlated_data.y)))
ax2.set_xlabel("Date")
ax2.set_ylabel("Event Count")
ax2.set_title("Correlated Event Counts")

plt.tight_layout()
fig.savefig("./figures/events_over_time.png")
plt.show()<jupyter_output><empty_output><jupyter_text>Model-fitting Independent Data<jupyter_code>gibbs_ind_ind = GibbsSampler()
gibbs_ind_cor = GibbsSampler(alpha_rho=1, beta_rho=1)
n_iter = 2000
gibbs_ind_ind.fit(independent_data, n_iter=n_iter)
pickle.dump(gibbs_ind_ind, open("gibbs_ind_ind.pkl", "wb")) 
# gibbs_ind_ind = pickle.load(open("gibbs_ind_ind.pkl", "rb"))
gibbs_ind_cor.fit(independent_data, n_iter=n_iter)
pickle.dump(gibbs_ind_cor, open("gibbs_ind_cor.pkl", "wb")) 
# gibbs_ind_cor = pickle.load(open("gibbs_ind_cor.pkl", "rb"))<jupyter_output><empty_output><jupyter_text>Correlated Data<jupyter_code>gibbs_cor_ind = GibbsSampler()
gibbs_cor_cor = GibbsSampler(alpha_rho=1, beta_rho=1)
gibbs_cor_ind.fit(correlated_data, n_iter=n_iter)
pickle.dump(gibbs_cor_ind, open("gibbs_cor_ind.pkl", "wb")) 
# gibbs_cor_ind = pickle.load(open("gibbs_cor_ind.pkl", "rb"))
gibbs_cor_cor.fit(correlated_data, n_iter=n_iter)
pickle.dump(gibbs_cor_cor, open("gibbs_cor_cor.pkl", "wb"))
# gibbs_cor_cor = pickle.load(open("gibbs_cor_cor.pkl", "rb"))<jupyter_output><empty_output><jupyter_text>Analysis<jupyter_code>burnin = int(n_iter/4)
n_days = independent_data.days.size
# Function for plotting E(Z) and sd(Z) for two models
def plot_mean_var(mod1, mod2, filename, proportion=1.):
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 3))

    first_day_idx = int((1-proportion)*n_days)
    ax1.plot(
        independent_data.days[first_day_idx:],
        independent_data.mu[first_day_idx:],
        label="Truth"
    )
    ax1.plot(
        independent_data.days[first_day_idx:],
        mod1.mu.values[burnin:, first_day_idx:].mean(axis=0),
        label="Independent Model",
#         linestyle="dashed",
    )
    ax1.plot(
        independent_data.days[first_day_idx:],
        mod2.mu.values[burnin:, first_day_idx:].mean(axis=0),
        label="Correlated Model",
#         linestyle="dashdot",
    )

    ax1.set_xlabel("Date")
    ax1.set_ylabel("Log Event Count")
    ax1.set_title("E(Z)")
    ax1.legend()

    ax2.plot(
        independent_data.days[first_day_idx:],
        independent_data.sigma[first_day_idx:],
        label="Truth"
    )
    ax2.plot(
        independent_data.days[first_day_idx:],
        mod1.sigma.values[burnin:, first_day_idx:].mean(axis=0),
        label="Independent Model",
#         linestyle="dashed",
    )
    ax2.plot(
        independent_data.days[first_day_idx:],
        mod2.sigma.values[burnin:, first_day_idx:].mean(axis=0),
        label="Correlated Model",
#         linestyle="dashdot",
    )
    ax2.set_xlabel("Date")
    plt.ylabel("Log Event Count")
    ax2.set_title("sd(Z)")
    ax2.legend()
    
    plt.tight_layout()
    plt.savefig(filename)
    plt.show()
plot_mean_var(gibbs_ind_ind, gibbs_ind_cor, "./figures/ind_data_mean_var.png", proportion=1)
plot_mean_var(gibbs_cor_ind, gibbs_cor_cor, "./figures/cor_data_mean_var.png")
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 3))

ax1.plot(independent_data.days, independent_data.y, label="Data")
ax1.plot(
    independent_data.days,
    np.exp(np.median(gibbs_ind_ind.mu.values[burnin:, :], axis=0)),
    label="Independent Model Fit",
)
ax1.plot(
    independent_data.days,
    np.exp(np.median(gibbs_ind_cor.mu.values[burnin:, :], axis=0)),
    label="Correlated Model Fit",
)

ax1.set_ylim((0, 1.1*np.max(independent_data.y)))
ax1.set_xlabel("Date")
ax1.set_ylabel("Event Count")
ax1.set_title("Independent Event Counts")
ax1.legend()

ax2.plot(correlated_data.days, correlated_data.y, label="Data")
ax2.plot(
    independent_data.days,
    np.exp(np.median(gibbs_cor_ind.mu.values[burnin:, :], axis=0)),
    label="Independent Model Fit",
)
ax2.plot(
    independent_data.days,
    np.exp(np.median(gibbs_cor_cor.mu.values[burnin:, :], axis=0)),
    label="Correlated Model Fit",
)

ax2.set_ylim((0, 1.1*np.max(correlated_data.y)))
ax2.set_xlabel("Date")
ax2.set_ylabel("Event Count")
ax2.set_title("Correlated Event Counts")
ax2.legend()

plt.tight_layout()
fig.savefig("./figures/model_fits.png")
plt.show()
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 3))
ax1.plot(gibbs_ind_cor.rho.values[burnin:])
ax1.set_xlabel("Index")
ax1.set_ylabel("ρ")
ax1.set_title("Independent Data")

ax2.plot(gibbs_cor_cor.rho.values[burnin:])
ax2.set_xlabel("Index")
ax2.set_ylabel("ρ")
ax2.set_title("Correlated Data")

plt.tight_layout()
fig.savefig("./figures/rho_trace_plots.png")
plt.show()
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 3))
i = 100
ax1.plot(gibbs_ind_ind.z.values[1:, i])
ax1.set_xlabel("Index")
ax1.set_ylabel(f"Z_{i+1}")
ax1.set_title("Trace Plot")

y_i = independent_data.y[i]
mu_i = gibbs_ind_ind.mu.values[burnin:, i].mean()
sigma_i = gibbs_ind_ind.sigma.values[burnin:, i].mean()
z_i_vals = np.linspace(mu_i - 2*sigma_i, mu_i + 2*sigma_i)
z_i_distribution = norm(loc=mu_i, scale=sigma_i)
density_scaling_i = z_i_distribution.cdf(np.log(y_i+1)) - z_i_distribution.cdf(np.log(y_i))
ax2.plot(z_i_vals, norm.pdf(z_i_vals, loc=mu_i, scale=sigma_i)/density_scaling_i)
ax2.set_xlabel(f"Z_{i+1}")
ax2.set_ylabel("Density")
ax2.hist(gibbs_ind_ind.z.values[burnin:, i], density=True, bins=4)
ax2.set_title("Estimated Density Plot with Sampled Values")
plt.tight_layout()
fig.savefig("./figures/z_plots.png")

plt.show()
def trace_coefs(model, filename):
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 3))
    p = 1
    ax1.plot(model.beta.values[1:, p])
    ax1.set_xlabel("Index")
    ax1.set_ylabel(f"β_{p+1}")
    ax1.set_title("Trace Plot")

    ax2.plot(model.alpha.values[1:, p])
    ax2.set_xlabel("Index")
    ax2.set_ylabel(f"α_{p+1}")
    ax2.set_title("Trace Plot")

    plt.tight_layout()
    fig.savefig(filename)
    plt.show()
trace_coefs(gibbs_ind_ind, "./figures/ind_alpha_beta_plots.png")
trace_coefs(gibbs_cor_cor, "./figures/cor_alpha_beta_plots.png")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1647578308265489, 'fraction_numerical': 0.031599282461708296, 'mean_word_length': 5.583106267029973, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 24, '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': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '199412', 'license': 'MIT', 'n_tokens_mistral': 3206, 'n_tokens_neox': 3044, 'n_words': 376, 'path': 'analysis.ipynb', 'repo_name': 'eastonhuch/discrete-lognormal-sampler', 'size': '7247'}