Many problems do not have easy analytical solution
Autocorrelation
State-Space Models
Mixed Models
Solution is simulated draws form the posterior…
Monte Carlo Sampling for Bayesian Estimation
Markov Chain Monte Carlo (MCMC)
The problems of MCMC
Hamiltonian Monte-Carlo (HMC)
HMC v. MCMC
Implementing HMC
King Markov and His Islands
King Markov and His Islands
How to move around Islands
Flip a coin. Heads, move left. Tails, move right.
Figure out the number of people on each island.
Assign a probability, p = # next island / # current island
Choose a random number. If number < p, move.
Rinse and repeat
What Metropolis MCMC Looks Like
mcmc_ex <-function(num_weeks =1e5, current=10, positions =rep(0, num_weeks)){for ( i in1:num_weeks ) {# record current position positions[i] <- current# flip coin to generate proposal proposal <- current +sample( c(-1,1) , size=1 )# now make sure he loops around the archipelagoif ( proposal <1 ) proposal <-10if ( proposal >10 ) proposal <-1# move? prob_move <- proposal/current current <-ifelse( runif(1) < prob_move , proposal , current ) } positions}
Metropolis MCMC in Action: 10 Weeks
Population = Island Number
Metropolis MCMC in Action: 50 Weeks
Population = Island Number
Metropolis MCMC in Action: 1000 Weeks
Population = Island Number
Metropolis MCMC For Models
Each island is a set of parameter choices
Each “population” is a posterior density
The path is a ‘chain’
Note the autocorrelation - we “thin” chains
Only use every ith sample so that there is no autocorrelation
MCMC In Practice for Models
from http://www.cnblogs.com/Nietzsche/p/4255948.html
Monte Carlo Sampling for Bayesian Estimation
Markov Chain Monte Carlo (MCMC)
The problems of MCMC
Hamiltonian Monte-Carlo (HMC)
HMC v. MCMC
Implementing HMC
MCMC is not a Panacea
MCMC is not a Panacea
How can MCMC Fail?
MCMC (particularly Metropolis) can get stuck
High correlation between parameters
Concentration of Measure problem
Start values can still be important
But, who can always choose the right start values?
One way we try and assess is fit with many chains and make sure they converge - This can lead to depression as things never converge, or converge in impossible values
Concentration of Measure Problem
When you sample at high dimensions, as more of the density is away from the mode, you sample further and further from the mode
Concentration of Measure Problem
This is a Real Problem - Perhaps Bigger Than We Knew
MCMC Algorithm Solution
Metropolis MCMC inefficient and prone to getting stuck
Many algorithms to come up with clever proposed moves to speed up
Gibbs sampling used for BUGS, JAGS, etc.
Still has same problems as Metropolis
Or… Abandon search and use more deterministic sampling
Hamiltonian MCMC
Monte Carlo Sampling for Bayesian Estimation
Markov Chain Monte Carlo (MCMC)
The problems of MCMC
Hamiltonian Monte-Carlo (HMC)
HMC v. MCMC
Implementing HMC
King Monty Of the Valley and Hamilton, His Advisor
Moving Around the Valley
Royal carraige pushed in random direction with random momentum.
As it moves up the hill, slows down until it reverses.
After a time, we stop it and record its location
Will spend more time in at bottom of the hill (most people)
arXiv:2006.16194v2
How Monty Samples
A Few Tuning/Algorithmic Details
We need to chose how large a step we have between each time the gradient is recalculated
Automated in the Warmup phase
Warmpup is NOT burn-in - no information
These are Leapfrog steps
We do NOT want our chain to come back to where it started
Honors Stanislaw Ulam, one of the originators of Monte-Carlo computing.
Stan is named after him!
Compiles a model, so it can take a while
Can specify number of chains and other parameters
And now our samples are already part of our model!
Careful, models can get large (in size) depending on number of parameters and samples
Data Prep for HMC
data(rugged)rugged$log_gdp <-log(rugged$rgdppc_2000)#Prep the datad <- rugged[complete.cases(rugged$log_gdp),]# Need indices >0d$cont_africa <- d$cont_africa +1# Only want the data we are using# (otherwise slows STAN down)d.trim <- d[ , c("log_gdp","rugged","cont_africa") ]
The Model…
int_mod <-alist(#likelihood log_gdp ~dnorm(mu, sigma),#Data generating process mu <- bR[cont_africa]*rugged + bA[cont_africa],#priors bR[cont_africa] ~dnorm(0,1), bA[cont_africa] ~dnorm(8,100), sigma ~dcauchy(0,2))
Warning: 1 of 1500 (0.0%) transitions ended with a divergence.
See https://mc-stan.org/misc/warnings for details.
Warning: 1103 of 1500 (74.0%) transitions hit the maximum treedepth limit of 10.
See https://mc-stan.org/misc/warnings for details.
Folk theorem of statistical computing: When you are having trouble fitting a model, it often indicates a bad model.
In many ways, HMC and MCMC failing will help you find issues in your own intuition in both model structure and how to build a good model, rather than just being a PITA.
