Simple normal regression model with inequality constraints on the parameters

If your constraint is:

a \le D\beta \le w

Then if you multiply all by the inverse of D:

D^{-1}a \le D^{-1}D\beta \le D^{-1}w

This simplifies to:

D^{-1}a \le \beta \le D^{-1}w

Stan has overloaded the division operator / for matrices/vectors to perform this multiplication by an inverse in more numerically stable and efficient way. So you would express your model in Stan as:

data {
  int<lower=0> N; // Number of observations
  int<lower=0> K; // number of explanatory variables.
  matrix[N,K] X; // Matrix of explanatory variables 
  matrix[K,K] D; // Matrix that represents inequality constraints
  vector[N] y; // Vector of dependent variables
  vector[K] a; // Vector of lower inequality constraints
  vector[K] w; // Vector of upper inequality constraints
}
transformed data {
  vector[K] D_div_a = D / a;
  vector[K] D_div_w = D / w;
}
parameters {
  vector<lower=D_div_a, upper=D_div_w>[K] b; 
  real<lower=0> sigma;
}
model {
  // More efficient likelihood for linear regression
  y ~ normal_id_glm(X, 0, b, sigma);
}
4 Likes