ordered_logistic.stan

data {
  int<lower=2> K; // number of categories
  int<lower=0> N; // number of observations
  int<lower=1> D; // number of predictors
  int<lower=1, upper=K> y[N];
  matrix[N, D] x;
}
parameters {
  vector[D] beta;
  ordered[K-1] c; // cutpoints
}
model {
  beta ~ normal(0, 1);
  c ~ normal(0, 5);
  y ~ ordered_logistic(x * beta, c);
}
generated quantities {
  vector[N] y_rep;
  for(n in 1:N) {
    y_rep[n] = ordered_logistic_rng(x[n, ] * beta, c);
  }
}