Stan_Turing

Stan_Turing.jl

cd(@__DIR__)
using Pkg
Pkg.activate("")

using CmdStan, Random, Distributions, MCMCChains

model = "
data { 
    // total number of data points in y
    int<lower=1> N;
    // number of groups
    int<lower=1> G;
    // number of predictors
    int<lower=1> P;
    // group index 
    int<lower=1> group[N]; 
    vector[N] y;
    matrix[N,P] X;
} 

parameters {
  real<lower=0> sigma;
  real<lower=0> s;
  real mu_beta0;
  vector[G] beta0;
  vector[P] betas;
}

model {
    real sigma_beta0;
    sigma ~ lognormal(0,1);
    mu_beta0 ~ normal(0, 1);
    s ~ inv_gamma(3, 2);
    sigma_beta0 = sqrt(s);
    beta0 ~ normal(0, sigma_beta0);
    for(i in 1:N){
        y[i] ~ normal(mu_beta0 + beta0[group[i]] + X[i,:]*betas, sigma);
    }
}
"

proj_dir = @__DIR__
# add your path here if not set up on machine's path
set_cmdstan_home!("/home/dfish/cmdstan-2.19.1")
seed = 1051
# Generate data.
Random.seed!(seed)
# number of observations per group
N = 10
# number of groups
G = 100
# number of coefficients
P = 3
X = randn(N*G, P)
z = repeat(1:G, inner=N)
beta = [-2, 0, 2]
intercept = 1
random_intercepts = rand(Normal(0,1), G)
sigma = .1
y = intercept .+ X * beta + random_intercepts[z] + randn(N*G) * sigma

stan_data = Dict("N"=>N*G, "P"=>P, "G"=>G, "X"=>X, "group"=>z, "y"=>y)

config = Sample(; num_samples=1000, num_warmup=1000, adapt=CmdStan.Adapt(;delta=.65))
stanmodel = Stanmodel(config, name="model", model=model, printsummary=false,
    output_format=:mcmcchains, random = CmdStan.Random(seed));
    config.algorithm.engine.max_depth = 10
@time rc, samples, cnames = stan(stanmodel, stan_data, proj_dir);




using Turing, ReverseDiff, Zygote
using Distributions
import Random

Turing.setadbackend(:reversediff)
# Turing.setadbackend(:zygote)

@model multilevel_with_random_intercept(y, X, z) = begin
  # number of unique groups 
  num_random_intercepts = length(unique(z))

  # number of predictors
  num_predictors = size(X, 2) 

  ### NOTE: Carefully chosen priors should be used for a particular application. ###

  # Prior for standard deviation for errors.
  sigma ~ LogNormal()

  # Prior for coefficients for predictors.
  beta ~ filldist(Normal(), num_predictors)

  # Prior for intercept.
  intercept ~ Normal()

  # Prior for variance of random intercepts. Usually requires thoughtful specification.
  s2 ~ InverseGamma(3, 2)
  s = sqrt(s2)

  # Prior for random intercepts.
  random_intercepts ~ filldist(Normal(0, s), num_random_intercepts)
  
  # likelihood.
  y .~ Normal.(intercept .+ X * beta + random_intercepts[z], sigma)
end

Random.seed!(seed) 
N = 10
G = 200
P = 3
X = randn(N*G, P)
z = repeat(1:G, inner=N)
beta = [-2, 0, 2]
intercept = 1
random_intercepts = rand(Normal(0,1), G)
sigma = .1
y = intercept .+ X * beta + random_intercepts[z] + randn(N * G) * sigma

# Sample via NUTS.
@time chain = sample(multilevel_with_random_intercept(y, X, z), NUTS(1000,.65), 1000, progress=true)