what’s the probability of a deadlock in solitaire?

Birth of MC methods at Los Alamos

Problem

• assume that \(x \sim X\) and that \(X\) is “unwieldy”
• we’d like to know some property \(F(X)\) of \(X\) that can be expressed as an expectation (where \(f(x)\) could be any useful transformation of single values \(x\)):

\[F(X) = \int P(X = x) \ f(x) \ \text{d}x\]

• examples of \(F(X)\) are the mean, median, variance, 95% HDI, or any cumulative probability, such as \(\int_{0.8}^\infty P(X = x) \ \text{d}x\)

dummy

solution: Monte Carlo sampling

• draw samples \(S = x_1, \dots x_n \sim X\) and compute:

\[F(S) = \frac{1}{N} \sum_{i = 1}^N f(x_i)\]

• if the samples are “good samples” (and if \(f\) is not crazy), this approximates the original expectation well:

\[F(S) \sim \mathcal{N}(\mathbb{E}_X, \frac{\text{Var}(f)}{N})\]

example 1: normal distribution

# get density values and samples
xVec = seq(-4, 4, length.out = 5000)     # vector of x-coordinates
myPlot = ggplot(tibble(x = rnorm(5000, mean = 0, sd = 1)), aes(x = x)) +  
  geom_density() + 
  geom_line(aes(x = xVec, y = dnorm(xVec, mean = 0, sd = 1)), color = "firebrick") +
  xlab("x") + ylab("(estimated) probability")
myPlot

example 2: simulating \(p\)-values

k_obs = 7
n_obs = 24
x_reps = 500000

lhs = map_dbl(1:x_reps, function(i) {
  k_hyp = rbinom(1, size = n_obs, prob = 0.5)
  dbinom(k_hyp, size = n_obs, prob = 0.5)
})

lh_obs = dbinom(k_obs, size = n_obs, prob = 0.5)

mean(lhs <= lh_obs) %>% show()

## [1] 0.064272

tibble(iteration = 1:x_reps, 
       p_value = cumsum(lhs <= lh_obs) / 1:x_reps) %>% 
  ggplot(aes(x = iteration, y = p_value)) + geom_line()

Markov chain

Markov property

\[ P(x_{n+1} \mid x_1, \dots, x_n) = P(x_{n+1} \mid x_n) \]

Markov Chain Monte Carlo methods

get sequence of samples \(x_1, \dots, x_n\) s.t.

1. sequence has the Markov property (\(x_{i+1}\) depends only on \(x_i\)), and
2. the stationary distribution of the chain is \(P\) (the target distribution).

dummy

consequences of Markov property

• non-independence -> autocorrelation
  (nuisance)
• easy proof that stationary distribution is \(P\)
  (reassuring)
• computationally efficient
  (great)
• we can work with non-normalized probabilities dummy
  (absolutely awesome)

island hopping

• set of islands \(X = \{x^s, x^b, x^p\}\)
• goal: hop around & visit every island \(x \in X\) proportional to its population \(P(x)\)
  • think: “samples” \(x \sim P\)
• problem: island hopper can remember at most 2 islands’ population
  • think: we don’t know the normalizing constant

Metropolis Hastings

• let \(f(x) = \alpha P(x)\) (e.g., unnormalized posterior)

• start at random \(x_0\), define probability \(P_\text{trans}(x_i \rightarrow x_{i+1})\) of going from \(x_{i}\) to \(x_{i+1}\)
  • proposal \(P_\text{prpsl}(x_{i+1} \mid x_i)\): prob. of considering jump to \(x_{i+1}\) from \(x_{i}\)
  • acceptance \(P_\text{accpt}(x_{i+1} \mid x_i)\): prob of accepting jump to proposal \(x_{i+1}\) \[P_\text{accpt}(x_{i+1} \mid x_i) = \text{min} \left (1, \frac{f(x_{i+1})}{f(x_{i})} \frac{P_\text{prpsl}(x_{i} \mid x_{i+1})}{P_\text{prpsl}(x_{i+1} \mid x_i)} \right)\]
  • transition \(P_\text{trans}(x_i \rightarrow x_{i+1}) = P_\text{prpsl}(x_{i+1} \mid x_i) \ P_\text{accpt}(x_{i+1} \mid x_i)\)

properties of MH

• motto: always up, down with probability \(\bfrac{f(x^{i+1})}{f(x^{i})}\)

• ratio \(\bfrac{f(x^{i+1})}{f(x^{i})}\) means that we can neglect normalizing constants

• \(P_\text{trans}(x^i \rightarrow x^{i+1})\) defines transition matrix -> Markov chain analysis!

• for suitable proposal distributions, it can be shown that:
  • stationary distribution exists (first left-hand eigenvector)
  • every initial condition converges to stationary distribution
  • stationary distribution is \(P\)

7-island hopping

7-island hopping, cont.

influence of proposal distribution

Gibbs sampling: introduction

• MH is a very general and versatile algorithm

• however in specific situations, other algorithms may be better applicable
  • e.g., when \(x \sim P\) is high-dimensional (think: many parameters)

• in such cases, the Gibbs sampler may be helpful

• basic idea: split the multidimensional problem into a series of simpler problems of lower dimensionality

Gibbs sampling: main idea

• assume \(X=(X_1,X_2,X_3)\) and \(p(x)=P(x_1,x_2,x_3)\)

• start with \(x^0=(x_1^0,x_2^0,x_3^0)\)

• at iteration \(t\) of the Gibbs sampler, we have \(x^{i-1}\) and need to generate \(x^{i}\)
  • \(x_1^{i} \sim P(x_1 \mid x_2^{i-1},x_3^{i-1})\)
  • \(x_2^{i} \sim P(x_2 \mid x_1^{i},x_3^{i-1})\)
  • \(x_3^{i} \sim P(x_3 \mid x_1^{i},x_2^{i})\)

• for a large \(n\), \(x^n\) will be a sample from \(P(x)\)
  • \(x^n_1\) will be a sample from the marginal distribution \(P(x_1)\):

\[ P(x_1) = \int P(\tuple{x_1, x_2, x_3}) \ \text{d} \tuple{x_2, x_3}\]

example: 2 coin flips

example: Gibbs jumps

example: Gibbs for 2 coin flips

example: MH for 2 coin flips

summary

• Gibbs sampling can be more efficient than MH

• Gibbs needs samples from conditional posterior distribution
  • MH is more generally applicable

• preformance of MH/Gibbs depends on proposal distribution
  • clever software helps determine when/how to do Gibbs and how to tune the proposal function (next class)

problem statements

convergence/representativeness

• we have samples from MCMC
• in the limit, samples must be representative of \(P\)
• how do we know that our meagre finite samples are representative?

efficiency

• ideally, we’d like the shortest samples that are representative
• how do we measure that we have “enough” samples?

general idea behind many practical solutions

• compare how similar several independent sample chains are

• look at the temporal development of single chains

packages to analyze MCMC chains

require('coda')
require('ggmcmc')

example MCMC data

# self-made MH 
# returns MCMC-list from package 'coda'
out = MH(f, iterations = 60000, 
         chains = 2, burnIn = 10000) 
head(out[[1]])      

## Markov Chain Monte Carlo (MCMC) output:
## Start = 1 
## End = 7 
## Thinning interval = 1 
##          variable
## iteration           mu    sigma
##         1 -0.005970447 1.051339
##         2 -0.005970447 1.051339
##         3 -0.005970447 1.051339
##         4 -0.005970447 1.051339
##         5 -0.005970447 1.051339
##         6 -0.005970447 1.051339
##         7 -0.005970447 1.051339

# cast the same information into 'ggmcmc' format
out2 = ggmcmc::ggs(out)
# same samples but now in a tibble
out2             

## # A tibble: 200,000 x 4
##    Iteration Chain Parameter    value
##        <int> <int> <fct>        <dbl>
##  1         1     1 mu        -0.00597
##  2         2     1 mu        -0.00597
##  3         3     1 mu        -0.00597
##  4         4     1 mu        -0.00597
##  5         5     1 mu        -0.00597
##  6         6     1 mu        -0.00597
##  7         7     1 mu        -0.00597
##  8         8     1 mu        -0.00597
##  9         9     1 mu        -0.00597
## 10        10     1 mu        -0.00597
## # ... with 199,990 more rows

trace plots in ‘ggmcmc’

• a trace plot plots the samples (separately for each chain) in the order in which they were gathered (by iteration)
  • did the chains converge to the same stable range of values?

ggs_traceplot(out2)

visual inspection of convergence

trace plots from multiple chains should look like:

examining sample chains: beginning

examining sample chains: rest

burn-in (= warm-up) & thinning

burn-in (warm-up)

remove initial chunk of sample sequence to discard effects of (random) starting position

thinning

only consider every \(i\)-th sample to remove autocorrelation

R hat

\(\hat{R}\)-statistics,

• a.k.a.:
  • Gelman-Rubin statistics
  • shrink factor
  • potential scale reduction factor
• idea:
  • compare variance within a chain to variance between chains
• in practice:
  • use software to compute it
  • aim for \(\hat{R} \le 1.1\) for all continuous variables

R hat in ‘coda’

gelman.diag(out)

## Potential scale reduction factors:
## 
##       Point est. Upper C.I.
## mu         10.73       66.0
## sigma       3.16       17.8
## 
## Multivariate psrf
## 
## 11.5

R hat in ‘ggmcmc’

ggs_Rhat(out2)

autocorrelation

autocorrelation in ‘coda’

autocorr.plot(out)

autocorrelation in ‘ggmcmc’

ggs_autocorrelation(out2)

effective sample size

• intuition:
  • how many samples are “efficient, actual samples” if we strip off autocorrelation
• definition:

\[\text{ESS} = \frac{N}{ 1 + 2 \sum_{k=1}^{\infty} ACF(k)}\]

effective sample size in ‘coda’

effectiveSize(out)

##      mu   sigma 
## 1361.64 1184.35

outlook

Tuesday

• introduction to Stan

Friday

• 1st practice session with BDA

to prevent boredom

• install Stan & peak into its documentation

• browse Gelman et al (2014) Part III for more on Bayesian computation