Constraints during sampling

General
, , , ,
1

Is there a way to incorporate, say, hyperplane constraints while sampling in stan? For example, suppose I want to obtain posterior samples of b in 5 dimensions such that Ab \geq c, where A is a 405 matrix and c is a 401 vector of reals and I put in some standard prior on b (suppose multivariate Gaussian).

2

As of Stan 2.24 you can use upper/lower bounds that are that same type and size as the given variable.

So if A is 40x5 matrix and b is 5x1 vector, then you write your parameters as:

parameters {
  matrix[40,5] A;
  vector[5] b;
  vector<upper=(A*b)>[40] c; 
}

But you’ll need to use cmdstan or one of the cmdstanR/cmdstanPy interfaces for this

3

Just noticed that you wrote Ab \geq c, whereas the constraint above implies Ab \gt c. You’ll probably want to write the constraint as:

vector<upper=(A*b + machine_precision())>[40] c;

So that Ab == c is admissible

4

Thanks a lot for your reply! However, in my example, neither A nor c are parameters - only b is a parameter. Say a special case could be b \geq 0, where c = 0 vector and A is the identity matrix. Are these constraints possible to implement in STAN?

5

So you have five parameters and 40 constraints. Are some of these constraints redundant or is the sampling region a 40-facet polytope?

Btw, it’s possible that Andrew’s solution still works. This is a bit hackish but I think in theory it’s correct.

data {
  matrix[40,5] A;
  vector[40] c; 
}
parameters {
  vector[5] b;
  vector<upper=(A*b)>[40] c_x; 
}
model {
  c ~ normal(c_x, 0.001); // small variance to make c_x and c equal
  ...
}

In practice the sampler may get stuck if the variance is too small.

2 Likes
6

Bit late coming back to this sorry, but if you’re working with a case where A is invertible, then you can take advantage of the fact that:

Ab \ge c

Is equivalent to:

b \ge A^{-1}c

(I think, been a while since I’ve dealt with matrix algebra properly!)

Then your parameter declaration for b would look like:

parameters {
  vector<lower=(inverse(A)*c - machine_precision())>[R] b;
}

But invertibility isn’t always possible, so @nhuurre’s suggestion is a good alternative