Sample back the data from trace

Hi, I have a simple regression model with data inputs and sampled coefficients. My question is once I have sampled the regression coefficients, is there a way to sample back the possible data? It’s probably easier to look at the code.

import arviz as az
import jax.numpy as jnp
from jax import ops
import jax.random as random
import numpy as np
import numpyro
import numpyro.distributions as dist
from numpyro.infer import MCMC, NUTS, Predictive, init_to_median
from numpyro import handlers

import matplotlib.pyplot as plt


N = 320
alpha, beta = -2., 0.4
scale = 1.2
rng = np.random.RandomState(32)
a = rng.standard_normal(N)
b = rng.standard_normal(N)
noise = rng.standard_normal(N)*scale
y = alpha*a + beta*b + noise

def model(a=None, b=None, y=None):
    alpha = numpyro.sample('alpha', dist.Normal(0, 0.5))
    beta = numpyro.sample('beta', dist.Normal(0, 0.5))
    scale = numpyro.sample('scale', dist.Exponential(1.))
    loc = alpha*a + beta*b
    obs = numpyro.sample('obs', dist.Normal(loc, scale), obs=y)

rng_sample, rng_post = random.split(random.PRNGKey(32), 2)
hmc = MCMC(
    NUTS(model, init_strategy=init_to_median),
), extra_fields=("z", "energy", "diverging"), y=y, a=a, b=b)
ppc = Predictive(model, hmc.get_samples())(rng_post, a=a, b=b)
idata = az.from_numpyro(hmc,  posterior_predictive=ppc, )
print( az.summary(idata) )

a_new = rng.standard_normal(N//2)
b_new = rng.standard_normal(N//2)
noise_new = 0.
y_new = alpha*a_new + beta*b_new + noise_new
ppc_new = Predictive(model, hmc.get_samples())(rng_post, a=a_new, b=b_new)

do1 =, data={'obs': y_new, 'a':a_new})
p1 = Predictive(do1, hmc.get_samples())(rng_post, a=a_new, b=XXXXX )

As expected ppc_new['obs'] gives values close to y_new.
However is it possible to provide y_new, a_new and get back samples of b_new ?