How to plot another line plot in the same plot of the density curve obtained using RStan

Interfaces
,
1

I am using “stan_dens(fit)” to plot the posterior distribution generated using RStan. I am doing this for the Beta Bernoulli model so I want to compare the results obtained with the closed-form solution. How can I plot the prior, posterior (Using Stan), and the closed-form posterior in a single plot? Currently, I am using the plot function to plot the lines plot for the prior and closed-form posterior in a single plot.

2

rstan_dens returns a ggplot so you should be able to just add additional layers to the plot. For example, if you have a data frame of posterior draws for the priors (with column names that are the same as the parameter names in the original stan_dens plot, so that density curves for the priors are drawn in the correct facet), you should be able to do something like this:

library(tidyverse)

p1 = stan_dens(fit)

# Assuming you have a "wide" data frame with prior draws that includes one 
#  column for each parameter
prior.draws.long = prior.draws
  pivot_longer(everything()) 

p1 +
  geom_density(data=prior.draws, colour="grey50")

Similarly, you can use your analytical solution to generate posterior densities for each parameter and put those in a data frame (something like what you get from posterior::as_draws_df(fit), but with the analytical values), and then add another set of density curves to the plot.

3

Thanks a lot, but I am completely new to both R and RStan. I tried following your solution but could not really implement it. Can you please suggest to me how can I plot f_y1, f_y2, fit and fit2 's distribution plots for theta in a single plot? Here is my code:

n <- 1000
theta<-0.4
y <- rbinom(n,size=1,prob=theta)
hist(y)

yhat<-mean(y)
alpha0 <- 10
beta0 <- 15
alpha_n<-sum(y)+alpha0
beta_n<-n-sum(y)+beta0

x = seq(0, 1, 1 / 1000)
f_y1 <- dbeta(x, shape1 = alpha0, shape2 = beta0)
f_y2 <- dbeta(x, shape1 = alpha_n, shape2 = beta_n)


library(rstan)

model <- stan_model('CoinToss.stan')
fit<-sampling(model,list(n=n,y=y,alpha0=alpha0,beta0=beta0),iter=1500,chains=8,warmup=500)
fit2<-vb(model,data=list(n=n,y=y,alpha0=alpha0,beta0=beta0))

Here is the model code :

data {
  int<lower=0> n;
  int<lower=0, upper=1> y[n];
  real alpha0;
  real beta0;
}
parameters {
  real<lower=0, upper=1> theta;
}
model {
  y ~ bernoulli(theta);
  theta ~ beta(alpha0, beta0);
}
4

Let me know if something like the following is what you had in mind. Note, however, that the code will be slightly different if the model has more than one parameter.

d = data.frame(x, f_y1, f_y2)
# Convert d to long format
d = d %>%
  pivot_longer(-x)

# Start with stan_dens plot and add a layer
stan_dens(fit2, fill=hcl(240,10,95), colour="grey60") +
  # Add f_y1 and f_y2 to stan_dens plot
  geom_line(data=d, aes(x, value, colour=name)) +
  # Adjust x-axis range
  coord_cartesian(xlim=c(0.1, 0.7)) +
  # Change default colors
  scale_colour_manual(name=NULL, values=c("blue", "red"))

Rplot04

# Create the entire plot from scratch
as_draws_df(fit2) %>%
  ggplot(aes(theta)) +
  geom_density(aes(fill="Model"), colour="grey60") +
  # Add f_y1 and f_y2 to plot
  geom_line(data=d, aes(x, value, colour=name)) +
  # Adjust x-axis range
  coord_cartesian(xlim=c(0.1, 0.7)) +
  # Change default colors
  scale_colour_manual(name=NULL, values=c("blue", "red", hcl(240,20,90))) +
  scale_fill_manual(name=NULL, values=hcl(240,20,90)) +
  labs(y=NULL) +
  theme_minimal() +
  theme(axis.line.x=element_line(size=1),
        axis.text.y=element_blank(),
        axis.ticks.x=element_line(),
        axis.ticks.length.x=unit(1,"mm"),
        panel.grid=element_blank())

Rplot05