Using posterior as likelihood

Hi, I’m trying to use a posterior distribution as a likelihood (first in a new model but eventually I want to merge these steps into one model). Is there a way to handle the posterior from model_01 as a likelihood in model_02?

My toy code looks like this:

#simulating some data
N = 1000
stimulus_simulation = np.random.normal(0, 1, N);

sigma_true = 0.6
obs_gstim_simulation = np.random.normal(stimulus_simulation, sigma_true);

sigma_hat_true = 0.8
response_gobs_mean = 1/(1 + sigma_true ** 2) * obs_gstim_simulation
response_gobs_sigma = np.sqrt(sigma_true ** 2 / (1 + sigma_true ** 2) + sigma_hat_true**2)
response_gobs_simulation = np.random.normal(response_gobs_mean, response_gobs_sigma);

#first model to have a posterior: P(x|x_tilde)
with pm.Model() as model_01:
#prior for x
x = pm.Normal(“x”, mu=0, sigma=1)
#prior for sigma
sigma = pm.InverseGamma(“sigma”, alpha=2, beta=1)
x_tilde = pm.Normal(“x_tilde”, mu=x, sigma=sigma, observed = obs_gstim_simulation)
alpha = 1

trace_01 = pm.sample(10000)

#second model to have a posterior for x_hat, using the posterior for x as a likelihood
with pm.Model() as model_02:
#x is observed
x = pm.Normal(“x”, mu=0, sigma=1, observed = stimulus_simulation)
#prior for sigma
sigma = pm.InverseGamma(“sigma”, alpha=2, beta=1)
#prior for x_tilde
x_tilde = pm.Normal(“x_tilde”, mu=x, sigma=sigma)
#this where I don’t know how to use the posterior from model_01
posterior_x = trace_01.posterior.x
#I assume I need a wrapper function for the posterior that would have the below logic
x_hat = pm.UserDefined(“x_hat”, posterior_x, observed = response_gobs_simulation)

trace_02 = pm.sample(10000)

1 Like

I think the current approach is to use pymc.Interpolated — PyMC 4.3.0+0.gea721e4a.dirty documentation

Basically create a random variable using the posterior samples (the smoothed histogram of the marginal posterior distribution), it might not work well especially the posterior is highly correlated

1 Like

There is also a new utility in pymc_experimental that uses a MvNormal approximation IIRC: API Reference — pymc_experimental 0.0.1 documentation