# Why I cannot generate random numbers having a truncated lognormal distribution?

My deduction is:

When the distribution is truncated, a normalization factor should be introduced: $$$$g(x) = \frac{C}{x\sigma\sqrt{2\pi}}e^{-\frac{1}{2}\left(\frac{\ln{x}-\mu}{\sigma}\right)^2},$$$$ where $$C$$ is the factor. Integrating $$g(x)$$ over the interval $$[x_1, x_2]$$ ($$x_1 \geq 0$$), we have $$$$\frac{C}{\sigma\sqrt{2\pi}}\int_{x_1}^{x_2}\frac{1}{x}e^{-\frac{1}{2}\left(\frac{\ln{x}-\mu}{\sigma}\right)^2}\mathrm{d}x = 1. \label{CDF_truncated}$$$$ To determine $$C$$, let $$t = \frac{\ln{x}-\mu}{\sigma}$$, so $$x = e^{\sigma t + \mu}$$ and $$\mathrm{d}x = \sigma e^{\sigma t + \mu}\mathrm{d}t$$. The interval becomes $$\left[t_1, t_2\right]$$ which is given by $$$$t_1 = \frac{\ln{x_1}-\mu}{\sigma} \text{ and } t_2 = \frac{\ln{x_2}-\mu}{\sigma}. \notag$$$$ Hence, the inregration can be rewritten as $$$$\frac{C}{\sigma\sqrt{2\pi}}\int_{t_1}^{t_2} \sigma\cdot e^\mu \cdot e^{-\sigma t - \mu}\cdot e^{-\frac{t^2}{2} + \sigma t}\mathrm{d}t = 1, \notag$$$$ $$$$\Rightarrow \frac{C}{\sqrt{2\pi}}\int_{t_1}^{t_2} e^{-\frac{t^2}{2}}\mathrm{d}t = 1. \notag$$$$

Next, let $$m = \frac{t}{\sqrt{2}}$$ and $$t = \sqrt{2}m$$. So, $$\mathrm{d}t=\sqrt{2}\mathrm{d}m$$. The Eintegration can be further simplified into $$$$\frac{C}{\sqrt{2\pi}}\int_{m_1}^{m_2} \sqrt{2}e^{-m^2}\mathrm{d}m = 1, \notag$$$$ $$$$\Rightarrow \frac{C}{\sqrt{\pi}}\int_{m_1}^{m_2} e^{-m^2}\mathrm{d}m = 1. \notag$$$$. So, the $$C$$ is $$$$C = \frac{-2}{[{erf}(m_1)-{erf}(m_2)]},$$$$ where $${erf}(\cdot)$$ is the error function. If $$x_1 = 0 \text{ and } x_2 = +\infty$$, then $$t_1 = -\infty \text{ and } t_2 = +\infty$$, and $$m_1 = -\infty \text{ and } m_2 = +\infty$$, and $${erf}(m_1) \text{ and } {erf}(m_1)$$ are -1 and 1, respectively, $$C = 1$$ which leads to a regular non-truncated PDF.

Now, let us generate random log-normal variables from the standard uniform distributions. $$$$p = C\int_{x_1}^{x}\frac{1}{x\sigma\sqrt{2\pi}}e^{-\frac{1}{2}\left(\frac{\ln{x}-\mu}{\sigma}\right)^2}\mathrm{d}x.$$$$ where $$p$$ is a random number within [0, 1]. In a similar way, $$$$p = \frac{C}{ \sqrt{\pi}}\int_{m_1}^{m}e^{-m^2}\mathrm{d}m, \notag$$$$

$$$$\Rightarrow m = {erfinv}\left({{erf}(m_1) - \frac{2 p }{C}}\right), \notag$$$$ and $$m = \frac{\frac{\ln{x}-\mu}{\sigma}}{\sqrt{2}}$$, $$erfinv(\cdot)$$ is the inverse of the error function.

I write the following MATLAB script to generate random lognormal variables, with given expection and variance. But I got NaN.

clc;
close all;
clear all;

p = rand;

mean = 6; %given expection
variance = 0.5; % given variance
x1 = 1.5;
x2 = 15;

mu_1 = log(mean * mean / ((variance + mean * mean)^0.5)); %determine mean
sigma_1 = (log(1 + (variance) / (mean * mean)))^0.5; %determine std. deviation

t1 = (log(x1) - mu_1) / sigma_1;
t2 = (log(x2) - mu_1) / sigma_1;

m1 = t1 / (2^0.5);
m2 = t2 / (2^0.5);

C = -2 / (erf(m1) - erf(m2));

m = erfinv(erf(m1) - 2 * p / C);



m =

NaN

what's wrong with my deduction? or my code?

Update------------------------------------------

according to the comments of @jbowman, I write the following MATLAB code:

clc;
close all;
clear all;

mean = 6; %given expection
variance = 9; % given variance
x1 = 1.5;
x2 = 15;

mu_1 = log(mean * mean / ((variance + mean * mean)^0.5)); %determine mean
sigma_1 = (log(1 + (variance) / (mean * mean)))^0.5; %determine std. deviation

% then, convert normal to std. normal
% Any point (x) from a normal distribution can be converted to the standard normal distribution (z) with the formula z = (x-mean) / standard deviation
l = (log(x1) - mu_1) / sigma_1;
u = (log(x2) - mu_1) / sigma_1;

pd = makedist('Normal', 'mu', 0, 'sigma', 1);

p_l = cdf(pd, l);
p_u = cdf(pd, u);

array_rand_num = [];

nf = 1000;

for i = 1:1:nf
x = rand;

x_prime = p_l + (p_u - p_l) * x;

z = icdf(pd, x_prime);

p = exp(z * sigma_1 + mu_1);

array_rand_num(i) = p;
end

figure(1)
nbins = 30;
histfit(array_rand_num, nbins, 'lognormal');
xticks(0:1:20);
pd = fitdist(array_rand_num', 'lognormal');
fitting_Ex = exp(pd.mu + pd.sigma^2 * 0.5)
fitting_Dx = exp(2 * pd.mu + pd.sigma^2) * (exp(pd.sigma^2) - 1)
%since the pdf is truncated, the fitting Ex and Dx are not exactly equal to
%the input ones

• Why not just generate random variables from the untruncated lognormal distribution, then throw away the ones that are below / above the lower / upper bounds respectively? Sep 9, 2021 at 0:27
• @jbowman This is indeed an alternative way, also, could be my last resort. I just want a more scientific method. Sep 9, 2021 at 0:36
• Rather than throw anything away, you can generate a suitably truncated uniform (by simple rescaling to the corresponding truncation points), apply the inverse Gaussian to convert to truncated normal, and then exponentiate (or if you have the inverse cdf of a lognormal, you may be able to do the last two steps in one and it might save a little overhead) Sep 9, 2021 at 1:41
• stats.stackexchange.com/questions/56747/…
– whuber
Sep 9, 2021 at 12:39

It is a little easier notationally to work directly with the standard Normal distribution than with the lognormal distribution and error functions, and the conversion from the random number generated from a truncated standard Normal distribution to one generated from the desired truncated Lognormal distribution should be clear.

Let us label the lower and upper bounds (converted to the standard Normal distribution appropriately) as $$(l, u)$$ respectively. Then define $$p_l = \Phi(l)$$ and $$p_u = \Phi(u)$$, the values of the cumulative density function at $$l$$ and $$u$$. Generate random numbers $$z$$ as follows:

1. $$x \leftarrow U(0,1)$$
2. $$x' \leftarrow p_l + (p_u - p_l)x$$
3. $$z \leftarrow \Phi^{-1}(x')$$

The second assignment creates $$x' \sim U(p_l, p_u)$$, and running this through the inverse CDF function of a standard Normal variate will produce a random number that has a standard Normal distribution truncated at $$(l, u)$$.

Of course, if you have access to a function that will calculate the inverse CDF of a Lognormal distribution with specified parameters, you can work directly with that in step 3, saving yourself any effort of converting from the Lognormal to the standard Normal and back again.

• So, that means: I have to convert the $X \sim LN(\mu_1, \sigma_1)$ to $Y\sim N(\mu_2, \sigma_2)$ and then to $Z \sim N(0, 1)$. $l = (ln(x1) - \mu_2) / \sigma_2$. $\Phi$ is the CDF of $N(0, 1)$. Am I right? Sep 9, 2021 at 1:26
• The parameters of the lognormal are (typically) the same as those of the Normal. More simply, $z = (\ln(x)-\mu) / \sigma$, where $x \sim LN(\mu, \sigma)$ and $z \sim N(0,1)$. Sep 9, 2021 at 1:29