STAT314/461 Semester 2, 2021
Assignment 4: Probability revision; rejection sampling;
importance sampling; Metropolis algorithm
Instructions:
1. Total Marks: This assessment will be worth up to 10% of your final
grade (7.5% for STAT461 students). Note that although the marks for
the questions below total to 25 the maximum mark that can be earned
for the assignment is 10. It is up to you which questions you submit for
assessment. Keep in mind the exam at the end of the course – working
through these questions is good practice / revision even if they don’t
directly contribute to your grade for the assignment. You are welcome
to submit all the questions
2. Due date: Due 11:55pm Wednesday October 6th 2018.
3. Your R code should be included with the submission as a separate
file(s). So if you are using R markdown to prepare your assignment
please include the .Rmd file as well as the html or pdf produced by
knitr. Otherwise include your .r file along with a pdf containing text
and graphics.
4. While you are welcome to discuss the assignment with other 314/461
students you must ultimately submit a document that persuades the
marker that you, individually, understand the material.
1. Probability practice
Suppose X and Y are discrete random variables that are exchangeable.
That is
Pr(X = a, Y = b) = Pr(X = b, Y = a),∀a, b.
Show that the marginal distributions of X and Y are identical. [ 2
marks]
( This may seem a bit obscure but we will find this result useful when
proving the convergence of the Metropolis-Hastings algorithm )
2. Likelihood and Rejection sampling Consider the simple but non-
conjugate model for the births data given in the R code contained in
births2021_assignment4.R
1
Yi
indep∼ Poisson(θNi), i = 1, . . . 16
θ ∼ logNormal(−4.5, 0.59) (1)
where Yi is the number of births for women aged 15-19 in the i
th region
and Ni is the number of women aged 15-19, in the i
th region. We
regard the population sizes as known constants. Under this model
E(Yi/Ni|θ) = θ, i = 1, . . . 16 That is, under the model, the expected
birth rate is identical for each of the regions and is equal to θ. The
parameter θ is sometimes referred to as the underlying rate. Under the
model the expected number of births in the ith region is E(Yi|θ) = θNi,
which varies over the regions on account of variations in population
size. The Poisson probability mass function in this case is
Pr(Yi = y|θ) = (θNi)y exp(−θNi)/y!
which can be computed in R with the dpois() function using
Pr(Yi = y|θ) = dpois(x = y, lambda = (θNi)),
with the addition of the
log=TRUE
option returning the log of the probability mass function.
logNormal() refers to the log−normal distribution which can be ac-
cessed using the lnorm family of functions in R; e.g. rlnorm(n =
1000,meanlog = 1, sdlog = 1) generates 1000 values from the log-
normal distribution with parameters (1, 1).
(a) Write a function in R to compute the log-likelihood, for any value
of θ based on the births data, and plot the log-likelihood for a
range of θ values spanning the MLE. [3 marks].
(Hint: the log-likelihood will be the sum of log-likelihood contri-
butions from each of the 16 regions, so you will need to include
a sum() statement at the end of your function in order to ensure
the overall log-likelihood for the data is returned). It would be
good practice to make the parameter θ the first argument of the
function.
(b) Write a R program to implement rejection sampling to approx-
imate the posterior for θ, using the prior as the approximating
density [3 marks]
(c) Write a R program to implement rejection sampling to approxi-
mate the posterior for θ, using a Gamma density as the approxi-
mating density [5 marks]
2
(d) Plot the posterior density for θ based on the sample obtained from
at least one of you rejection sampler methods with the approxi-
mating density shown also. [2 marks]
3. Importance sampling
(a) Write a progam to implement an importance sampler to simulate a
t3(0, 1) density, using a N(0, 1) as the approximating density. [3 marks]
(b) Compute the effective Monte Carlo sample size Neff [2 marks]
(c) Contrast the importance sampling estimates of mean, standard devia-
tion, and proportion less than -2.35 and contrast these with the true
values of the corresponding features of the t3(0, 1) distribution [2 marks]
(d) In light of your answers to (b) and (c) comment on the adequacy of
the choice of a N(0, 1) as the approximating density for the importance
sampler. [1 mark]
4. Metropolis algorithm for posterior inference
Write a R program to implement a Metropolis algorithm to obtain a
sample from the posterior distribution of the parameters of a logistic
regression model fitted to data from the Income Survey. You can use
the binary outcome Y = (income > $1, 250 per week) discussed in
lectures or another variable of your choosing. However, your model
should include at least two predictors and “make sense.” For example
I don’t want to see models predicting male sex from age and region.
A shell program has been provided to help get you going. This code is
in:
MH_examples2021_logisticregression_shell.r
This program reads in the Income Survey data and functions for com-
puting the log-likelihood and log- unnormalised posterior for a logistic
regression model as well as functions to help with setting priors for
logistic regression model parameters. In the main code you will need
to replace they symbol “…” with appropriate choices. In the course of
so-doing you will
(a) specify a proper prior for the parameters of your model.
(b) make clear how you choose your initial values
(c) specify your jumping distribution (aka proposal distribution)
3
(d) decide whether your Metropolis algorithm has converged to the
target posterior distribution and produce evidence that the algo-
rithm has in fact converged.
(e) produce some posterior inferences for one parameter of interest in
your model: posterior density plot for the odds ratio correspond-
ing to the parameter of interest; 90% credible interval for the odds
ratio; the posterior probability that the odds ratio exceeds one.
[10 marks]
4