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

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

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"))
```

4 Likes