Hello. Not sure where this fits, and this may be an easy question. Suppose I run a hierarchical model with random intercepts. I want to generate a plot of the population distribution, accounting for parameter uncertainty on *mu* and *sigma*, the population parameters. How do I do that?

I can extract all the posterior draws using the **posterior** package in a data frame. Every pair of *mu* and *sigma* posterior draws would give me a population distribution. Is there a way to plot all the population distributions in a graph? Does the **bayesplot** package have a command for this? Or do I need a separate command in R? Any feedback is appreciated.

Probably a density plot similar to the style of plots in `pp_check()`

would be effective. You could have a highly opaque density line for the posterior median, and maybe overlay the density from 100 random samples or something. Perhaps Draw a function as a continuous curve — geom_function • ggplot2 would be a simple way to do it?

I posted the same on Reddit and got a response from a bayerplot dev. I am sharing it here and I will test this out later. But I believe this can achieve what I asked:

Bayesplot dev here 👋, but bayesplot will be no help for you. ggdist is probably your best bet.

`geom_function()`

does not let you treat function arguments as plotting aesthetics, so you can only plot one function’s curve at a time.

That said, we can fake what `geom_function()`

does by computing the densities ourselves beforehand and using `geom_line()`

.

```
library(dplyr)
library(ggplot2)
mus <- rnorm(10, 4, 1)
sigmas <- abs(rnorm(10, .5, .1))
x <- data.frame(
draw = seq_along(mus),
mu = mus,
sigma = sigmas
)
# Find the overall x limits of the plot
x_lim <- range(c(
qnorm(p = .0001, x$mu, x$sigma),
qnorm(p = .9999, x$mu, x$sigma)
))
# Draw the curves' densities at these x values
x_grid <- seq(x_lim[1], x_lim[2], length.out = 1000)
# Compute the densities on the x-grid for each curve
x <- x |>
group_by(draw, mu, sigma) |>
summarise(
x = x_grid,
y = dnorm(x, mu, sigma)
) |>
ungroup()
ggplot(x) +
geom_line(
aes(x = x, y = y, group = draw)
)
```