You can probably do this but I’m not understanding how C, in the form you have, can be used. You have to make some assumptions to do this and I’m not sure those are met with C. Because Stan needs to know each constraint on the previous beta’s to work. The assumption is that for row i all the previous beta’s (j < i) are known and then we use those to construct the constraint on beta i.
For example, if the constraints are
c1 + c3 - 2 * c2 > 0
c3 + c5 - 2 * c4 > 0
In Stan we want
c1 unconstrained
c2 unconstrained
c3 > 2 * c2 - c1
c4 unconstrained
c5 > 2 * c4 - c3
However, the constraints on c3 and c5 will induce c1, c2, c4 to be non-uniform.
You can code this in the parameters block as
parameters {
real c1, c2;
real<lower=2 * c2 - c1> c3;
real c4;
real<lower=2 * c4 - c3> c5;
}
But that becomes unwieldy to do for many different combinations.
Let’s put the lower bounds in a matrix
M = | -1 2 0 0 0 |
| 0 0 -1 2 0 |
I_c = {3, 5} // says row[1] is for index 3, row[2] is for index 5
I_uc = {1, 2, 4}
parameters {
vector[K] x_raw;
}
transformed parameters {
real log_abs_jac = 0;
{
int R = rows(M);
int K = cols(M);
vector[K] x;
x[I_uc] = x_raw[I_uc];
for (r in 1:R) {
x[I_c[r]] = M[r] * x + exp(x_raw[I_c[r]]);
log_abs_jac += x_raw[I_c[r]]];
}
}
}
What does M look like if you just want the value to be constrained to 0 itself? Well, you just make a row of all 0s.
For e.g., let x[4] > 0. We now have
M = | -1 2 0 0 0 |
| 0 0 0 0 0 |
| 0 0 -1 2 0 |
I_c = {3, 4, 5}
I_uc = {1, 2}
I’m sure you know that having many constraints could lead you out of a feasible set.
I’d code this as an _lp function like:
functions {
vector positive_constrain_contrasts_lp(vector y, matrix M, array[] int I_c,
array[] int I_uc) {
int R = rows(M);
int K = cols(M);
vector[K] x;
x[I_uc] = y[I_uc];
for (r in 1 : R) {
x[I_c[r]] = M[r, 1 : I_c[r] - 1] * x[1 : I_c[r] - 1] + exp(y[I_c[r]]);
}
target += y[I_c];
return x;
}
}
data {
int p;
int cn;
matrix[cn, p] M;
array[cn] int I_c;
array[p - cn] int I_uc;
}
parameters {
vector[p] beta_raw;
}
transformed parameters {
vector[p] beta = positive_constrain_contrasts_lp(beta_raw, M, I_c, I_uc);
}
model {
beta ~ std_normal();
}
And you can test in R as
mod <- cmdstanr::cmdstan_model('test_c_fun.stan')
out1 <- mod$sample(
data = list(
p = 5,
cn = 2,
M = matrix(c(-1, 2, 0, 0, 0,
0, 0, -1, 2, 0), byrow = T, nr = 2),
I_c = c(3, 5),
I_uc = c(1, 2, 4)
),
parallel_chains = 4
)
out1$summary("beta")
beta <- out1$draws("beta", format = "matrix")
table(beta[, 1] + beta[, 3] - 2 * beta[, 2] > 0)
table(beta[, 3] + beta[, 5] - 2 * beta[, 4] > 0)
out2 <- mod$sample(
data = list(
p = 5,
cn = 3,
M = matrix(c(-1, 2, 0, 0, 0,
0, 0, 0, 0, 0,
0, 0, -1, 2, 0), byrow = T, nr = 3),
I_c = c(3, 4, 5),
I_uc = c(1, 2)
),
parallel_chains = 4
)
out2$summary("beta")
beta2 <- out2$draws("beta", format = "matrix")
table(beta2[, 1] + beta2[, 3] - 2 * beta2[, 2] > 0)
table(beta2[, 4] > 0)
table(beta2[, 3] + beta2[, 5] - 2 * beta2[, 4] > 0)