In ChristofSeiler/curvature: Convergence Diagnosis of Hamiltonian Monte Carlo

This vignette shows how to calculate sectional curvature from an rstan object for the multivariate $t$ distribution with 50 degrees of freedom in 100 dimensions.

We need these packages.

library("rstan")
library("curvature")

Stan Model

Stan code for $t$ distribution.

model_code <- "
  data {
    int<lower=1> nu;
    int<lower=1> K;
    vector[K] mu;
    cov_matrix[K] Sigma;
  } 
  parameters {
    vector[K] y;
  } 
  model {
    y ~ multi_student_t(nu, mu, Sigma);
  } 
"
model = stan_model(model_code = model_code)

Simulation

Prepare parameters for simulations.

d = 100 # dimension of distribution
mu = rep(0,d) # mean
Sigma = diag(1,d) # covariance
nu = 50 # degrees of freedom
T0 = 1000 # burn-in
T = 100 # after burn-in
data = list(nu = nu,
            mu = mu, 
            Sigma = Sigma, 
            K = d)
fit = sampling(model, 
               data = data,
               chains = 1, 
               iter = T+T0, 
               warmup = T0,
               algorithm = "HMC",
               control = list(metric = "unit_e"))
res = curvature(fit = fit)
plot(res)
summary(res)

Session Info

sessionInfo()

ChristofSeiler/curvature documentation built on May 28, 2019, 12:17 p.m.