Hello,

I have been trying for some time to fit a custom Gaussian dynamic linear system over a set of latent variables and the model just doesn’t converge even in the presented minimal example. I tried the QR decomposition, Cauchy reparametrization, changing the sampler settings (adapt delta, number of draws, thee max), what not…

My predictive checks look great, but the chain diagnostic is a complete disaster!

I am really out of options here so I hope someone will an idea how to improve the model or point me to a mistake that I am probably missing.

```
data {
int W; // Number of time points
int N; // Number of subjects
int Xdim;
int exo_q_num;
real U[N,W,exo_q_num];
int<lower=0> idx_itc_obs[N,W]; // Indices of time points with data
int<lower=1> P_itc; // number of parameters
int<lower=0> T_max_itc; // Max number of trials across subjects across time points
int<lower=0> Tr_itc[N, W]; // Number of trials for each subj for each time point
real<lower=0> delay_later[N, W, T_max_itc]; // Delay later - subj x time x trials
real<lower=0> amount_later[N, W, T_max_itc]; // Amount later - subj x time x trials
real<lower=0> amount_sooner[N, W, T_max_itc]; // Amount sooner - subj x time x trials
int<lower=-1, upper=1> choice_itc[N, W, T_max_itc]; // choice itc - 0 for instant reward, 1 for delayed reward - subj x time x trials
}
transformed data {
real Q = 1;
real cauchy_alpha = 10; // cauchy scale parameter hyperprior
int num_par = P_itc;
}
parameters {
real<lower=0> sigma_x; // for time point 1, sigma for mu_prior_x
real<lower=0> sigma_a; // sigma for A
real<lower=0> sigma_b; // sigma for B
real<lower=0> sigma_c; // sigma for C
real<lower=0> sigma_mu; // sigma for mu_d
real<lower=0> sigma_v; // for week 1, sigma for X
vector<lower=0>[num_par] sigma_r;
vector[num_par] mu_d; // constant offset
real mu_prior_x; // for week 1, X mean
real X[N,W,Xdim];
vector[Xdim] A; // dynamics matrix (here I create just the diagonal vector)
matrix[Xdim, exo_q_num] B;
matrix[num_par, Xdim] C;
real itc_k_pr[N,W];
real itc_beta_pr[N,W];
}
transformed parameters {
real<lower=0> itc_k[N,W];
real<lower=0> itc_beta[N,W];
for (n in 1:N) {
itc_k[n,] = exp(itc_k_pr[n,]);
itc_beta[n,] = exp(itc_beta_pr[n,]);
}
}
model {
// put sigma priors
sigma_v ~ cauchy(0,cauchy_alpha);
sigma_r ~ cauchy(0,cauchy_alpha);
sigma_x ~ cauchy(0,cauchy_alpha);
sigma_mu ~ cauchy(0,cauchy_alpha);
sigma_a ~ cauchy(0,cauchy_alpha);
sigma_b ~ cauchy(0,cauchy_alpha);
sigma_c ~ cauchy(0,cauchy_alpha);
mu_prior_x ~ normal(0,sigma_x); // prior on X mean
mu_d ~ normal(0,sigma_mu);
// put priors an A, B, C
A ~ normal(0,sigma_a);
to_vector(B) ~ normal(0,sigma_b);
to_vector(C) ~ normal(0,sigma_c);
{
for (s in 1:N) {
for (w in 1:W) {
vector[num_par] params_pr;
params_pr[1] = itc_k_pr[s,w];
params_pr[2] = itc_beta_pr[s,w];
if (w == 1) { // first time point
to_vector(X[s,w,]) ~ normal(mu_prior_x,sigma_v);
}
else {
to_vector(X[s,w,]) ~ normal(diag_matrix(A) * to_vector(X[s,w-1,]) + B * to_vector(U[s, w-1,]), Q);
}
params_pr ~ normal(mu_d + C * to_vector(X[s,w,]),sigma_r);
if (idx_itc_obs[s,w] != 0) { // if time point exists
vector[Tr_itc[s,w]] ev_later = to_vector(amount_later[s,w,:Tr_itc[s,w]]) ./ (1 + itc_k[s,w] * to_vector(delay_later[s,w,:Tr_itc[s,w]])/7);
vector[Tr_itc[s,w]] ev_sooner = to_vector(amount_sooner[s,w,:Tr_itc[s,w]]);
choice_itc[s,w,:Tr_itc[s,w]] ~ bernoulli_logit(itc_beta[s,w] * (ev_later - ev_sooner));
}
}
}
}
}
```

I use cmdstanpy with cmdstan 2.24.1

Any suggestion will be welcomed!