Hamiltonsche Mechanik

Die Idee des HMC-Algorithmus kommt aus der physikalischen Statistik. Die in der Statistik benutzte (Radon-Nikodým-)Dichte \(f\) entspricht der physikalischen (Gibbs-)Energie \(\phi\):

\[ f(x) \propto \exp(-\phi(x)) \]

Die Bewegung eines Objekts durch einen Raum lässt sich durch Differentialgleichungen mit Ort \(x\) (potentielle Energie \(U(x)\)) und Impuls \(v\) (kinetische Energie \(K(v)\)) beschreiben:

\[ E(x,v)=U(x)+K(v) \]

• Im HMC-Algorithmus ist das “Objekt” die Zufallsvariable \(x\), dass sich durch seinen Wertebereich bewegt.
• Über den Zusammenhang von Energie und Dichte kann man die Idee auch für Monte Carlo benutzen.

Hamiltonian MC

• Für den HMC-Algorithmus braucht man eine zusätzliche zufällige Hilfsvariable \(v\) (“Geschwindigkeit”)

• Priori von \(v\) (Normalverteilung) ist nötig

• In den Iterationen des HMC-Algorithmus zieht man zufällig aus der Priori

• Vorstellung: Bewegung eines Pucks auf der Parameter-Oberfläche, der immer wieder in zufällige Richtungen geschubst wird.

• Für MCMC müssen Hamilton-Gleichungen zeitdiskretisiert werden, mit Schrittgröße \(\epsilon\)

Vor- und Nachteile von HMC

• Funktioniert nur bei differenzierbaren Funktionen
• Einzelne Iteration ist sehr viel aufwendiger als bei M-H oder Gibbs
• I.d.R. geringere Abhängigkeit, großer Abstand zwischen Ziehungen
• Hohe Akzeptanzrate
• Probleme bei isolierten lokalen Minima (bzw. Maxima)
• Tuning ist notwendig (Schrittweite \(\epsilon\))

Literatur: Radford M. Neal: MCMC using Hamiltonian dynamics, in: Steve Brooks, Andrew Gelman, Galin Jones, and Xiao-Li Meng. Handbook of Markov Chain Monte Carlo. Chapman & Hall / CRC Press, 2011

RStan

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

Modell

• Die Modell-Spezifikation wird in einer Datei bzw. als Text abgespeichert.
• Die Modell-Sprache ähnelt der BUGS/JAGS-Sprache
• Variablen können integer, real oder vector sein
• Für die Variablen können Grenzen angegeben werden: <lower=0> oder <upper=1>, auch <lower=0, upper=1>
• Die Daten werden hier als Bernoulli-verteilt definiert, mit der Logit-Transformation auf die unbakannte Wahrscheinlichkeit

linmodel.stan = "
    data {
      int<lower=0> N;
      int<lower=0> M;
      matrix[N, M] x;
      vector[N] y;
    }
    parameters {
      vector[M] beta;
      real<lower=0> sigma;
    }
    model {
      y ~ normal(x * beta, sigma);
      for(i in 1:M)
        beta[i] ~ normal(0, 1000);
      sigma ~ chi_square(2);
    }
"

linreg <- stan(model_code = linmodel.stan, data = standata, 
    model_name = "Schweiz", iter = 2000, warmup = 1000)

Ergebnisse

Ein einfacher plot-Befehl zeigt uns Median, 80%- und 95%-Kredibilitätsintervalle der Parameter:

plot(linreg)

## ci_level: 0.8 (80% intervals)

## outer_level: 0.95 (95% intervals)

Mehr Details

Mit print erhalten wir mehr Details.

print(linreg)

## Inference for Stan model: swiss.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##            mean se_mean    sd    2.5%     25%     50%     75%   97.5% n_eff
## beta[1]   67.10    0.25 10.73   46.01   59.86   67.29   74.37   87.72  1782
## beta[2]   -0.17    0.00  0.07   -0.31   -0.22   -0.17   -0.13   -0.03  2031
## beta[3]   -0.27    0.01  0.26   -0.79   -0.44   -0.26   -0.10    0.23  2082
## beta[4]   -0.87    0.00  0.18   -1.21   -0.99   -0.87   -0.75   -0.52  2461
## beta[5]    0.10    0.00  0.03    0.03    0.08    0.10    0.13    0.17  2432
## beta[6]    1.08    0.01  0.38    0.34    0.83    1.08    1.33    1.83  2204
## sigma      7.07    0.01  0.72    5.84    6.56    7.01    7.48    8.65  2600
## lp__    -117.98    0.06  2.02 -123.08 -118.99 -117.58 -116.57 -115.19  1349
##         Rhat
## beta[1] 1.00
## beta[2] 1.00
## beta[3] 1.00
## beta[4] 1.00
## beta[5] 1.00
## beta[6] 1.00
## sigma   1.00
## lp__    1.01
## 
## Samples were drawn using NUTS(diag_e) at Fri May 10 14:02:46 2024.
## 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).

STAN kann intern mit mehreren parallelen Ketten rechnen. Sprich: der Algorithmus wird auf mehreren Rechenkernen gleichzeitig ausgeführt, die Ziehungen der verschiedenen Ketten werden dann zusammengeworfen.

ACF

Der ACF-Plot zeigt, dass die einzelnen Ziehungen im HMC-Algorithmus sehr schwach korreliert sind. Die Ketten sind wieder in unterschiedlichen Farben dargestellt.

coda::acfplot(beta)

Hier verwenden wir das coda-Paket, welches Funktionen zur Diagnostik von MCMC-Ziehungen bereitstellt. Genaueres dazu später.