Previously, we introduced Bayesian Inference with R using the Markov Chain Monte Carlo (MCMC) techniques. The first set of exercises gave insights on the Bayesian paradigm, while the second set focused on well-known sampling techniques that can be used to generate a sample from the posterior distribution .
We saw in the first set of exercises that, for simple problems, the posterior distribution can be found analytically from the prior and the likelihood (e.g. beta distribution). Therefore, we will use distributions where this is not so easy and MCMC sampling is needed.
We will use concepts of the Extreme Value Theory (EVT) with the Generalized Extreme Value distribution and focus on the Gumbel distribution, which is a distribution with only 2 parameters: location and scale.
Don’t worry: the goal is not to study EVT here, but to use it for Bayesian inference. We will use the
evd package with the
dgev() family of functions (
pgev(), etc). Note that when the shape parameter is set to 0, you retrieve the Gumbel distribution.
The solutions for this set of exercises can be found here.
Consider a Gumbel distribution with location and scale parameters set to 30 and 1.5 (resp.) and generate a random sample size of 1,000. Plot the kernel density of this sample together with the theoretical density. What do you see?
seed 1234 to obtain the same results as the solutions.)
Relying on the
dgev() function from the
evd package, write the log-likelihood function of the Gumbel distribution. Then, compute the log-likelihood value for the sample of Exercise 1.
(Note: Working with log-values will simplify the computations.)
Consider a non-informative prior in terms of a large variance normal prior for each of the two parameters, individually centered on the true parameter value. Take an SD of 10. Consider the log-priors, then sum up the two individual log-priors to obtain the “total” log-prior.
Using the code you wrote for Exercises 2 and 3, write a function that computes the log-posterior.
(Note: this is then equal to the sum of the log-prior and the log-likelihood.)
Compute the log-posterior value of the generated sample.
Now you have all the tools to build the Metropolis algorithm.
a. Take the mean and the SD of the generated sample of Exercise 1 as starting values for the location, and the scale parameters (resp.). Take 10,000 iterations.
b. Initialize the output matrix, which will contain the generated chains for the parameters.
c. Run the algorithm by taking normal distributions for proposals (these are symmetric distributions) and take 0.6 and 0.6 as respective standard deviations.
Note: When the proposal distribution is not symmetric, the sampler will be named Metropolis-Hastings algorithm. Metropolis algorithm is a special case of the Metropolis-Hastings.
Check quickly if the chains look stationary, and state whether the Metropolis sample has (potentially) converged or not.
The acceptance rate of the Metropolis algorithm (i.e. the proportion of times that a proposal is accepted in the algorithm) is crucial for the efficiency of the sampler. Complete the previous code in order to get the acceptance rate of this Metropolis sampler.
Now, finish this series of exercises by gathering everything you did
earlier into one function.
BONUS: add a parameter for the burn-in inside the function.
It is recommended to target an acceptance rate of 20%.
Then, use this function to find better values for the proposal standard deviations.