Simulate from truncated multivariate normal distribution in Stan

Modeling
1

Hello Stan team,

I am looking for a way to simulate a truncated Multivariate normal distribution in STAN, i.e. something like multi_normal_rng(mu, Sigma), but with constraints on the bounds, in that each underlying univariate variable has its own bounded domain (or more specifically, some have 0 upper bound and -inf lower bound, and some have 0 lower bound and +inf upper bound). And of course Sigma is not diagonal :-)
I would need to do it generated_quantities() block.

Anyone having brilliant ideas on how to possibly do it ?

Thanks a lot,

Philippe

1 Like
2

You could do it “normally” in a model block (use mcmc to sample it). And just follow normal mcmc checks if posterior is valid.

Other way to do this is rejection sampling and in some cases it works.

1 Like
3

Hi Ari,

Thanks for your swift response. In fact I was looking for simulating without MCMC as I was planning to use it within a Laplace Approximation procedure. This procedure needs to use STAN in a optimizing() call to find the mode and the inverse Hessian matrix of the marginal posterior of some parameters, en subsequently simulate the truncated multivariate normal is this is the actual conditional distribution of the remaining parameters.

However, I realized soon after sending the current post, that simulating from a non independent multivariate truncated normal distribution, of respectable dimensions, is in itself a complex problem. Hence your answer makes sense :-) I will need to look for simplifying routes.
Thank you,
Philippe.

4

Did you see this How to truncate a multivariate normal distribution? (sidebar, how to access old stan mailing lists posts?)

The relevant code has the generated quantities truncated mvn

functions {
  vector[] make_stuff(vector mu, matrix L, vector b, vector s, vector u) {
    int K = rows(mu); vector[K] d; vector[K] z; vector[K] out[2];
    for (k in 1:K) {
      int km1 = k - 1;
      if (s[k] != 0) {
        real z_star = (b[k] - 
                      (mu[k] + ((k > 1) ? L[k,1:km1] * head(z, km1) : 0))) / 
                      L[k,k];
        real v; real u_star = Phi(z_star);
        if (s[k] == -1) {
          v = u_star * u[k];
          d[k] = u_star;
        }
        else {
          d[k] = 1 - u_star;
          v = u_star + d[k] * u[k];
        }
        z[k] = inv_Phi(v);
      }
      else {
        z[k] = inv_Phi(u[k]);
        d[k] = 1;
      }
    }
    out[1] = z;
    out[2] = d;
    return out;
  }
}
data {
  int<lower=2> K;             // number of dimensions
  vector[K] b;                // lower or upper bound
  
  // s[k] ==  0 implies no constraint; otherwise 
  // s[k] == -1 -> b[k] is an upper bound 
  // s[k] == +1 -> b[k] is a lower bound
  vector<lower=-1,upper=1>[K] s;
  
  vector[K] mu;
  cholesky_factor_cov[K,K] L;
}
parameters {
  vector<lower=0,upper=1>[K] u;
}
model {
  target += log(make_stuff(mu, L, b, s, u)[2]); // Jacobian adjustments
  // implicit: u ~ uniform(0,1)
}
generated quantities {
  vector[K] y = mu + L * make_stuff(mu, L, b, s, u)[1];
}
5

Hi spinkney,

No I did not. Thanks for forwarding this.

Cheers,

Philippe