I’ve been struggling with getting an ordered probit model working in stan for the past couple of weeks, and am not sure where I’m going wrong. I am trying to reproduce some of the problems in Kruschke’s Doing Bayesian Data Analysis, trying to keep it as simple as possible - so the 1 group ordinal model from Chapter 23 in the book, similar to the R-Jags code here: https://github.com/davidkretch/dbda/blob/master/DBDA2Eprograms/Jags-Yord-Xnom1grp-Mnormal.R. I’ve looked for the Stan implementations of these models and I can’t find those either.

I’m working in python, and can get a simple model working in Jags pretty easily but my Stan translation gives me a mess - it runs but the samples make no sense. I have searched for complete Stan examples of ordinal probit response models or regression and only find partial solutions which haven’t worked when I’ve tried. I’m including the full python code here on the smallest self-contained problem I could find to help debug these problems, and tried to make the programs as similar to each other. Not sure of the best way to debug problems like this.

Thanks for any advice!

Here’s the pyjags code:

```
data={'y': np.array([3, 3, 3, 3, 1, 3, 2, 2, 4, 2]),
'N': 10,
'K': 5,
}
# the cutoff points
K=data['K']
c = np.arange(1.5, K)
c = np.ma.asarray(c)
c[1:-1] = np.ma.masked # these are the missing ones
data['c']=c
model_code= """
model {
for ( k in 2:(K-2) ) { # 1 and K-1 are fixed, not stochastic
c[k] ~ dnorm( k+0.5 , 0.25)
}
sigma ~ dunif( K/1000 , K*10 )
mu ~ dnorm( (1+K)/2 , 1/(K)^2 )
for ( i in 1:N ) {
pr[i,1] <- pnorm( c[1] , mu , 1/sigma^2 )
for ( k in 2:(K-1) ) {
pr[i,k] <- max( 0 , pnorm( c[ k ] , mu , 1/sigma^2 )
- pnorm( c[k-1] , mu , 1/sigma^2 ) )
}
pr[i,K] <- 1 - pnorm( c[K-1] , mu , 1/sigma^2 )
y[i] ~ dcat( pr[i,1:K] )
}
}
"""
import pyjags
model = pyjags.Model(model_code, data=data,
adapt=500,
chains=4)
samples = model.sample(500)
```

Here’s the pystan program:

```
import numpy as np
data={'y': np.array([3, 3, 3, 3, 1, 3, 2, 2, 4, 2]),
'N': 10,
'K': 5,
}
# the cutoff points
K=data['K']
# not needed in stan - done through init function
c = np.arange(1.5, K)
c = np.ma.asarray(c)
c[1:-1] = np.ma.masked # these are the missing ones
# print('c:\t\t{}'.format(c))
# data['c']=c
data
model_code= """
functions {
real mymax(real x,real y) {
if (x>y) {
return x;
} else {
return y;
}
}
}
data {
int<lower=2> K;
int<lower=0> N;
int<lower=1,upper=K> y[N];
}
parameters {
ordered[K-1] c;
real<lower=0> sigma;
real mu;
}
model {
vector[K] pr[N];
for ( k in 2:(K-2) ) { # 1 and K-1 are fixed, not stochastic
c[k] ~ normal( k+0.5 , 0.25);
}
sigma ~ uniform( K/1000 , K*10 );
mu ~ normal( (1+K)/2 , 1/(K)^2 );
for (i in 1:N) {
pr[i,1] = Phi( (c[1]-mu)/sigma );
for ( k in 2:(K-1) ) {
pr[i,k] = mymax( 0.0 , Phi( (c[k]-mu)/sigma ) - Phi( (c[k-1]-mu)/sigma ) );
}
pr[i,K] = 1 - Phi( (c[K-1]-mu)/sigma );
y[i] ~ categorical(pr[i]);
}
}
"""
def initf():
return {
'mu':mean(data['y']),
'sigma':2,
'c':arange(1.5,data['K']),
}
initf()
import pystan
model = pystan.StanModel(model_code=model_code)
fit=model.sampling(data=data,iter=500,init=initf)
samples=fit.extract()
```