How to add constraints to Stan model


I am attempting to add constraints to my Stan model, specifically for a simplex vector (w, r, m): simplex [m] theta_2 [M, K]; where m=3, M = 2, and K = 7. I aim to impose constraints such as w(2) < w(4) < w(6) and w(3) > w(5) > w(7). I intend to apply similar constraints to ‘r’ and ‘m’ as well. How can I implement this in the parameter block? The following is my Stan code for the parameter block.

// The input data is a vector 'y' of length 'N'.
data {
  int N; // number of observations (rally lengths)
  int M; // number of players
  int t1 [N]; // number of touches
  int t2 [N]; // number of touches
  int t3 [N]; // number of touches
  int K;
  int touch [N];
  matrix [N, 2] dat; // touches, indicator, servereID and receiverID 
  int ServerID [N];
  int ReceiverID [N];
  int <lower = 2> m; 
  row_vector [3] beta1;
  row_vector [3] beta2;


parameters {
     simplex [m] theta_1 [M];
     simplex [m] theta_2 [M, K];

// The model to be estimated.
model {
  // Define priors

 theta_1 ~ dirichlet(beta1);
  for(i in 1:M){
  for(j in 1:K){
    theta_2[i,j] ~ dirichlet(beta2);
   dat ~ tenismodel(theta_1, theta_2, ServerID, ReceiverID, t1, t2, t3, touch);

Hi, @niro. Our new syntax is

array[N] int t1;
array[M] simplex[m] theta_1;

though I’d urge you not to use both M and m as constants.

I didn’t understand specifically what you want to do as I couldn’t tell what w is supposed to be, but there’s no way to add constraints to a simplex in addition to its being a simplex. You can do the transform manually and add constraints and it shouldn’t be too much algebra/calculus to get it right.

Thanks, @Bob_Carpenter Bob for letting me know about the new syntax.

I want to know, how I can add constraint to my simplex [m] theta_2 [M, K];. W is one parameter in this theta_2[M,K}.. As an example, if M = 1 and K = 7 and one parameter inside this simplex is “W”, I want to make constraint for W such that, w1(2) < w1(4) < w1(6) and w1(3) > w1(5) > w1(7).

Can I do something like the bellow?.

parameters {
  simplex[m] theta_1[M];
  simplex[m] theta_2[M, K];

model {
  // Define priors
  theta_1 ~ dirichlet(beta1);
  theta_2 ~ dirichlet(beta2);

  // Add specific constraints for W in theta_2 

    for (i in 1:m) {
      // Constraints: w(2) < w(4) < w(6) and w(3) > w(5) > w(7)
      theta_2[1, 2, i] < theta_2[1, 4, i];
      theta_2[1, 4, i] < theta_2[1, 6, i];
      theta_2[1, 3, i] > theta_2[1, 5, i];
      theta_2[1, 5, i] > theta_2[1, 7, i];

  // Add likelihood and other model components here
  dat ~ tenismodel(theta_1, theta_2, ServerID, ReceiverID, t1, t2, t3, touch);

No, you can’t specify constraints in this way. To specify these relatively exotic constraints, you need to come up with a continuously differentiable one-to-one mapping between some parameters you declare and the set of simplexes that satisfy your constraint, and then implement this tranformation in either the transformed parameters or model blocks, and include an appropriate Jacobian adjustment as part of your prior specification.

For a simplex w such that w[2] < w[4] < w[6] and w[3] > w[5] > w[7], one solution would be to create an unordered 3-simplex S whose elements represent w[1], w[2] + w[4] + w[6], and w[3] + w[5] + w[7], and then two ordered 3-simplexes T and U, and then let w[1] = S[1], w[6,4,2] = S[2] * T, and w[3, 4, 5] = S[3] * U.

The Jacobians here will all be readily computable.

Edit: see here (particularly posts downthread) for details of how to construct ordered simplexes Ordered simplex constraint transform - #18 by martinmodrak

1 Like

Thanks @jsocolar. I will look into this more. It seems it is tricky to introduce such constraints to the model But I will give it a try. Thanks for the explanation.

It is tricky, which is why it’s not built in! You’re not the first person to ask about constraining simplexes. I believe the last person to ask wanted an ordered simplex, so that the probabilities are in ascending order. Here you may need a different simplex transform than the stick-breaking transform we use. For example, if you use something like a softmax transform, it’s easier to impose the constraints, but you have to be careful with identifiability.

Thanks, @Bob_Carpenter. I am still trying to solve this. Thanks again for your explanations.