72 Markov chain Monte Carlo algorithms for Mixture Models - M4L1HW2

Caution

Section omitted to comply with the Honor Code

--- title : 'Markov chain Monte Carlo algorithms for Mixture Models - M4L1HW2' subtitle : 'Bayesian Statistics: Mixture Models' categories: - Bayesian Statistics keywords: - Mixture Models - Homework --- ::::: {.content-visible unless-profile="HC"} ::: {.callout-caution} Section omitted to comply with the Honor Code ::: ::::: ::::: {.content-hidden unless-profile="HC"} ## HW - Simulation of Lifetime of Fuses ::: {.callout-note} ### Instructions \index{dataset!fuses} Data on the lifetime (in years) of fuses produced by the ACME Corporation is available in the file `data/fuses.csv`. $$ f(x)=w \lambda \exp \left\{−λx\right\}+(1−w) \frac{1}{\sqrt{2\pi}\tau x} \exp\left\{− \frac{1}{2τ^{2}} (\log(x)−μ)^{2} \right\}, \quad x > 0. $$ {#eq-mixture} The first component, which corresponds to an exponential distribution with rate $\lambda$, is used to model low-quality components with a very short lifetime. The second component, which corresponds to a log-Gaussian distribution, is used to model normal, properly-functioning components. You are asked to modify the implementation of the MCMC algorithm contained in the Reading "Sample code for MCMC example 1" so that you can fit this two-component mixture distributions instead. You then should run your algorithm for 10,000 iterations after a burn-in period of 1,000 iterations and report your estimates of the posterior means, rounded to two decimal places. Assume the following priors: $w∼Uni[0,1], λ∼Exp(1), μ∼Normal(0,1)$ and $τ^2∼IGam(2,1)$. ::: ::: {.callout-note} ### Grading overview The code you generate should follow the same structure as "Sample code for MCMC example 1". In particular, focus on a Gibbs sampler that alternates between the full conditionals for $\omega$, $\lambda$,$\mu$,$\tau^2$ and the latent component indicators $c_1,...,c_n$ Peer reviewers will be asked to check whether the different pieces of code have been adequately modified to reflect the fact that 1. parameters have been initialized in a reasonable way, 2. each of the two full conditional distributions associated with the sampler are correct, and 3. the numerical values that you obtain are correct. To simplify the peer-review process, assume that component 1 corresponds to the exponential distribution, while component 2 corresponds to the log-Gaussian distribution. ::: full conditional for $w_i$ is given by $$ w_i \mid \cdots \sim Beta(1 + \sum_{i=1}^n \mathbb{I}(c_i=1), 1 + n - \sum_{i=1}^n \mathbb{I}(c_i=1)) $$ full conditional for $lambda_i$ is given by $$ \lambda_i \mid \cdots \sim Gamma\left(\sum_{j=1}^n \mathbb{I}(c_j=1), \sum_{j=1}^n x_j \mathbb{I}(c_j=1)\right) $$ Full conditional for $\mu$ $$ \mu \mid \cdots \sim N\left(\frac{\sum_{i:c_i=2} x_i}{{\frac{m(c)}{τ^2} + 1}} + 0, \frac{1}{\frac{m(c)}{τ^2} + 1}\right) $$ Full conditional for $\tau$ $$ \tau^2 \mid \cdots \sim IGam(2+n-m(c), 1+ \frac{1}{2} \sum_{i:c_i=2} (logx_i - μ)^2) $$ ```{r} #| label: fig-load-data #| fig-cap: "Histogram of the fuses failure times" # Load the nest size data rm(list=ls()) library(MCMCpack) set.seed(81196) # So that results are reproducible (same simulated data every time) fuses <- read.csv("data/fuses.csv",header=FALSE) #colnames(fuses) <- c("n") x <- fuses$V1 n <- length(x) # Number of observations # how many rows in the data nrow(fuses) # how many zeros in x sum(x==0) # almost half of the data is zeros! par(mfrow=c(1,1)) xx.true = seq(0, max(x), by=1) hist(x, freq=FALSE, xlab="Fuses", ylab="Density", main="Empirical distribution of fuses failure times") ``` ```{r} #| label: lbl-fuse-mix-1 # MCMC algorithm for fitting a mixture of exponential and log-normal distributions # The actual MCMC algorithm starts here ## MCMC iterations of the sampler iterations <- 11000 burn <- 1000 ## Initialize the parameters, latent vars & Priors KK = 2 # Number of components w = 0.1 # fewer zeros w = rbeta(1, 1, 9) # Uniform prior on w alpha = 2 beta = 1 mu = mean(log(x)) tau = sd(log(x)) lambda = 20 / mean(x) aa = c(1, 1) # Uniform prior on w prob = aa/sum(aa) # latent indicators cc = sample(1:KK, n, replace = TRUE, prob = prob) # Storing the samples cc.out = array(0, dim=c(iterations, n)) w.out = numeric(iterations) lambda.out = numeric(iterations) tau.out = numeric(iterations) #mu.out = array(0, dim=c(iterations, KK)) mu.out = numeric(iterations) logpost = rep(0, iterations) # MCMC iterations for (s in 1:iterations) { # Full conditional for cc v = rep(0,2) for(i in 1:n){ v[1] = log(w) + dexp(x[i], lambda, log=TRUE) v[2] = log(1-w) + dlnorm(x[i], mu, tau, log=TRUE) v = exp(v - max(v))/sum(exp(v - max(v))) cc[i] = sample(1:2, 1, replace=TRUE, prob=v) } #print(cc) # Sample the weights # Full conditional for w w = rbeta(1, 1+sum(cc==1), 1+sum(cc==2)) # Full conditional for lambda lambda = rgamma(1, 1 + sum(cc==1), 1 + sum(x[cc==1])) # Full conditional for mu mean.post = (sum(log(x[cc==2]))/tau^2 + 0)/(sum(cc==2)/tau^2 + 1) std.post = sqrt(1/(sum(cc==2)/tau^2 + 1)) mu = rnorm(1, mean.post, std.post) # Full conditional for tau #tau = sqrt(1/rgamma(1, 2 + sum(cc==2)/2, 1 + 0.5*sum((log(x[cc==2]) - mu)^2))) tau2 = 1/rgamma(1, 2 + sum(cc==2)/2, 1 + 0.5*sum((log(x[cc==2]) - mu)^2)) tau = sqrt(tau2) # Store samples cc.out[s,] = cc w.out[s] = w lambda.out[s] = lambda mu.out[s] = mu tau.out[s] = tau } ``` The posterior means, rounded to two decimal places, are $E{w} \approx 0.10$, $E{λ} \approx 2.29$, $E{μ} \approx 0.79$ and $E{τ} \approx 0.38$. ```{r} #| label: lbl-posterior-means # Posterior summaries w.post = w.out[-(1:burn)] lambda.post = lambda.out[-(1:burn)] mu.post = mu.out[-(1:burn)] tau.post = tau.out[-(1:burn)] cat("Posterior mean of w:", round(mean(w.post),2), "\n") cat("Posterior mean of lambda:", round(mean(lambda.post),2),"\n") cat("Posterior mean of mu:", round(mean(mu.post),2),"\n") cat("Posterior mean of tau:", round(mean(tau.post),2),"\n") ``` add unit tests for the posterior means ```{r} #| label: lbl-test-posterior-means #| error: true library(testthat) testthat::test_that("Posterior means are correct", { expect_equal(round(mean(w.post),2), 0.10) expect_equal(round(mean(lambda.post),2), 2.29,tolerance = 0.15) expect_equal(round(mean(mu.post),2), 0.79) expect_equal(round(mean(tau.post),2), 0.38) }) ``` :::::