Hello,

I’m trying to fit a distribution to my data while incorporating prior knowledge. Based on my limited knowledge of Bayesian modeling, this would be equivalent to obtaining the posterior predictive distribution. This is my model:

```
data = np.array([1, 1, 1, 62, 56, 40, 7600, 2, 7, 1, 43, 98, 76])
with pm.Model() as example:
x = pm.Data('x', data, mutable=True)
alpha = pm.TruncatedNormal('alpha', mu=1.02, sigma=1, lower=1.05, upper=1.52)
likelihood = pm.Pareto('likelihood', alpha=alpha, m=x.min(), observed=x)
with example:
trace = pm.sample(idata_kwargs=dict(log_likelihood=True),)
pm.sample_posterior_predictive(trace, extend_inferencedata=True)
```

I plotted the mean posterior predictive using az.plot_ppc(), but the plot looks jagged regardless of the num_pp_samples used (2nd image below). Looking at the Axis object returned by plot_ppc(), only 512 points are plotted while the long tails of a Pareto distribution cover a huge range of points.

While looking for alternatives, I thought, isn’t everything in my model feeding into my likelihood variable? Wouldn’t the likelihood variable post-sampling be equivalent to the posterior predictive? I plotted pm.logp(likelihood).exp(), but it looks significantly different from az.plot_ppc() (plot below). Perhaps there’s a conceptual difference between the two methods that I’m not picking up on.

Here’s how I plotted the two distributions for reference:

```
# Plot density via likelihood variable's logp function
domain = np.linspace(0, data.max()+10, 1_000_000)
probs = pm.logp(likelihood, domain).exp().eval()
plt.plot(domain, probs, label='pm.logp()')
# Plot density via Arviz's plot_ppc
ppc = az.plot_ppc(trace, observed=False, num_pp_samples=4_000, show=False); plt.close()
x0, dens0 = ppc.lines[-1].get_xdata(), ppc.lines[-1].get_ydata() # lines[-1] corresponds to mean posterior predictive
plt.plot(x0, dens0, label='az.plot_ppc()')
# Plot observations for reference
plt.scatter(data, data*0, color='black', label='Observations'); plt.xlabel('x'); plt.ylabel('Density'); plt.legend()
```

To summarize, there are two points I’m not clear on:

- Why is the plot given by az.plot_ppc() different from the one by pm.logp()?
- Is there a way to make the Arviz plot smoother?

Appreciate any help with this!