Hey! First off: thanks to all the contributors to Pyro for allowing laymen like me to engage in the world of probabilistic programming. I have a question regarding making predictions.

- What tutorial are you running?

The Bayesian Regression tutorial. - What version of Pyro are you using?

0.3.0 - Please link or paste relevant code, and steps to reproduce.

When doing predictions, I’ve seen multiple approaches. In the tutorials, TracePredictive is used. However, I’ve also seen Jupyter Notebooks and blogposts in which the authors just sample the Guide n times and use those to predict. Which one is better? When applying both to my own dataset My r^2 increases with 0.15 when using the latter approach.

**Code using the predictive distribution:**

```
get_marginal = lambda traces, sites:EmpiricalMarginal(traces, sites)._get_samples_and_weights()[0].detach().cpu().numpy()
def summary(traces, sites):
marginal = get_marginal(traces, sites)
site_stats = {}
for i in range(marginal.shape[1]):
site_name = sites[i]
marginal_site = pd.DataFrame(marginal[:, i]).transpose()
describe = partial(pd.Series.describe, percentiles=[.05, 0.25, 0.5, 0.75, 0.95])
site_stats[site_name] = marginal_site.apply(describe, axis=1) \
[["mean", "std", "5%", "25%", "50%", "75%", "95%"]]
return site_stats
def wrapped_model(x_data, y_data):
pyro.sample("prediction", Delta(model(*test_data.tensors)))
posterior = svi.run(*train_data.tensors)
# posterior predictive distribution we can get samples from
trace_pred = TracePredictive(wrapped_model,
posterior,
num_samples=1000)
post_pred = trace_pred.run(x_test.tensors[0], None)
post_summary = summary(post_pred, sites= ['prediction'])
y = post_summary["prediction"]
predictions = pd.DataFrame({
"y_mean": y["mean"],
"y_std": y["std"],
"y_perc_5": y["5%"],
"y_perc_95": y["95%"],
"true_los": y_data.squeeze(),
})
r2_score(predictions['true_los'], predictions['y_mean'])
```

Out: 0.306

**Sampling the guide:**

```
preds = []
for i in range(1000):
sampled_reg_model = guide(None, None)
pred = sampled_reg_model(X_test_tensor).data.numpy().flatten()
preds.append(pred)
all_preds = np.stack(preds).mean(0)
r2_score(y_test, all_preds)
```

Out: 0.472

Help would be greatly appreciated!

Thanks.