--- title: "ABC and Convergence of MCMC" author: "Jacob Carey" date: \today output: beamer_presentation: toc: true theme: "Szeged" slide_level: 2 fig_caption: no --- # Approximate Bayesian Computation ## Posterior Approximation >- Likelihood: $P(X | \theta)$ >- Prior: $P(\theta)$ >- Posterior: $P(\theta | X) \propto P(X | \theta) P(\theta)$ ## Example - Normal Model with Known Variance - Likelihood: $X \sim \text{Normal}(\mu)$ ($\sigma^2$ known) - Prior: $\mu \sim \text{Normal}(\mu_0, \sigma^2_0)$ - Posterior: $P(\mu | X)$ ## Metropolis-Hastings (Type of MCMC Algorithm) >- Generic (nearly) always works instance of MCMC >- Easy to implement >- Draw $\theta_{t+1} \sim F(\theta | \theta_t, \delta)$ ($F$ is proposal distribution) >- Accept $\theta_{t+1}$ with probability $\alpha=P(X | \theta_{t+1})P(\theta_{t+1})/(P(X | \theta_t)P(\theta_t))$ ## MCMC Overview - Draw $S$ samples from $P(\theta | X)$ - Repeat $L$ times (creating $L$ chains of $S$ samples each) - Compare each set of $S$ samples to check convergence (next section!) - Combine posterior samples from chains, determine posterior summaries (e.g. mean of samples for posterior mean), publish ## Likelihood >- Sometimes, likelihood is computationally difficult >- Sometimes, it's not really possible to formulate a likelihood >- This is where Approximate Bayesian Computation (ABC) can help >- Also called Likelihood-free method >- Note: is generally slower than MH ## ABC Alg. 1 - Given observation $y$, repeat until $S$ points have been accepted >- Draw $\theta_i \sim p(\theta)$ >- Simulate $X_i \sim p(x | \theta_i)$ >- Accept $\theta_i$ if $X_i = y$ ## ABC Alg. 2 - Given observation $y$, repeat until $S$ points have been accepted >- Draw $\theta_i \sim p(\theta)$ >- Simulate $X_i \sim p(x | \theta_i)$ >- Accept $\theta_i$ if $\rho(S(X_i), S(y)) < \epsilon$ ## ABC Alg. 2 Example See code ## ABC Alg. 3 (MCMC-ABC) >- Draw $\theta_{t+1} \sim F(\theta | \theta_t, \delta)$ ($F$ is proposal distribution) >- Simulate $X_i \sim p(x | \theta_i)$ >- If $\theta_i$ if $\rho(S(X_i), S(y)) < \epsilon$ >- Accept $\theta_{t+1}$ with probability $\alpha=P(\theta_{t+1})/P(\theta_t)$ >- If rejected at any point, retain previous value ## ABC Alg. 3 Example See code ## ABC Further - Many improvements that can be made - One of which is assigning a prior to epsilon - Others involve SMC - See review by Beaumont (2010) for more information # Convergence ## MCMC Overview - Draw $S$ samples from $P(\theta | X)$ - Repeat $L$ times (creating $L$ chains of $S$ samples each) - Compare each set of $S$ samples to check convergence (this section!) - Combine posterior samples from chains, determine posterior summaries (e.g. mean of samples for posterior mean), publish ## Gelman-Rubin Diagnostic (Estimated Potential Scale Reduction Factor) - Compute m independent Markov chains - Compares variance of each chain to pooled variance - Provides estimate of how much variance could be reduced by running chains longer http://astrostatistics.psu.edu/RLectures/diagnosticsMCMC.pdf ## R implementation (Convergence) {r gelm_diag} library(coda); set.seed(1) x <- mcmc.list(mcmc(rnorm(1000)), mcmc(rnorm(1000)), mcmc(rnorm(1000))) gelman.diag(x)  ## R implementation (Not converged) {r gelm_bad} x <- mcmc.list(mcmc(rnorm(1000, 1)), mcmc(rnorm(1000, 2)), mcmc(rnorm(1000, 1))) gelman.diag(x)  ## Gelman-Rubin with CoGAPS {r cogaps, message=FALSE, warning=FALSE, results='hide'} library(CoGAPS) data(SimpSim) A <- vector("list", 3) atoms <- vector("list", 3) for (l in 1:3) { out <- gapsRun(SimpSim.D, SimpSim.S, nFactor=3, numSnapshots=1000) A[[l]] <- out$ASnapshots atoms[[l]] <- mcmc(out$atomsASamp) }  {r cogaps-clean, echo=FALSE} atoms_mcmc <- mcmc.list(atoms) A_11 <- mcmc.list(mcmc(A[][1, 1, ]), mcmc(A[][1, 1, ]), mcmc(A[][1, 1, ]))  ## Checking Convergence {r convergence} gelman.diag(atoms_mcmc) gelman.diag(A_11)