ggmcmc.R

shige — Dec 28, 2013, 9:11 PM

library(ggmcmc)

Loading required package: plyr
Loading required package: reshape2
Loading required package: ggplot2

library(coda)

Loading required package: lattice

N <- 500
X <- runif(N)
y <- 3 + 5 * X + rnorm(N)

D <- list(N=N, X=X, y=y)

md <- '
data {
  int<lower=0> N;
  real X[N];
  real y[N];
}
parameters {
  real theta[2];
  real<lower=0> sigma;
}
model {
  for (n in 1:N) {
    y[n] ~ normal(theta[1] + theta[2] * X[n], sigma);
  }
}
'

library(rstan)

Loading required package: Rcpp
Loading required package: inline

Attaching package: 'inline'

The following object is masked from 'package:Rcpp':

    registerPlugin

rstan (Version 2.1.0, packaged: 2013-12-27 18:21:10 UTC, GitRev: 548aa7bbbb89)

Attaching package: 'rstan'

The following object is masked from 'package:coda':

    traceplot

fit <- stan(model_code=md, data=D)

TRANSLATING MODEL 'md' FROM Stan CODE TO C++ CODE NOW.
COMPILING THE C++ CODE FOR MODEL 'md' NOW.
SAMPLING FOR MODEL 'md' NOW (CHAIN 1).

Iteration:    1 / 2000 [  0%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration:  600 / 2000 [ 30%]  (Warmup)
Iteration:  800 / 2000 [ 40%]  (Warmup)
Iteration: 1000 / 2000 [ 50%]  (Warmup)
Iteration: 1200 / 2000 [ 60%]  (Sampling)
Iteration: 1400 / 2000 [ 70%]  (Sampling)
Iteration: 1600 / 2000 [ 80%]  (Sampling)
Iteration: 1800 / 2000 [ 90%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)
Elapsed Time: 0.4 seconds (Warm-up)
              0.49 seconds (Sampling)
              0.89 seconds (Total)

SAMPLING FOR MODEL 'md' NOW (CHAIN 2).

Iteration:    1 / 2000 [  0%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration:  600 / 2000 [ 30%]  (Warmup)
Iteration:  800 / 2000 [ 40%]  (Warmup)
Iteration: 1000 / 2000 [ 50%]  (Warmup)
Iteration: 1200 / 2000 [ 60%]  (Sampling)
Iteration: 1400 / 2000 [ 70%]  (Sampling)
Iteration: 1600 / 2000 [ 80%]  (Sampling)
Iteration: 1800 / 2000 [ 90%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)
Elapsed Time: 0.39 seconds (Warm-up)
              0.44 seconds (Sampling)
              0.83 seconds (Total)

SAMPLING FOR MODEL 'md' NOW (CHAIN 3).

Iteration:    1 / 2000 [  0%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration:  600 / 2000 [ 30%]  (Warmup)
Iteration:  800 / 2000 [ 40%]  (Warmup)
Iteration: 1000 / 2000 [ 50%]  (Warmup)
Iteration: 1200 / 2000 [ 60%]  (Sampling)
Iteration: 1400 / 2000 [ 70%]  (Sampling)
Iteration: 1600 / 2000 [ 80%]  (Sampling)
Iteration: 1800 / 2000 [ 90%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)
Elapsed Time: 0.4 seconds (Warm-up)
              0.4 seconds (Sampling)
              0.8 seconds (Total)

SAMPLING FOR MODEL 'md' NOW (CHAIN 4).

Iteration:    1 / 2000 [  0%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration:  600 / 2000 [ 30%]  (Warmup)
Iteration:  800 / 2000 [ 40%]  (Warmup)
Iteration: 1000 / 2000 [ 50%]  (Warmup)
Iteration: 1200 / 2000 [ 60%]  (Sampling)
Iteration: 1400 / 2000 [ 70%]  (Sampling)
Iteration: 1600 / 2000 [ 80%]  (Sampling)
Iteration: 1800 / 2000 [ 90%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)
Elapsed Time: 0.42 seconds (Warm-up)
              0.52 seconds (Sampling)
              0.94 seconds (Total)

s <- mcmc.list(lapply(1:ncol(fit), function(x) mcmc(as.array(fit)[,x,])))

S <- ggs(s)

ggs_histogram(S)

ggs_traceplot(S)

ggs_density(S)

ggs_compare_partial(S)

ggs_running(S)

ggs_autocorrelation(S)

ggs_crosscorrelation(S)

ggs_caterpillar(S)