Add lower and upper bounds to ordered vector

[spinkney]

The proposal is to have the following

ordered<lower = a>[N] x1;
ordered<upper = b>[N] x2;
ordered<lower = a, upper = b>[N] x3;

For scalar a, b this is straightforward. When a, b are vectors they have to be ordered vectors themselves. In the case of having both a lb and ub then each element of a must be less than each element of b.

Here is Stan code for the scalar case

  vector ordered_lb_lp (vector y, real lb) {
    int N = rows(y);
    vector[N] x;
    
    x[1] = exp(y[1]) + lb;
    target += exp(y[1]);
    
    for (i in 2:N) {
      x[i] = exp(y[i]) + x[i - 1];
      target += exp(y[i]);
    }
    
    return x;
  }

  vector ordered_ub_lp (vector y, real ub) {
    int N = rows(y);
    vector[N] x;
    
    x[1] = ub - exp(y[1]);
    target += exp(y[1]);
    
    for (i in 2:N) {
      x[i] = x[i - 1] - exp(y[i]);
      target += exp(y[i]);
    }
    
    return x;
  }
  
  vector ordered_lb_ub_lp (vector y, real lb, real ub) {
    int N = rows(y);
    vector[N] x;
    
    x[1] = lb + (ub - lb) * inv_logit(y[1]);
    target += log(ub - lb) + log_inv_logit(y[1]) + log1m_inv_logit(y[1]);
    
    for (i in 2:N) {
      x[i] = x[i - 1] + (ub - x[i - 1]) * inv_logit(y[i]);
      target += log(ub - x[i - 1]) + log_inv_logit(y[i]) + log1m_inv_logit(y[i]);
    }
    
    return x;
    
  }

[nhuurre]

Yes, this is a good idea.

When a, b are vectors they have to be ordered vectors themselves.

That requirement is not mathematically necessary, failing it just makes the constraint partially redundant. Regardless, it is rather difficult to construct a smooth constraining function for vector bounds. Consider this:

data {
  real<lower=0, upper=0.99> x;
}
transformed data {
  ordered[2] a = [0,1-x]';
  ordered[2] b = [x,1]';
}
parameters {
  ordered<lower=a,upper=b>[2] y;
}

For x=0 it's just ordered<lower=0,upper=1>[2] y but for x > 0.5 the vector behaves like two independent parameters

parameters { // x=0.9
  real<lower=0, upper=0.1> y1;
  real<lower=0.9, upper=1> y2;

So I think the lower and upper bounds should always be scalars.

---

That example ordered_ub_lp is reverse-ordered. Ordered would be

  vector ordered_ub_lp (vector y, real ub) {
    int N = rows(y);
    vector[N] x;
    
    x[N] = ub - exp(y[N]);
    target += exp(y[N]);
    for (i in 1:N-1) {
      x[N-i] = x[N-i-1] - exp(y[N-i]);
      target += exp(y[N-i]);
    }
    
    return x;
  }

And I think ordered_lb_ub_lp should offset the unconstrained parameter so that zero-initialized value corresponds to an evenly spaced vector in the constrained space. (This makes it identical to cumulative sum of simplex constraint)

    for (i in 2:N) {
      x[i] = x[i - 1] + (ub - x[i - 1]) * inv_logit(y[i] - log(N+1-i));
      target += log(ub - x[i - 1]) + log_inv_logit(y[i] - log(N+1-i)) + log1m_inv_logit(y[i] - log(N+1-i));
    }

[WardBrian]

I think this would be no problem at the language level