<- 10L # number of observations
nobs <- 8 # average waiting time
ybar <- 1000L # number of draws
B # Un-normalized log posterior: scaled log likelihood + log prior
<- function(x){
upost dgamma(x = x, shape = nobs + 1L, rate = nobs*ybar, log = TRUE) +
log(2) + dt(x = x, df = 1, log = TRUE)}
<- rust::ru(logf = upost,
post_samp n = B,
d = 1, # dimension of parameter (scalar)
init = nobs/ybar)$sim_vals # initial value of mode
Exercise 3
Exercise 3.1
Consider the Laplace family of distribution, \(\mathsf{La}(\nu, \tau)\), with density \[\begin{align*} g(x; \nu, \tau) = \frac{1}{2\tau} \exp\left(- \frac{|x-\nu|}{\tau}\right), \qquad \nu \in \mathbb{R}, \tau > 0 \end{align*}\] as a candidate distribution for rejection sampling from \(\mathsf{No}(0,1)\).
- Provide an inversion sampling algorithm to generate from \(\mathsf{La}(\nu, \tau)\).
- Can you use the proposal to generate from a standard Gaussian? for Student-\(t\) with 1 degree of freedom? Justify your answer.
- Consider as proposal a location-scale version of the Student-t with \(\nu=3\) degrees of freedom. Find the optimal location and scale parameters and the upper bound \(C\) for your choice.
- Use the accept-reject to simulate 1000 independent observations and compute the empirical acceptance rate.
Exercise 3.2
We revisit Exercise 1.3, which used a half-Cauchy prior for the exponential waiting time of buses.
The ratio-of-uniform method, implemented in the rust
R package, can be used to simulate independent draws from the posterior of the rate \(\lambda\). The following code produces
Estimate using the Monte Carlo sample:
- the probability that the waiting time is between 3 and 15 minutes
- the average waiting time
- the standard deviation of the waiting time
Next, implement a random walk Metropolis–Hastings algorithm to sample draws from the posterior and re-estimate the quantities. Compare the values.
Exercise 3.3
Repeat the simulations in Example 3.6, this time with a parametrization in terms of log rates \(\lambda_i\) \((i=1,2)\), with the same priors. Use a Metropolis–Hastings algorithm with a Gaussian random walk proposal, updating parameters one at a time. Run four chains in parallel.
- Tune the variance to reach an approximate acceptance rate of 0.44.
- Produce diagnostic plots (scatterplots of observations, marginal density plots, trace plots and correlograms). See
bayesplot
orcoda
. Comment on the convergence and mixing of your Markov chain Monte Carlo. - Report summary statistics of the chains.