In 2013, I requested a multivariate normal integral feature, but it looks like this isn’t going to happen, at least not any time soon.

So, I was thinking it might make sense to try to get the same result by modeling a latent multivariate normal variable and summing it over different regions. I wrote a Stan model today, and I got it to compile after some debugging, but when I try to fit it, it seems to hang (or maybe it’s just taking a very long time to do anything). It’s been a while since I wrote a Stan model from scratch, so I’m sure there are lots of things about this that are suboptimal (for the past few years, I have been using PyMC 2 to fit the models that led me to request the multivariate normal integral feature - my PyMC models are probably also suboptimal, but they work, so…).

Any advice would be greatly appreciated. I would love to get this stuff working in Stan if possible. Ultimately, I want to scale this up to fit a fairly large number of matrices simultaneously (for fairly complicated multilevel structure in my data). Here’s the model I wrote:

```
data {
int<lower=0> n_stim; // number of stimulus categories
int<lower=0> n_resp; // number of response categories
int<lower=0> max_N; // number of simulated latent perceptual effects
int<lower=0> cm[n_stim, n_resp]; // confusion matrix
real<lower=-1,upper=1> signs[n_resp, 2]; // assumes 2D
}
parameters {
vector[2] mu[n_stim]; // perceptual distribution means
matrix[2,2] L[n_stim]; // Cholesky factored perceptual correlation matrices
}
transformed parameters {
vector[2] mu_p[n_stim]; // prior means for mu
matrix[2,2] mu_S; // prior variance for mu
real<lower=0> tm[n_stim, n_resp]; // tallied latent variables
vector[n_resp] pm[n_stim]; // predicted response probabilities
vector[2] pe[n_stim,max_N]; // latent perceptual effects
mu_S[1,1] = 2;
mu_S[1,2] = 0;
mu_S[2,1] = 0;
mu_S[2,2] = 2;
for (i in 1:n_stim) {
mu_p[i,1] = -1*signs[i,1];
mu_p[i,2] = -1*signs[i,2];
for (j in 1:n_resp) {
for (k in 1:max_N) {
tm[i,j] += (signs[j,1]*pe[i,k,1] < 0) && (signs[j,2]*pe[i,k,2] < 0);
}
pm[i,j] = tm[i,j]/sum(tm[i]);
}
}
}
model {
for (i in 1:n_stim) {
pe[i] ~ multi_normal_cholesky(mu[i], L[i]);
cm[i] ~ multinomial(pm[i]);
L[i] ~ lkj_corr_cholesky(2.0);
mu[i] ~ multi_normal(mu_p, mu_S);
}
}
generated quantities {
matrix[2,2] S[n_stim];
for (i in 1:n_stim) {
S[i] = multiply_lower_tri_self_transpose(L[i]);
}
}
```

And here is Python code to generate data and fit it with pystan:

```
import numpy as np
import pystan as ps
import pickle
from glob import glob
signs = np.array([[1,1],[1,-1],[-1,1],[-1,-1]])
mu = np.array([[-1,-1],[-1,1],[1,-1],[1,1]])
R = np.array([[[1,.5],[.5,1]],
[[1,-.5],[-.5,1]],
[[1,-.5],[-.5,1]],
[[1,.5],[.5,1]]])
N = 100
cm = np.zeros((4,4), dtype=int)
for i in range(4):
y = np.random.multivariate_normal(mu[i], R[i], size=N)
for j in range(4):
cm[i,j] = np.sum( (signs[j,0]*y[:,0] < 0) & (signs[j,1]*y[:,1] < 0) )
data_dict = {'n_stim':4, 'n_resp':4, 'max_N':N, 'cm':cm, 'signs':signs}
pkl_files = glob('*pkl')
model_pkl = 'simple_grt_mod.pkl'
if model_pkl not in pkl_files:
model = ps.StanModel(file='simple_grt.stan')
with open(model_pkl,'wb') as f:
pickle.dump(model,f)
else:
with open(model_pkl,'rb') as f:
model = pickle.load(f)
fit_model = True
if fit_model:
n_iter = 5000
warmup = 2500
thin = 1
fit = model.sampling(data=data_dict, iter=n_iter, warmup=warmup, thin=thin, n_jobs=2)
fit_pkl = 'simple_grt_fit.pkl'
smp_pkl = 'simple_grt_smp.pkl'
with open(fit_pkl,'wb') as f:
pickle.dump(fit,f)
with open(smp_pkl,'wb') as f:
pickle.dump(fit.extract(),f)
```