@Simon, I really appreciate your effort in helping me. However, when I created the plot to check the prior distribution, I still observed a flat prior. I made some modifications to your code so that I could visualize both the prior and posterior predictive distributions:
with pm.Model() as mod:
alpha = pm.HalfNormal("alpha", 1) #intercept
beta = pm.HalfNormal("beta", 1) #slope
mu = pm.Deterministic("mu", alpha + beta*x) #location
sigma = pm.HalfNormal("sigma", 1) #sd
y_hat = pm.LogNormal("y", mu=mu, sigma=sigma, observed=y)
trace = pm.sample(idata_kwargs={"log_likelihood": True})
prior_pred = pm.sample_prior_predictive(samples=10)
post_pred = pm.sample_posterior_predictive(trace,extend_inferencedata=True)
As you can see, the prior predictive mean remains a horizontal line. My primary concern is as follows: In cases where I lack the time or resources to conduct new laboratory tests, or when the available data is limited in size, will my prior information be sufficiently robust to yield reasonable parameter estimates? Consequently, in my real-world scenario, the prior distribution plays a pivotal role.
I made slight adjustments to your example to make the data more closely resemble my actual data:
import pymc as pm
import arviz as az
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
#really simplistic simulation of data distributed as two correlated LogNormals
r = 0.8
alpha = 4
sigma = 0.8
x0 = np.sort(pm.draw(pm.LogNormal.dist(alpha, sigma), draws=100, random_seed=27))
x = np.sort(pm.draw(pm.LogNormal.dist(alpha, sigma), draws=100, random_seed=77))
y = r*x + np.sqrt(1-r**2)*x0
# Plot the data
_, ax = plt.subplots(1,2, figsize=(8, 4))
ax[0].plot(x, y, 'C0.')
ax[0].set_xlabel('x')
ax[0].set_ylabel('y', rotation=0)
az.plot_kde(y, ax=ax[1])
ax[1].set_xlabel('y')
ax[1].set_ylabel('p(y)')
plt.title('Real Data')
plt.tight_layout()
# Create the model
with pm.Model() as mod:
alpha_m = pm.HalfNormal("alpha", alpha) #intercept
beta_m = pm.HalfNormal("beta", alpha) #slope
mu = pm.Deterministic("mu", alpha_m + beta_m*x) #location
sigma = pm.HalfNormal("sigma", 170) #sd
y_hat = pm.LogNormal("y", mu=mu, sigma=sigma, observed=y)
trace = pm.sample(idata_kwargs={"log_likelihood": True})
prior_pred = pm.sample_prior_predictive(samples=10)
post_pred = pm.sample_posterior_predictive(trace,extend_inferencedata=True)
# Plot results
y_pred_post = post_pred.posterior_predictive['y'].mean(dim=["chain", "draw"])
_, ax = plt.subplots(4, 2, figsize=(10, 12))
plt.subplots_adjust(hspace=0.5)
# row 0
az.plot_bpv(post_pred,kind="p_value",ax=ax[0,0])
az.plot_bpv(post_pred, kind="u_value", ax=ax[0,1])
# row 1
ax[1,0].set_title("mean")
az.plot_bpv(post_pred, kind="t_stat", t_stat="mean",ax=ax[1,0])
ax[1,1].set_title("standard deviation")
az.plot_bpv(post_pred, kind="t_stat", t_stat="std", ax=ax[1,1])
# row 2
ax[2,0].set_title("Prior Predictions Distribution")
az.plot_ppc(prior_pred, ax=ax[2,0],alpha=0.01,colors=["gray","blue",'red'], mean=True, legend=True,group='prior')
# ax[2,0].set_xlim([-1,1000])
ax[2,1].set_title("Observed Data Distribution")
ax[2,1].hist(y, bins=20,color='gray',alpha=0.7, density=True)
sns.kdeplot(y, ax=ax[2, 1], color='blue', linestyle='--')
# row 3
ax[3,0].set_title("Posterior Predictions Distribution")
az.plot_ppc(post_pred, ax=ax[3,0],alpha=0.01,colors=["gray","blue",'red'], mean=True, legend=True,group='posterior')
ax[3,1].set_title("Observed x Predicted Data Distribution")
ax[3,1].scatter(x,y,alpha=0.4,color='blue',label='Observed')
ax[3,1].scatter(x,y_pred_post,alpha=0.4,color='gray',label='Predicted')
ax[3,1].legend()
# plot trace
az.plot_trace(trace)
plt.tight_layout()
# plot parameters
def plot_pretty_posteriors(trace):
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
# Plot for 'Intercept'
az.plot_posterior(trace, var_names='alpha', color="#7A68A6",kind='hist',hdi_prob=0.95,bins=25, ax=axes[0], textsize=12)
axes[0].set_title('Intercept', fontsize=16)
# Plot for 'Slope'
az.plot_posterior(trace, var_names='beta', color="#A60628",kind='hist',hdi_prob=0.95,bins=25, ax=axes[1], textsize=12)
axes[1].set_title('Slope', fontsize=16)
# Plot for 'sigma'
az.plot_posterior(trace, var_names='sigma', color="#348ABD",kind='hist',hdi_prob=0.95,bins=25, ax=axes[2], textsize=12)
axes[2].set_title('sigma', fontsize=16)
plt.tight_layout()
plt.show()
plot_pretty_posteriors(trace)
Then, I compared it with my previous model posted at the beginning of this discussion using a leave-one-out validation:
and it seems the first model is performing better.
I just need to have the prior predictive mean following a halfnormal and not a horizontal line, but I didn’t figure out how yet 
