Exercise: Program a simple grid approximation to obtain a posterior distribution

Please download a copy of BayesGrid.Rmd and the data file arrival.csv into the same folder and fill in the relevant sections below.

To download the files follow the link and right-click anywhere on the page and the click “save as”.

The data Data

Use this data about arrivals times relative to the start of a lecture.

my_data = 
  read.table(
    "arrival.csv",
    header = T,
    sep = ",")
head(my_data)
##   age time
## 1  31   -1
## 2  19    0
## 3  24   -1
## 4  23   -1
## 5  24   -5
## 6  21   -1

Here is the plot of the data

with(my_data,
     plot(age, time))

age = my_data$age
hist(age)

The Model

For now, we are only interested in estimating the mean age (age) using the following model. The table shows the incomplete model expressed with equations, quap code fragments, and how one would calculate quantities directly in R.


What Notation quap code R code
Likelihood \(a_i \sim Normal(\mu,\sigma)\) age ~ dnorm(mu, sigma) L[i] = dnorm(age[i], mu, sigma)
Prior \(\mu \sim Normal(25,2)\) mu ~ dnorm(25 , 2) p_mu = dnorm(25 , 2)
Prior \(\sigma \sim Uniform(0,2)\) sigma ~ dunif(0,2) p_sigma ~ dunif(0,2


Your first task is to complete the model to do this. Please add the missing information in the table above. You need to add the parameter values to the prior for \(\mu\) and you have to choose a distribution and parameter(s) for the standard deviation \(\sigma\).

The posterior disribution

The posterior distribution is calculated as follows:

\[ \overset{Posterior}{P(\mu,\sigma|w)} = \frac{\prod_i \overset{Likelihood}{Normal(a_i|\mu,\sigma)} \cdot \overset{Prior}{Normal(\mu|25,2) Uniform(\sigma|0,2)}} {\overset{Evidence}{\int\int \prod_i Normal(a_i|\mu,\sigma) \cdot Normal(\mu|25,2)Uniform(\sigma|0,2)d\mu d\sigma}} \]


Please add the missing information, using what you just added to the table, wherever are three dots \(...\) or it reads \(YourDist\).

Calculate the likelihood for a data point

Let’s assume these parameter values

  • mu = 22
  • sd = 3

Please calculate the likelihood for the first data point in the vector age in the next code block. If you are uncertain about how to proceed, check the table above

# Likelihood for the first data point
L = dnorm(age[1], 22, 3)
L
## [1] 0.001477283

Now calculate the probability of the data point given the parameter values for each recorded age You do not need a loop to do this!

# Likelihoods for all individual data points
L2 = dnorm(age, 22, 3)
L2
##  [1] 0.0014772828 0.0806569082 0.1064826685 0.1257944092 0.1064826685
##  [6] 0.1257944092 0.1064826685 0.0546700249 0.0014772828 0.0806569082
## [11] 0.1257944092 0.1064826685 0.0331590463 0.0806569082 0.1257944092
## [16] 0.0546700249 0.1257944092 0.1064826685 0.0331590463 0.0806569082
## [21] 0.0806569082 0.0005140930 0.1064826685 0.0179969888 0.0546700249
## [26] 0.0087406297 0.0005140930 0.1257944092 0.0546700249 0.0001600902
## [31] 0.0806569082 0.0331590463 0.1257944092 0.0806569082 0.0546700249
## [36] 0.0806569082 0.0331590463 0.1329807601 0.0331590463 0.1257944092
## [41] 0.0331590463 0.0806569082 0.0806569082 0.0806569082 0.0546700249
## [46] 0.1064826685 0.0005140930 0.0037986620 0.1329807601 0.0179969888

And now calculate the joint probability of all recorded ages. (Read the help for the prod function, i.e. type ?prodin the console)

# Joint likelihood for all individual data points
prod(L2)
## [1] 2.930732e-72
joint_likel<- prod(L2)

Now calculate the prior probability of the the parameters mu = 22 and sd = 3. You should use the prior distributions you specified above for this.

# Prior probability for mu = 22
p_mu = dnorm(22,3)
p_mu
## [1] 1.624636e-79
# Prior probability for sd = 3
p_sigma = dexp(3)
p_sigma
## [1] 0.04978707

And the joint prior probability for mu = 22 and sd = 2.

# Joint prior probability for mu = 22 and sd = 3
joint_prob = p_mu*p_sigma

By putting the code you wrote until now, you should be able to calculate the numerator of the equation above, i.e. likelihood times prior probability (density), for parameters mu = 22 and sd = 3.

# Joint probability of data and parameters for  mu = 22 and sd = 3 
# given the prior distributions for mu and sd

prod(joint_prob, joint_likel)
## [1] 2.370548e-152

Grid approximation

Now that you have calculated the joint likelihood of data and parameters given the prior distributions for mu and sd for one set of parameters, you can do this for a larger number of parameters combinations, which we can generate by constructing a grid where each point is a combination of one mu and one sd parameter value. Continue either with Grid approximation, version 1 or Grid approximation, version 2.

Grid approximation, version 1

First set up the grid. You can use the the code example below. Feel free to change the range of the parameters you explore by changing the first and second number in the seq commands

# Here we are making our grid
# expand.grid produces all combinations 
# of the N (here 2) variables submitted
# and punts them in a matrix with N columns
params = 
  expand.grid(
    mu = seq(20,30,.1),
    sd = seq(2.5,4.5, .025)
  )
params = data.frame(params)
head(params)
##     mu  sd
## 1 20.0 2.5
## 2 20.1 2.5
## 3 20.2 2.5
## 4 20.3 2.5
## 5 20.4 2.5
## 6 20.5 2.5
tail(params)
##        mu  sd
## 8176 29.5 4.5
## 8177 29.6 4.5
## 8178 29.7 4.5
## 8179 29.8 4.5
## 8180 29.9 4.5
## 8181 30.0 4.5

Now you can calculate the unnormalized posterior for all parameter combinations. It is easiest to use a for loops for this.

Store the result for each parameter combination in the the variable UP. Set “eval = TRUE” at the beginning of the code block, so tha the code is executed when you knit the document!

# Here comes your code in which you calculate the 'likelihood * prior expression'
# (see the formula or your code above) for each combination of parameter values.
params$UP = NA
for (k in 1:nrow(params)) {
  params$UP[k] = 
}

Calculate the evidence

Use the data.frame params that you just created, specifically one column in it, to calculate the evidence. If you are uncertain, check the equation above

# Calculate the evidence from the unnormalized posterior

Normalize to get a proper posterior distribution.

Now you can use the evidence to add a one more column to the params data frame, let’s call it PP, which has the posterior probability for each parameter combination.

# Calculate posterior probabilities

Plot the posterior

Make a plot of the posterior

# Plot posterior probabilities
# The easiest is to use ggplot here, following this pattern
# ggplot(params, aes(x = mu, y = sd, z = PP)) + geom_contour_filled()
# which assumes that you have called you column for posterior
# probabilities "PP"
library(ggplot2)

Grid approximation, version 2

Here is a template to generate a grid of parameters. Your main task here is to add

  • the upper and lower boundaries for the parameter values you want to explore in the grid
  • the code to calculate the unormalized (or standardized) posterior.

You have to set eval = TRUE at the beginning of the code block so that the code is run when you knit the document.

N_grid_points <-100

# Make a set of grid values for the two parameters
# replace the __ with values of your choice
my_mu<-seq(from=__, to=__, length.out=N_grid_points)
my_sigma<-seq(from=0, to=__, length.out=N_grid_points)

# Make an empty placeholder for the posterior
posterior_unstandardized<-matrix(NA,ncol=N_grid_points, nrow=N_grid_points)

# Calculate the posterior in each grid point
for (i in 1:N_grid_points)
  for (j in 1:N_grid_points)
  {
    posterior_unstandardized[i,j]<- # your code comes here
  }

Next, you need to calculate the evidence and standardize the posterior. Don’t forget to change “eval = FALSE” to “eval = TRUE”.

posterior =

And here is some code to plot the results

Plot the results. Don’t forget to change “eval = FALSE” to “eval = TRUE”.

# Plot (2D) image of the posterior
image(my_mu,my_sigma,posterior,  axes=FALSE);
axis(1, at = seq(20, 30, by = 2))
axis(2, at = seq(0, 5, by = 1))

# Plot (3D) histogram
library(plot3D)
hist3D(z=posterior, border="black")