library(ggplot2) # example use # # plot.for.class(0, 10, 8.5, 25) # # this creates a plot with # # (a) the prior, N(0, 10) # (b) 25 data points drawn from N(8.5, 1) # (c) the posterior p(\mu | data) and # (d) the true mean (a vertical line at x=8.5) plot.for.class <- function(prior.mu, prior.var, true.mu, ndata) { x <- rnorm(ndata, mean=true.mu, sd=1) p <- prior.posterior.plot(prior.mu, prior.var, x) p <- p + geom_vline(x=true.mu, color="green") p } # compute the posterior mean and variance of an unknown mean, based on # a gaussian prior and unit-variance normal likelihood. gaussian.posterior <- function(prior.mu, prior.var, data) { # compute the length of the data n <- length(data) # compute the posterior mean post.mu <- (prior.mu/prior.var + sum(data)) / (1/prior.var + n) # compute the posterior variance post.var <- 1 / (1 / prior.var + n) # return the posterior mean and the variance c(mu=post.mu, var=post.var) } # make a nice plot of the prior, posterior, and data prior.posterior.plot <- function(prior.mu, prior.var, data, .n=1000) { # compute the grid at which we evaluate the distributions range.data <- range(data) interval <- range.data[2] - range.data[1] from <- min(prior.mu - 4 * sqrt(prior.var), min(data) - interval * 0.10) to <- max(prior.mu + 4 * sqrt(prior.var), max(data) + interval * 0.10) grid <- seq(from, to, length = .n) # compute the posterior parameters post <- gaussian.posterior(prior.mu, prior.var, data) # make functions for the prior and posterior dpost <- function (x) dnorm(x, mean= post["mu"], sd=sqrt(post["var"])) dprior <- function (x) dnorm(x, mean=prior.mu, sd=sqrt(prior.var)) # make a data frame that applies the prior and posterior to the grid df <- data.frame(x=grid, prior=dprior(grid), post=dpost(grid)) # plot the prior and posterior p <- ggplot() p <- p + geom_line(data=df, aes(x=x, y=prior), colour="red") # prior p <- p + geom_line(data=df, aes(x=x, y=post), colour="blue") # posterior p <- p + geom_rug(data=data.frame(x=data), aes(x=x)) # data p <- p + labs(x="", y="") # p <- p + theme_bw() p }