Sampling from a user-defined truncated distribution

Modeling
,
1

Hello Stan users,

I wanted to create a function which will sample from N(mu, sigma) T[0,1]
with some help I manage to create the function and it worked well. Here is the code:

functions {
real normal_lub_rng(real mu, real sigma, real lb, real ub) {
  real p_lb = normal_cdf(lb, mu, sigma);
  real p_ub = normal_cdf(ub, mu, sigma);
  real v = uniform_rng(p_lb, p_ub);
  real y = mu + sigma * inv_Phi(v);
  return y;
}
}

generated quantities {
real u[N];
real alpha_u;
real<lower=0> sigma_u;

alpha_u = normal_rng(0.20,0.10);
sigma_u = gamma_rng(0.15,1);
for (i in 1:N) 
u[i] = normal_lub_rng(alpha_u, sigma_u,0.01,1);

}

In the above function, The mean mu is real. But now I want it to be a vector so I came up with the code

functions {
real normal_lub_rng(real mu, real sigma, real lb, real ub) {
 real p_lb = normal_cdf(lb, mu, sigma);
real p_ub = normal_cdf(ub, mu, sigma);
  real v = uniform_rng(p_lb, p_ub);
real y = mu + sigma * inv_Phi(v);
  return y;
}
}




data {
int<lower=0> N; // Number of obs
}

generated quantities {
real u[N];
real alpha_u[N] ;
real<lower=0> sigma_u;


sigma_u = gamma_rng(0.15,1);
for (i in 1:N) {
alpha_u[i] = normal_rng(0.20,0.10);
u[i] = normal_lub_rng(alpha_u[i], sigma_u,0.01,1);
}
}

But I am getting an error here

> SYNTAX ERROR, MESSAGE(S) FROM PARSER:
Exception: Exception: uniform_rng: Upper bound parameter is 1, but must be greater than 1 (in ‘unknown file name’ at line 8)

I am not sure what I am doing wrong. Can some one please help.

Cheers

2

The problem is that the gamma(0.15, 1) distribution is highly peaked at zero and generates sigma_u values that are as low as 10^{-8}. Such low value of sigma_u means that p_lb and p_ub are zero within rounding error–in other words, there’s practically no chance of a draw from untruncated normal being in the interval. Such extreme truncation is implausible, and you probably need to rethink the model.

If you do need extreme truncation then you’ll have to special-case the tails of the distribution, possibly something along the lines of

real normal_lub_rng(real mu, real sigma, real lb, real ub) {
  print("mu: ", mu);
  print("sigma: ", sigma);
  if (lb < mu) {
    real p_lb = normal_cdf(lb, mu, sigma);
    real p_ub = normal_cdf(ub, mu, sigma);
    if (p_ub == 0) 
      return ub;
    real v = uniform_rng(p_lb, p_ub);
    real y = mu + sigma * inv_Phi(v);
    return y;
  } else {
    real p_lb = exp(normal_lccdf(lb | mu, sigma));
    real p_ub = exp(normal_lccdf(ub | mu, sigma));
    if (p_lb == 0) 
      return lb;
    real v = uniform_rng(p_ub, p_lb);
    real y = mu - sigma * inv_Phi(v);
    return y;
  }
}
3

Thanks very much for the help, I think I understand what you are saying.
Just one thing when you say :

…there’s practically no chance of a draw from untruncated normal being in the interval.

do you meant to say …

there’s practically no chance of a draw from truncated normal being in the interval

?

4

No, I meant what I said but I see it may be a bit confusing.
By definition, the draws from the truncated distribution are always inside the interval.
Stan doesn’t have a truncated normal rng so you must simulate the truncated distribution using the function normal_cdf() which is the cumulative probability of the untruncated distribution. This fails when the untruncated distribution has too low a probability in the truncation interval.

1 Like
5

Thaks for taking time to respond. It helps a great deal.