I’ve been playing with the R package limSolve’s xsample function recently. For those not aware, it provides a MCMC approach to simulating from a distribution with equality and inequality constraints.

I was able to implement the same approach in Stan for equality constraints (most of the hard work is done outside of Stan). For some variable x, it expresses the problem in terms of x = x0 + Z * q, where x0 is the result of a least squares, equality-constrained problem and Z is the matrix that provides a basis for the null space (using MASS).

So far I haven’t been able to properly implement the inequality constraints in the general case (I was able to get something that produced what I expected where there are N independent inequality constraints, but with lots of divergences). I would be interested if anyone can get that to work.

The R code below creates constraints E and F that that restrict the sum of the variable x to be 1 and the sum of the first and second x’s to be 0.3.

```
library(rstan)
rstan_options(auto_write = TRUE)
library(limSolve)
library(MASS)
wd <- "C:/..." #location of where file is saved
setwd(wd)
N <- 10
E <- matrix(rep(1, N), 1, N)
E <- rbind(E, matrix(c(c(1, 1), rep(0, N - 2)), 1, N))
F <- matrix(1, 1, 1)
F <- rbind(F, matrix(0.3, 1, 1))
x0pre <- lsei(E = E, F = F)
x0 <- matrix(x0pre$X, N, 1)
tol <- sqrt(Machine$double.eps)
Z <- Null(t(E));
Z[abs(Z)<tol] <- 0
stan_data <- list(N = N, E = E, N_Z = dim(Z)[2], Z = Z, x0 = as.numeric(x0))
fit_eq <- stan(file = 'equality_constr.stan', data = stan_data,
verbose = FALSE, iter = 500, chains = 4)
e_eq_x <- extract(fit_eq, pars = "x", include = TRUE, permuted=TRUE)
e_eq_x <- e_eq_x$x
rowSums(e_eq_x)
rowSums(e_eq_x[, 1:2])
```

This is the accompanying stan code.

```
data {
int<lower=0> N;
int<lower=0, upper=N> N_Z;
matrix[N, N_Z] Z;
vector[N] x0;
}
parameters {
vector<lower=-1, upper=1>[N_Z] q;
}
transformed parameters {
vector[N] x;
x = x0 + Z * q;
}
model {
}
****
```