I am trying to reproduce this model from PyMC into Stan (from this publication):
with pm.Model(coords=coords) as new_hbm:
if sampler == 'MCMC':
x_input = pm.Data('x_input1', X_data)
y_input = pm.Data('y_input', y_data)
group_idx = pm.Data(f'{hgroup}_idx', idx, dims='Obs')
elif sampler == 'JAX':
x_input = X_data.values
y_input = y_data
group_idx = idx
sd_dist = pm.HalfNormal.dist(1.0)
# get back standard deviations and rho:
chol, corr, stds = pm.LKJCholeskyCov("chol", n=X_data.shape[1]+1, eta=4.0, sd_dist=sd_dist, compute_corr=True)
# Model global priors (Hyperpriors)
alpha = pm.Normal('globl_alpha', mu=0, sigma=1.0)
mu_beta1 = pm.Normal('global_beta', mu=0, sigma=1.0, dims='Var0')
# Model group level priors
mu_z = pm.Normal("mu_z", mu=0, sigma=1.0, dims=(hgroup,'Var1'))
alpha_beta = pm.Deterministic(f'group_beta_{hgroup}', pm.math.dot(chol, mu_z.T).T, dims=(hgroup, 'Var1'))
#Model
mean = alpha + alpha_beta[group_idx,0] + (mu_beta1 + alpha_beta[group_idx,1:]*x_input).sum(axis=1)
#Model Error
sigma = pm.HalfNormal('sigma', 1.0)
# Define likelihood
y_pred = pm.Normal('y_pred', mu=mean, sigma=sigma, observed=y_input, testval=1, dims="Obs")
This is my Stan code
// PARTIAL POOLING
data{
int<lower=0> N; // number of training instances (samples)
int<lower=0> K; // number of input variables (features)
int<lower=0> J; // number of landcover classes (group-level)
int<lower=0, upper=J> landcover_idx[N]; // indicator variable for each landcover class (u)
matrix[N, K] X; // input matrix (X)
vector[N] y; // target vector (Y)
}
parameters{
// vector[J] alpha; // LC-specific intercepts
vector[K] beta; // LC-specific regression coefficients
real<lower=0> sigma;
// modelling covariance structure of random effects
matrix[K, J] z;
cholesky_factor_corr[K] L_Omega;
vector<lower=0,upper=pi()/2>[K] tau_unif;
}
transformed parameters {
matrix[J, K] u; // random effects
vector<lower=0>[K] tau; // prior scale
//
for (k in 1:K){
tau[k] = 2.5 * tan(tau_unif[k]);
};
u = (diag_pre_multiply(tau, L_Omega) * z)';
}
model {
beta ~ normal(0, 2.5);
sigma ~ normal(0, 1);
to_vector(z) ~ std_normal();
L_Omega ~ lkj_corr_cholesky(2);
vector[N] y_hat;
y_hat = rows_dot_product(rep_matrix(beta, N)' + u[landcover_idx], X);
y ~ normal(y_hat, sigma);
}
The python code for creating the data and sampling from it using cmdstanpy
import pandas as pd
import numpy as np
from cmdstanpy import cmdstan_path, CmdStanModel, CmdStanMCMC
def make_dummy_data():
variables = [f"SM_lag_{i}" for i in range(6)]
variables += [f"PCP_lag_{i}" for i in range(6)]
variables += [f"VCI_lag_{i}" for i in range(6)]
N_samples = 100
X_train = pd.DataFrame({
var: np.random.random(N_samples)
for var in variables
})
y_train = np.random.random(len(X_train))
lc_idx = np.repeat([0, 1, 2], 20)[:len(X_train)]
K = len(X_train.columns)
data = dict(
N=X_data.shape[0],
K=K,
X_train=X_data.values,
y_train=y_data.values,
landcover_idx=lc_idx,
J=len(np.unique(lc_idx))
)
return data
if __name__ == "__main__":
data = make_dummy_data()
# build model
stan_file = "model_hierarchical.stan"
stan_model = CmdStanModel(stan_file=stan_file)
stan_model.compile()
# fit model
model_fit: CmdStanMCMC = stan_model.sample(
data=data,
chains=4,
parallel_chains=4,
seed=1111,
show_progress=True,
)
The model compiles but I get the following error:
In [147]: %run stan_test.py
INFO:cmdstanpy:found newer exe file, not recompiling
INFO:cmdstanpy:found newer exe file, not recompiling
INFO:cmdstanpy:CmdStan start procesing
chain 1 | | 00:00 StatuERROR:cmdstanpy:Chain [3] error: error during processing Operation not permitted | 00:00 Status
ERROR:cmdstanpy:Chain [1] error: error during processing Operation not permitted | 00:00 Status
ERROR:cmdstanpy:Chain [2] error: error during processing Operation not permitted | 00:00 Status
ERROR:cmdstanpy:Chain [4] error: error during processing Operation not permitted
chain 1 |███████████████████████████████████████████████████████████████████████████████████████████████| 00:00 Sampling completed
chain 2 |███████████████████████████████████████████████████████████████████████████████████████████████| 00:00 Sampling completed
chain 3 |███████████████████████████████████████████████████████████████████████████████████████████████| 00:00 Sampling completed
chain 4 |███████████████████████████████████████████████████████████████████████████████████████████████| 00:00 Sampling completed
INFO:cmdstanpy:CmdStan done processing.
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
~/github/LEARNING/bayes_hierarchical_model/stan_test.py in <module>
90
91 # fit model
---> 92 model_fit: CmdStanMCMC = stan_model.sample(
93 data=data,
94 chains=4,
~/miniconda3/envs/cmdstanpy/lib/python3.9/site-packages/cmdstanpy/model.py in sample(self, data, chains, parallel_chains, threads_per_chain, seed, chain_ids, inits, iter_warmup, iter_sampling, save_warmup, thin, max_treedepth, metric, step_size, adapt_engaged, adapt_delta, adapt_init_phase, adapt_metric_window, adapt_step_size, fixed_param, output_dir, sig_figs, save_latent_dynamics, save_profile, show_progress, show_console, refresh, time_fmt, force_one_process_per_chain)
1083 msg, runset.__repr__()
1084 )
-> 1085 raise RuntimeError(msg)
1086
1087 mcmc = CmdStanMCMC(runset)
RuntimeError: Error during sampling:
Exception: matrix[multi] row indexing: accessing element out of range. index 0 out of range; expecting index to be between 1 and 1 (in '/Users/tommylees/github/LEARNING/bayes_hierarchical_model/model_hierarchical.stan', line 45, column 2 to column 71)
Exception: matrix[multi] row indexing: accessing element out of range. index 0 out of range; expecting index to be between 1 and 1 (in '/Users/tommylees/github/LEARNING/bayes_hierarchical_model/model_hierarchical.stan', line 45, column 2 to column 71)
Exception: matrix[multi] row indexing: accessing element out of range. index 0 out of range; expecting index to be between 1 and 1 (in '/Users/tommylees/github/LEARNING/bayes_hierarchical_model/model_hierarchical.stan', line 45, column 2 to column 71)
Exception: matrix[multi] row indexing: accessing element out of range. index 0 out of range; expecting index to be between 1 and 1 (in '/Users/tommylees/github/LEARNING/bayes_hierarchical_model/model_hierarchical.stan', line 45, column 2 to column 71)
Exception: matrix[multi] row indexing: accessing element out of range. index 0 out of range; expecting index to be between 1 and 1 (in '/Users/tommylees/github/LEARNING/bayes_hierarchical_model/model_hierarchical.stan', line 45, column 2 to column 71)
Exception: matrix[multi] row indexing: accessing element out of range. index 0 out of range; expecting index to be between 1 and 1 (in '/Users/tommylees/github/LEARNING/bayes_hierarchical_model/model_hierarchical.stan', line 45, column 2 to column 71)
Exception: matrix[multi] row indexing: accessing element out of range. index 0 out of range; expecting index to be between 1 and 1 (in '/Users/tommylees/github/LEARNING/bayes_hierarchical_model/model_hierarchical.stan', line 45, column 2 to column 71)
Exception: matrix[multi] row indexing: accessing element out of range. index 0 out of range; expecting index to be between 1 and 1 (in '/Users/tommylees/github/LEARNING/bayes_hierarchical_model/model_hierarchical.stan', line 45, column 2 to column 71)Command and output files:
RunSet: chains=4, chain_ids=[1, 2, 3, 4], num_processes=4
cmd (chain 1):
['/Users/tommylees/github/LEARNING/bayes_hierarchical_model/model_hierarchical', 'id=1', 'random', 'seed=1111', 'data', 'file=/var/folders/q3/0lmt64ld10s14_0n0vxpt_m00000gp/T/tmpk68_iqww/z8xlwpnn.json', 'output', 'file=/var/folders/q3/0lmt64ld10s14_0n0vxpt_m00000gp/T/tmpk68_iqww/model_hierarchical-20211222000529_1.csv', 'method=sample', 'algorithm=hmc', 'adapt', 'engaged=1']
retcodes=[1, 1, 1, 1]
per-chain output files (showing chain 1 only):
csv_file:
/var/folders/q3/0lmt64ld10s14_0n0vxpt_m00000gp/T/tmpk68_iqww/model_hierarchical-20211222000529_1.csv
console_msgs (if any):
/var/folders/q3/0lmt64ld10s14_0n0vxpt_m00000gp/T/tmpk68_iqww/model_hierarchical-20211222000529_0-stdout.txt
I’m not entirely sure what’s going on here, I am using the model defined here: Custom hierarchical logistic reg has larger std errors than stan_glmer results
as a reference for defining the cholseksy factorisation.
Any help would be super appreciated! Happy Christmas!