Recap and simulation validation

STA602 at Duke University

Goal

The goal of this lab is to (1) comprehensively recap some of the main ideas introduced to date and (2) to think about how you can validate certain ideas via simulation.

Exercise I

  • Read about the Pareto distribution using ?Pareto::rPareto().
    • Note: You might need to install the package first, install.packages("Pareto")
  • Simulate data using rPareto for some fixed \(\alpha\), \(t\).
  • Evaluate the likelihood you derive in exercise 1a) on the exam and compare to Pareto::dPareto() from the package using the simulated data.

Exercise II

Code exercise 1c from the exam.

Exercise III

Gibbs sampling preview

Setup

  • Let \(\sigma^2 = 1\), \(a = 3\), \(b = 3\).
  • For purposes of simulation, let the true vector \(w\) and \(\theta\) be given below:
set.seed(602)
N <- 10
w <- rgamma(N, 3, 3)
thetaTrue <- .3
sigma2 <- 1
  • Simulate 10 data points using the true \(w\) and \(\theta\).
set.seed(602)
y <- rnorm(N, mean = thetaTrue, sd = sqrt(sigma2/w))

Next, alternate sampling from each full conditional density some large number of times, as derived in exercise 2 from the exam. In other words,

  1. sample from \(\theta | \boldsymbol{w}, \boldsymbol{y}\), next
  2. using that new value of \(\theta\), sample \(w_i | \theta, \boldsymbol{w}_{-i}, \boldsymbol{y}\), for \(i\) in \(1,\ldots, 10\)
  3. Repeat, starting at step 1.

Exercise III (part ii)

Plot the empirical density function of the marginal posterior distribution of \(\theta\), i.e. plot a density/histogram plot of \(\{\theta^{(1)}, \ldots, \theta^{(10000)}\}\). Also plot the marginal posterior density of each \(w_i\). Comment on what you see. We’ll get back to why this works in just a few classes.

Exercise IV

  • Come up with a way to empirically verify exercise 2b from the exam. Discuss with your neighbor and implement.