STAT3405/STAT4066 Department of Mathematics and Statistics

Task 1. This task explores the Metropolis-Hastings (MH) algorithm used in the second com- puter lab a bit further.  However, for this assignment, we will use the likelihood based on the data, discussed in the sixth lecture, from a survey that was done on bicycle and other vehicular traffic in the neighbourhood of the campus of the University of California, Berkeley.  As in the computer lab, we will use Zellner’s prior:

f(p) ∝ pp (1 − p)1−p

for the parameter p in our binomial model for the survey data..

In the second lab, we used an MH algorithm with a proposal distribution that was inde- pendent 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 algorithm2 . 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).

An obvious difficulty arises when we try to use a random walk MH algorithm when 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 a transformation to a value

p(b)r(e)o(t)p(w)o(e)s(e)a(n)l on thi(−∞ a)s(n)t(d)ran(∞.)sform(We c)e(a)d(n)scale and th(then add so)en(me)tr(n)a(o)n(is)sf(e)o(t)r(o)m(t)t(h)h(e)is(tr)proposa(ansform)l(e)back to th(d current s)e(t)i(a)nt(te)erva(to o)l(b)0(a)1(n)

Specifically,  in this exercise we will use the probit  transformation,  which  is  the inverse function of the cumulative distribution function (CDF) of the standard normal distribution:

probit(p) = Φ−1 (p).  This transformation is readily implemented in R:

probit  <- function(p) qnorm(p)

The inverse of the probit transformation is, of course, the CDF of the standard normal distri- bution, h(x) = Φ(x). This function takes as argument a value on the real line and transforms it (back) to the interval (0, 1). This function is also easily implemented in R:

invprobit  <- function(x) pnorm(x)

Finally, it can be shown that the proposal distribution for this process (i.e. (1) transforming the current state p(t−1)  to the real line using the probit transformation to obtain, say, u(t−1) , (2) adding random noise coming from a normal distribution with mean 0 and variance σ 2 to u(t−1)  to obtain, say, w(t), and (3) transforming w(t)  back to (0, 1) using the inverse probit transformation (i.e. the CDF of the standard normal distribution) to obtain the proposal x(t)) has the following PDF:

q(x|p) = 2πσ(1)2 φ(prob(1)it(x)) exp ( − 2σ(1)2   {probit(x) − probit(p)}2 ） o(if)th(0)erwise(< x <) 1 and 0 < p < 1

where

φ(z) = e2 ,    z ∈ R

is the probability density function of the standard normal distribution (implemented as dnorm() in R).

While this formula looks impressive,  the corrective term for using this particular non- symmetric proposal distribution in the MH algorithm is actually quite simple:

φ(probit (x(t),,

φ(probit (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 σ = σ 2  = How would you describe the Markov chain produced by this setting of σ?

(c) What is the acceptance rate of your MH algorithm if you use σ = 3?  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?

Hint: To answer the last three question you may want to look at the trace plot of the chains and also look at a plot of their auto-correlation functions.

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.

Assume we model the observations yi , i = 1, . . . , 10, as being realisations from independent binomial distributed random variables Yi  with parameters mi  and p.  Regard the mis as fixed and given. To finalise the specification of our Bayesian model assume that we use a Beta(α,β) prior for p, the parameter of interest.

(a)  Determine the Bayesian estimates ˆ(p)k that you obtain when selecting the hyper-parameters

B(α)aye(an)s(d)ian(β)ies(n)ti(t)m(he)at(p)e(r)s(io)kaa(s)ga(α)ins(=)t(β)log(=)ly(k) th(=)ip(1)ot(0),h(1)s(2), t(.)o(.) .e(2)s(0)ubmitte(. Produ)d(c)e a plot of the

(b)  Discuss the plot briefly (i.e. less than a paragraph). Which priors Beta(γk,γk) would you call non-informative? Which priors would you consider to be mildly informative? Which would you call strongly informative?