Visualizing the population distributions for the random intercepts and/or slopes?

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)
)
``````