Hello Stanizens,

I’m looking for some help this time with regard to some plot engineering.

I’m trying to make a plot consisting of panels of prior and posterior density pairs, with vertical lines of the true \theta values (the data generating process is known in the case of this inference).

Following various introductory guides, I’ve been able to make a plot of just the posterior densities extracted from a CmdStanR fit object. I run code along the lines of below,

```
library(cmdstanr)
library(bayesplot)
fit <- readRDS(...)
posterior <- fit$draws()
post_plots <- mcmc_dens(fit$draws(c("u_M", "a_SD", "a_DS", "a_M", "a_MSC", "k_S_ref", "k_D_ref", "k_M_ref", "Ea_S", "Ea_D", "Ea_M")))
```

for the following quick and unpolished plot:

After this step, I’m a little at a loss for how to

- plot the priors on top of the posteriors and
- plot vertical lines representing the true theta value

across all of the parameter panels in straightforward fashion.

I do have the prior predictive samples in the `generated quantities`

section of the code:

```
generated quantities {
...
// Obtain prior predictive samples.
real u_M_prior_pred = normal_lb_ub_rng(u_M_prior_dist_params[1], u_M_prior_dist_params[1] * prior_scale_factor, u_M_prior_dist_params[2], u_M_prior_dist_params[3]);
real a_SD_prior_pred = normal_lb_ub_rng(a_SD_prior_dist_params[1], a_SD_prior_dist_params[1] * prior_scale_factor, a_SD_prior_dist_params[2], a_SD_prior_dist_params[3]);
real a_DS_prior_pred = normal_lb_ub_rng(a_DS_prior_dist_params[1], a_DS_prior_dist_params[1] * prior_scale_factor, a_DS_prior_dist_params[2], a_DS_prior_dist_params[3]);
real a_M_prior_pred = normal_lb_ub_rng(a_M_prior_dist_params[1], a_M_prior_dist_params[1] * prior_scale_factor, a_M_prior_dist_params[2], a_M_prior_dist_params[3]);
real a_MSC_prior_pred = normal_lb_ub_rng(a_MSC_prior_dist_params[1], a_MSC_prior_dist_params[1] * prior_scale_factor, a_MSC_prior_dist_params[2], a_MSC_prior_dist_params[3]);
real k_S_ref_prior_pred = normal_lb_ub_rng(k_S_ref_prior_dist_params[1], k_S_ref_prior_dist_params[1] * prior_scale_factor, k_S_ref_prior_dist_params[2], k_S_ref_prior_dist_params[3]);
real k_D_ref_prior_pred = normal_lb_ub_rng(k_D_ref_prior_dist_params[1], k_D_ref_prior_dist_params[1] * prior_scale_factor, k_D_ref_prior_dist_params[2], k_D_ref_prior_dist_params[3]);
real k_M_ref_prior_pred = normal_lb_ub_rng(k_M_ref_prior_dist_params[1], k_M_ref_prior_dist_params[1] * prior_scale_factor, k_M_ref_prior_dist_params[2], k_M_ref_prior_dist_params[3]);
real Ea_S_prior_pred = normal_lb_ub_rng(Ea_S_prior_dist_params[1], Ea_S_prior_dist_params[1] * prior_scale_factor, Ea_S_prior_dist_params[2], Ea_S_prior_dist_params[3]);
real Ea_D_prior_pred = normal_lb_ub_rng(Ea_D_prior_dist_params[1], Ea_D_prior_dist_params[1] * prior_scale_factor, Ea_D_prior_dist_params[2], Ea_D_prior_dist_params[3]);
real Ea_M_prior_pred = normal_lb_ub_rng(Ea_M_prior_dist_params[1], Ea_M_prior_dist_params[1] * prior_scale_factor, Ea_M_prior_dist_params[2], Ea_M_prior_dist_params[3]);
}
```

Ideally, I’d like to do things in vectorised fashion in something along the lines of `post_plot + priors_plot + vertical_lines_plot`

, but perhaps I need to do this iteratively along the lines of the following pseudocode,

```
total_param_count <- 11
plot_list <- vector(mode = "list", length = total_param_count)
param_names_list = c("u_M", ...)
names(plot_list) <- param_names_list
for (i in 1:total_param_count) {
plot_list[i] <- custom_function_for_single_prior_posterior_plot(prior_draws[param_names_list[i]], posterior_draws[param_names_list[i]], ...)
}
```

?

Has anyone had experience doing a plot like this, and if so, would you happen to have example code you can share or point me in the right direction? Ultimately, I’m trying to get to a plot that resembles the following:

On a less-involved note, how does one tune the alpha of the fill in `mcmc_dens`

? `mcmc_dens`

by default does not come with an alpha argument.

Thank you all kindly for your patient help.