Model scoring on out-of-samples predictions

I am trying to score a pymc model on different datasets. Here is how I am generating the trace -

with model_factory(df):
    trace = pm.sampling.jax.sample_numpyro_nuts(

I need to sum up all the values of obs variable after taking mean. Here is how I am feeding the test dataset. Is this the right way to get predictions on test datasets

with model_factory(df):
    pm.set_data({"input_data": test_df})
    trace_new = pm.sample_posterior_predictive(trace)
trace_new.posterior_predictive['obs'].mean(dim=["chain", "draw"]).sum().values

Yes, this looks correct. But it’s a bit of a waste to reduce your posteriors to a point estimate. If you’re going to do that, why bother running the expensive and time consuming MCMC algorithm in the first place?

You should consider computing a model score for each posterior draw, and look at the posterior distributions over the scoring functions.