---
title: "Hamiltonian MC und STAN"
author: "Volker Schmid"
date: "12. Juni 2017"
output: beamer_presentation
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Hamiltonsche Mechanik

Stochastische (Radon-Nikodým-)Dichte $f$ entspricht physikalischer (Gibbs-)Energie $\phi$:

$$ p(x) \propto \exp(-\phi(x)) $$

**Hamiltonsche Mechanik**: 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) $$

## Hamiltonian MC (Hybrid MC)

* Benutze Hilfsvariable $v$ 
* Zusätzliche zufällige Komponente (ansonsten bleibt $E(x,v)$ konstant): Ziehe zufällig aus der Priori von $v$ (Normalverteilung)

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$ (einfach: Euler-Methode, besser: Leapfrog).

## Euler/Leapfrog

![20 Iterationen Euler vs. Leapfrog](hmc-leapfrog.png)

## 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 ($\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](https://arxiv.org/pdf/1206.1901.pdf)

## STAN (hier: RStan)  

```{r message=FALSE}
#install.packages('rstan')
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
```

Beispiel Lineare Regression
$$y_i \sim N(x_i'\beta,\sigma^2)$$
$$ \beta_j \sim N(0, 1000)$$
$$ \sigma^2 \sim \chi^2(2)$$

## STAN Modell 

    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);
    }
    
## STAN Beispiel (1)

```{r}
# Simulate some data
N <- 1000
M <- 3
x <- cbind(rep(1,N), rpois(N,5), runif(N,2,3))
truebeta <- c(-14, 0, 5)
truesigma <- 1
y <- rnorm(N, x%*%truebeta, truesigma)
```

## STAN Beispiel (2)

```{r cache=TRUE}
stan_data <- list(N=N, M=M, y=y, x=x)
limo <- stan(file = 'limo.stan', 
  model_name='LinearesModell', data=stan_data, 
  iter=1000, chains=4)
```

## STAN Beispiel Ergebnisse 
```{r}
print(limo)
```

## STAN Beispiel Ergebnisse 

```{r}
beta<-As.mcmc.list(limo,pars="beta")
beta0<-As.mcmc.list(limo,pars="beta[1]")
coda::traceplot(beta0)
```

## STAN Beispiel Ergebnisse 
```{r}
coda::densplot(beta0)
```

## STAN Beispiel Ergebnisse 
```{r}
  beta1<-As.mcmc.list(limo,pars="beta[2]")
  plot(beta1)
```

## STAN Beispiel Ergebnisse 
```{r}
  beta2<-As.mcmc.list(limo,pars="beta[3]")
  plot(beta2)
```

## STAN Beispiel Ergebnisse 
```{r}
coda::acfplot(beta)
```

## STAN Beispiel Ergebnisse 
```{r}
sigma<-As.mcmc.list(limo,pars="sigma")
plot(sigma)
```

## STAN Beispiel Ergebnisse 
```{r}
coda::acfplot(sigma)
```

## STAN Beispiel Ergebnisse 
```{r}
coda::gelman.diag(As.mcmc.list(limo))
```