Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
R code for creating plot showing the effect of different priors
library(tidyverse)
library(viridis)
set.seed(92414)
### Draw samples from prior distributions --------------------------------------
# Non-Information: N(-2, 10)
# Weakly Informative: N(-2, 2)
# Informative: N(-2, 1)
# Prior distributions (for plotting)
plot_data <- data_frame(x = seq(-10, 10, by = 0.01)) %>%
mutate(
`Non-Informative` = dnorm(x, mean = -2, sd = 10),
`Weakly Informative` = dnorm(x, mean = -2, sd = 2),
Informative = dnorm(x, mean = -2, sd = 1)
) %>%
gather(key = "prior_type", value = "prior", -x) %>%
mutate(prior_type = factor(prior_type, levels = c("Non-Informative",
"Weakly Informative",
"Informative")))
# Prior samples
ni_prior <- data_frame(prior_type = "Non-Informative",
prior = rnorm(n = 10000, mean = -2, sd = 10))
wi_prior <- data_frame(prior_type = "Weakly Informative",
prior = rnorm(n = 10000, mean = -2, sd = 2))
in_prior <- data_frame(prior_type = "Informative",
prior = rnorm(n = 10000, mean = -2, sd = 1))
prior_samples = bind_rows(ni_prior, wi_prior, in_prior) %>%
mutate(prior_type = factor(prior_type, levels = c("Non-Informative",
"Weakly Informative",
"Informative")))
### Calculate likelihood -------------------------------------------------------
# Likelihood for sample data
sample_data <- rnorm(n = 1000, mean = 5, sd = 2.5)
plot_data <- plot_data %>%
mutate(likelihood = dnorm(x = mean(sample_data), mean = x, sd = 2.5))
# Likelihood for prior draws
prior_samples <- mutate(prior_samples,
prior_likelihood = dnorm(x = mean(sample_data),
mean = prior, sd = 2.5))
### Calculate posteriors -------------------------------------------------------
posterior_samples <- prior_samples %>%
group_by(prior_type) %>%
sample_n(size = 10000, replace = TRUE, weight = prior_likelihood)
### Create plot ----------------------------------------------------------------
ggplot(plot_data, aes(x = x)) +
facet_wrap(~ prior_type, nrow = 1) +
geom_line(aes(y = prior, color = "Prior"), linetype = "dashed") +
geom_line(aes(y = likelihood, color = "Likelihood"), linetype = "dashed") +
geom_density(data = posterior_samples,
aes(x = prior, color = "Posterior", fill = "Posterior"),
alpha = 0.5, size = 1, show.legend = FALSE) +
scale_color_manual(values = c("Prior" = viridis(3)[2],
"Likelihood" = viridis(3)[1],
"Posterior" = viridis(3)[3]),
breaks = c("Prior", "Likelihood", "Posterior"),
guide = guide_legend(override.aes = list(
linetype = "solid", alpha = 1, size = 5
))) +
scale_fill_manual(values = c("Posterior" = viridis(3)[3])) +
coord_cartesian(xlim = c(-10, 10)) +
labs(x = NULL, y = "Density", color = NULL) +
theme_bw() +
theme(legend.position = "bottom")
@wjakethompson

This comment has been minimized.

Copy link
Owner Author

commented Oct 11, 2018

priors_plot

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.