Department of Mathematics and Statistics STAT3405/STAT4066
Important: This assignment is assessed. Your work for this assignment must be submitted
via LMS by 5pm on Friday, 16 April 2021.
The expectation for a submission are:
• The questions are answered in complete sentences. Marks will be awarded for the cor-
rectness of the answers and that they are given in complete sentences.
The answers to the questions should be submitted via LMS.
• That numerical answers are rounded to an appropriate number of digits.
In particular, R code developed for Task 1(a) should be submitted as an attachment in
your LMS submission. But so should be the BUGS or Stan code that you use to solve
Task 2
Submission of code must be in plain text format, do not use any binary format.
Marks will be awarded for
– the correctness of the code, i.e. how easy it is to read it1 ; and
– how easy it is to get the code to run, if necessary.
You may receive comments on the efficiency of your code, but there are no marks for
efficiency.
Unless special consideration is granted, any student failing to submit work by the deadline
will receive a penalty for late submission (as described in the unit outline).
Plagiarism: You are encouraged to discuss assignments with other students and to solve
problems together. However, the work that you submit must be your sole effort (i.e. not
copied from anyone else). If you are found guilty of plagiarism you may be penalised. You
are reminded of the University’s policy on ‘Academic Conduct’ and ‘Academic Misconduct’
(including plagiarism):
http://www.student.uwa.edu.au/learning/resources/ace/conduct
Various material at the following URL might be helpful too:
https://www.uwa.edu.au/students/study-success/studysmarter
1You might want to have a read of http://adv-r.had.co.nz/Style.html or https://mc-stan.org/docs/
2_22/stan-users-guide/stan-program-style-guide.html for discussions on how to write code in a readable
manner, if you have never come across style guides before.
Semester 1, 2021, page 1 Set 1 Due date: Friday, 2021-04-16, 5:00pm
Department of Mathematics and Statistics STAT3405/STAT4066
Task 1. This task explores the Metropolis-Hastings (MH) algorithm used in the second com-
puter lab a bit further. In that lab, we used an MH algorithm with a proposal distribution that
was independent of the current location of the chain. While such MH algorithm can be used,
they are typically quite inefficient if the proposal distribution is not “close” to the distribution
one wants to sample from.
Another popular way of implementing an MH algorithm is to use a random walk MH
algorithm. In such an implementation the proposed new state is obtained by adding some
random noise (typically using a normal distribution with mean zero) to the current state. In
other words, the proposal distribution is a normal distribution with mean equal to the current
state and some variance (which serves can serve as a tuning parameter for the MH algorithm).
The difficulty with using a random walk MH algorithm in the example used in the computer
lab is that we are modelling a probability, which is always between 0 and 1. Adding random
noise to the current state could lead to a proposal that is outside the interval [0, 1]! We can
circumvent this difficulty by first transforming the current state (a value between 0 and 1) via
the logit transformation to the corresponding log-odds (a value between −∞ and ∞). We can
then add some noise to the transformed current state to obtain a proposal on the log-odds scale
and then transform this proposal back to the interval [0, 1].
Specifically, the logit transformation is the function logit(p) = log
(
p
1−p
)
which is readily
implemented in R:
logit <- function(p) log(p/(1-p))
The inverse of the logit transformation is the function h(lo) = 11+exp(−lo) , which takes a value
on the real line (some log-odds) and transforms it back to the interval [0, 1]. This function is
also easily implemented in R:
invlogit <- function(lo) 1/(1+exp(-lo))
Finally, it can be shown that the proposal distribution for this process (i.e. (1) transform-
ing the current state p(t) to the log-odds scale to obtain, say, u(t), (2) adding random noise
coming from a normal distribution with mean 0 and variance σ2 to u(t) to obtain, say, w(t),
and (3) transforming w(t) back to (0, 1) using the inverse logit transformation to obtain the
proposal y(t)) has the following PDF:
q(y|p) =
1√2piσ2 1y(1−y) exp
(
− 1
2σ2
[
log
(
y
1−y
)
− log
(
p
1−p
)]2)
if 0 < y < 1 and 0 < p < 1
0 otherwise
While this looks impressive, the corrective term for using this particular non-symmetric pro-
posal distribution in the MH algorithm is actually quite simple:
q
(
p(t−1)|y(t))
q
(
y(t)|p(t−1)) = y(t)
(
1− y(t))
p(t−1)
(
1− p(t−1))
Now
(a) Implement the random walk MH algorithm as described above. Appropriate code for this
part of the task should be submitted.
(b) What is the acceptance rate of your MH algorithm if you use σ = 1100 (i.e σ
2 = 110000)?
How would you describe the Markov chain produced by this setting of σ?
Semester 1, 2021, page 2 Set 1 Due date: Friday, 2021-04-16, 5:00pm
Department of Mathematics and Statistics STAT3405/STAT4066
(c) What is the acceptance rate of your MH algorithm if you use σ = 10? How would you
describe the Markov chain produced by this setting of σ?
(d) For what value of σ do you obtain an acceptance rate of around 70%? How does the
corresponding Markov chain look like?
Task 2. In the sixth lecture we discussed data from a survey was done of bicycle and other
vehicular traffic in the neighbourhood of the campus of the University of California, Berkeley.
The observed data from this survey for fairly busy streets with a bike route are:
# cycles (yi) 8 35 31 19 38 47 44 44 29 18
# other vehicles 29 415 425 42 180 675 620 437 47 462
# total (mi) 37 450 456 61 218 722 664 481 76 480
For these data:
(a) Fit a model that models yi, i = 1, . . . , 10, as being a realisation from a binomial dis-
tributed random variable Yi with parameters mi and p. Regard the mis as fixed and
given, and put your favourite non-informative prior on p.
What is your Bayesian estimate for p?
(b) Fit a model that models yi, i = 1, . . . , 10, as being a realisation from a binomial dis-
tributed random variable Yi with parameters mi and pi. Regard the mis as fixed and
given, and put your favourite non-informative prior on the pis.
What are your Bayesian estimate for the pis?
For this model, also look at the range of the pis, i.e. study the parameter
r = max
i=1,...,10
pi − min
i=1,...,10
pi
Include a plot of the posterior density of r in your submission. What is your Bayesian
estimate for r? What is the standard deviation of the posterior distribution of r?
(c) Given your results in part (b), which model do you think is more appropriate for these
data? Discuss briefly (i.e. less than a paragraph).
Semester 1, 2021, page 3 Set 1 Due date: Friday, 2021-04-16, 5:00pm
学霸联盟