Hamiltonian Monte CarloStan

  • Published on
    23-Aug-2014

  • View
    531

  • Download
    0

DESCRIPTION

Hamiltonian Monte CarloStan

Transcript

Hamiltonian Monte Carlo Stan 2014-07-28 4 Hamiltonian Monte Carlo (HMC)! Stan Hamiltonian Monte Carlo MCMC (Markov chain Monte Carlo) Hybrid Monte Carlo MCMC + deterministic simulation method! MCMC Metropolis-Hastings 0 +/ HMC 0 +/ http://mc-stan.org/ Stan Sampling through adaptive neighborhoods (Gelman et al. 2013) Stanisaw Ulam (19091984) http://en.wikipedia.org/wiki/Stanislaw_Ulam#mediaviewer/File:STAN_ULAM_HOLDING_THE_FERMIAC.jpg Stan 20128 1.0 201212 1.1 20133 1.2 201310 2.0 201312 2.1 20142 2.2 20146 2.3 20147 2.4 Stan 3 RStan! R! PyStan! Python! CmdStan! Stan Hamiltonian Monte Carlo, No-U-Turn Sampler C++ BUGS HMC int real vector row_vector matrix *1010.3 BUGS model! {! for (i in 1:N) {! Y[i] ~ dbin(q[i], 8)! logit(q[i]) # read data! > d ! 1010.3 > model > data ! > fit print(fit, digits = 2)! Inference for Stan model: model.! 3 chains, each with iter=10100; warmup=100; thin=10; ! post-warmup draws per chain=1000, total post-warmup draws=3000.! ! mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat! beta 0.04 0.01 0.33 -0.60 -0.18 0.04 0.27 0.69 3000 1! s 3.04 0.01 0.37 2.41 2.78 3.02 3.27 3.80 2807 1! lp__ -443.66 0.18 9.50 -463.66 -449.84 -443.33 -437.14 -425.83 2805 1! ! Samples were drawn using NUTS(diag_e) at Tue Jun 17 15:45:43 2014.! For each parameter, n_eff is a crude measure of effective sample size,! and Rhat is the potential scale reduction factor on split chains (at ! convergence, Rhat=1). > traceplot(fit) OpenBUGS 1010.3 chain: 3 iteration: 10000 burn-in (warmup): 100 thin: 10 Mac Pro (2.8 GHz Quad-Core Intel Xeon) OS X 10.9.3 OpenBUGS 3.2.2 RStan/CmdStan 2.2.0 OpenBUGS OpenBUGS RStan CmdStan 0 75 150 225 300 compile run *OpenBUGSWine LinuxLinuxOpenBUGSStan (2) Zero-Inated Negative Binomial Model 0 1 2 3 4 5 6 7 8 9 11 13 15 17 y count 0 10 20 30 40 0 : , : p(y| , , ) = (1 )+ NegBin(0| , ) (y = 0) NegBin(y| , ) (y > 0) data {! int N;! int y[N];! }! parameters {! real theta;! real alpha;! real beta;! }! model {! // priors! theta ~ uniform(0, 1);! alpha ~ uniform(0, 100);! beta ~ uniform(0, 100);! for (i in 1:N) {! if (y[i] == 0) {! increment_log_prob(log((1 - theta)! + theta * (beta / (beta + 1))^alpha));! } else {! increment_log_prob(bernoulli_log(1, theta)! + neg_binomial_log(y[i], alpha, beta));! }! }! }! generated quantities {! real mu;! mu