How to avoid the errors of rejecting initial value (an ordered variable)

Modeling
1

I want to set the varying ordered vector, the model like this:

data {
  int<lower=1> N;             // Number of observations
  int<lower=1> K;             // Number of ordinal categories
  int<lower=1> D;
  array[N] int<lower=1, upper=K> y; // Observed ordinals
  matrix[N,D] x;
  int g[N];
  int<lower=1> P;        //P different threshold 
}

parameters {
  vector[D] beta;       
 
  ordered[K-1] alpha;   // (Internal) cut points
  vector[K-1] yita[P-1];
  real<lower=0,upper=10> sigma;
}
transformed parameters{
  ordered[K - 1] thresh[P]; 
  thresh[1] = alpha;
  for (i in 2:P){
    thresh[i] = thresh[i-1]+ yita[i-1];
  }
}
model {
  beta~ normal(0,10);

  for (i in 1: (P-1))
    yita[i]~ normal(0,sigma);
    alpha ~  normal(0,10);

  for (i in 1:N)
     y[i]~ ordered_logistic(x[i]*beta,thresh[g[i]]);
}

I want to construct varying thresh. By adding a random variable yita. But how can I control the order? When sampling I have an error like this:

Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: model26c023f6ff1__namespace::log_prob: thresh[sym1__] is not a valid ordered vector. The element at 3 is 4.3765, but should be greater than the previous element, 4.45567 (in 'string', line 67, column 2 to column 27)
Chain 1: Rejecting initial value:

Sometimes the four chains can find the initial values, and I can get a convergent result. And sometimes I have only two chains that can get initial values. I think if there is a method to control I can get efficient sampling every time.

2

You need to declare the relevant constraints in the parameters block; the constraints in transformed parameters are only checked after-the-fact which is too late.
On the other hand, sampling statements can operate on any expression (because they don’t really “sample” but simply compute a log-probability density) (but nonlinear transformations need to be mindful of jacobians; in this case the difference is linear so it just works)

parameters {
  vector[D] beta;
  ordered[K - 1] thresh[P];
  real<lower=0,upper=10> sigma;
}
model {
  thresh[1] ~ normal(0, 10);
  for (i in 1: (P-1)) {
    (thresh[i+1] - thresh[i]) ~ normal(0,sigma);
  }
  ...
}
2 Likes
3

Very Thanks! it really works!