Sampling from truncated lognormal

Modeling
1

How do I generate samples from a truncated lognormal distribution in stan?

This gives me Exception: validate transformed params errors.

generated quantities{
  real<lower = 0, upper = 20> y_hat_level[NUMBER_OF_TRAINING_POINTS];

  for(training_point in 1:NUMBER_OF_TRAINING_POINTS){
    y_hat_level[training_point] = lognormal_rng(mu[training_point], sigma);
  }

thanks

1 Like
2

I think rejection sampling is the easiest to code, though things can get stuck in infinite loops if you’re unlucky:

generated quantities{
  real<lower = 0, upper = 20> y_hat_level[NUMBER_OF_TRAINING_POINTS];

  for(training_point in 1:NUMBER_OF_TRAINING_POINTS){
    real generated = lognormal_rng(mu[training_point], sigma);

    while(generated > 20) {
      generated = lognormal_rng(mu[training_point], sigma)
    }

    y_hat_level[training_point] = generated;
  }
}

The other option is computing an inverse CDF of the truncated distribution. Then if you generate a random number [0, 1] you can.

Yet another option is add y_hat_level as a parameter and add a statement to the model block to sample them with MCMC.

y_hat_level ~ lognormal(mu, sigma);

I suspect in this case you’d get mixing problems though. You could probably sample this on the log scale and then do an exp transformation and do a constrained non-centered parameterization there (Non-centered parameterisation with boundaries).

2 Likes
3

This also works (be sure to plot the output of the rng using draws and not mean). See the manual on truncated random number generation (18.10).

According to wikipedia the quantile function of the lognormal needs the inverse error function but Stan doesn’t have that. However, the math works out to just use inv_Phi() because

\begin{aligned} \Phi(x) &= \frac{1}{2\pi} \int_{-\infty}^x e^{-t^2}dt \\ &= \frac{1}{2} + \frac{1}{2}\text{erf}\bigg( \frac{x}{\sqrt{2}} \bigg) \end{aligned}

that implies that

\Phi^{-1}(x) = \sqrt{2} \, \text{inv_erf}(2x - 1) 
  real lognormal_quantile (real p, real mu, real sigma){
    return exp( mu + sigma * inv_Phi(p) );
  }
  
   real lognormal_trunc_rng(real mu, real sigma, real lb, real ub){
     real p_ub = lognormal_cdf(ub, mu, sigma);
     real p_lb = lognormal_cdf(lb, mu, sigma);
     real p = uniform_rng(p_lb, p_ub);
    
    return lognormal_quantile(p, mu, sigma);
  }

Where I’ve used

log_norm_loc <- function (mu, sigma) {
  log(mu^2 / sqrt(mu^2 + sigma^2))
}

log_norm_scale <- function (mu, sigma) {
  log( 1 + ( sigma^2 / mu^2 ) )
}
mu <- 15
sigma <- 10
> log_norm_loc(mu, sigma)
[1] 2.524188
> log_norm_scale(mu, sigma)
[1] 0.3677248

# lb = 0, ub = 20 in the stan program
mod_trunc <- mod$sample(
  data = list(N = 500,
              mu = log_norm_loc(mu, sigma),
              sigma = log_norm_scale(mu, sigma),
              p_uni = runif(500)),
  chains = 2,
  parallel_chains = 2
  
)
hist(mod_trunc$draws("y_stan"))

image
image967×792 9.98 KB

4 Likes