What does truncated distribution mean?

Question

In a research article about sensitivity analysis of an ordinary differential equation model of a dynamic system, the author provided the distribution of a model parameter as Normal distribution (mean=1e-4, std=3e-5) truncated to the range [0.5e-4 1.5e-4]. He then uses samples from this truncated distribution for simulations of the model. What does it mean to have a truncated distribution and sample from this truncated distribution?

I could come up with two ways of doing this:

• Sample from a Normal distribution but ignore all the random values falling outside the specified range before simulations.
• Somehow obtain a special "Truncated Normal" distribution and get samples from it.

Are these valid and equivalent approaches?

I believe in the first case, if one were to plot the experimental cdf/pdf of the sample, it would not look like a Normal distribution because the curves do not extend to $\pm\infty$.

Answer 1

To truncate a distribution is to restrict its values to an interval and re-normalize the density so that the integral over that range is 1.

So, to truncate the $N(\mu, \sigma^{2})$ distribution to an interval $(a,b)$ would be to generate a random variable that has density

$$ p_{a,b}(x) = \frac{ \phi_{\mu, \sigma^{2}}(x) }{ \int_{a}^{b} \phi_{\mu, \sigma^{2}}(y) dy } \cdot \mathcal{I} \{ x \in (a,b) \} $$

where $\phi_{\mu, \sigma^{2}}(x)$ is the $N(\mu, \sigma^2)$ density. You could sample from this density in a number of ways. One way (the simplest way I can think of) to do this would be to generate $N(\mu, \sigma^2)$ values and throw out the ones that fall outside of the $(a,b)$ interval, as you mentioned. So, yes, those two bullets you listed would accomplish the same goal. Also, you are right that the empirical density (or histogram) of variables from this distribution would not extend to $\pm \infty$. It would be restricted to $(a,b)$, of course.

Answer 2

Simulating from the normal $\mathcal{N}(\mu,\sigma^2)$ distribution until the outcome falls within an interval $(a,b)$ is fine when the probability $$ \varrho = \int_a^b \varphi_{\mu,\sigma^2}(x)\,\text{d} x $$ is large enough. If it is too small, this procedure is too costly since the average number of draws for one acceptance is $1/\varrho$.

As described in Monte Carlo Statistical Methods (Chapter 2, Example 2.2), as well as in my arXiv paper, a more efficient way to simulate this truncated normal is to use an accept-reject method based on an exponential $\mathcal{E}(\alpha)$ distribution.

Consider, without loss of generality, the case $\mu = 0$ and $\sigma = 1$. When $b=+\infty$, a potential instrumental distribution is the translated exponential distribution, $\mathcal{E} (\alpha,{ a})$, with density $$ g_{\alpha}(z) = \alpha e^{- \alpha(z - {a})} \; \mathbb{I}_{z \geq {a }} \;. $$ The ratio $$ p_{a,\infty}(z)/g_{\alpha}(z) \propto e^{- \alpha(z - a )}e^{-z^{2}/2} $$ is then bounded by $\exp(\alpha^{2}/2 - \alpha{a })$ if $\alpha > a$ and by $\exp(- a^{2}/2)$ otherwise. The corresponding (upper) bound is $$ \begin{cases} 1/\alpha \; \exp (\alpha^{2}/2 - \alpha{a }) & \hbox{if } \alpha > a , \cr 1/\alpha \; \exp (- a^{2}/2) & \hbox{otherwise.} \cr \end{cases} $$ The first expression is minimized by \begin{equation} \alpha^{*} = \frac{1}{2}a + \frac{1}{2} \sqrt{a^2 + 4}\;,\qquad (1) \end{equation} whereas $\tilde\alpha = a $ minimizes the second bound. The optimal choice of $\alpha$ is therefore (1).

Answer 3

Sample from a Normal distribution but ignore all the random values falling outside the specified range before simulations.

This method is correct, but, as mentionned by @Xi'an in his answer, it would take a long time when the range is small (more precisely, when its measure is small under the normal distribution).

As any other distribution, one could use the inversion method $F^{-1}(U)$ (also called inverse transform sampling), where $F$ is the (cumulative function of the) distribution of interest and $U\sim\text{Unif}(0,1)$. When $F$ is the distribution obtained by truncating some distribution $G$ on some interval $(a,b)$, this is equivalent to sample $G^{-1}(U)$ with $U\sim\text{Unif}\bigl(G(a),G(b)\bigr)$.

However, and this is already mentionned by @Xi'an's in a comment, for some situations the inversion method requires a very precise evaluation of the quantile function $G^{-1}$, and I would add it also requires a fast computation of $G^{-1}$. When $G$ is a normal distribution, the evaluation of $G^{-1}$ is rather slow, and it is not highly precise for values of $a$ and $b$ outside the "range" of $G$.

Simulate a truncated distribution using importance sampling

A possibility is to use importance sampling. Consider the case of the standard Gaussian distribution ${\cal N}(0,1)$. Forget the previous notations, now let $G$ be the Cauchy distribution. The two above mentionned requirements are fulfilled for $G$ : one simply has $\boxed{G(q)=\frac{\arctan(q)}{\pi}+\frac12}$ and $\boxed{G^{-1}(q)=\tan\bigl(\pi(q-\frac12)\bigr)}$. Therefore, the truncated Cauchy distribution is easy to sample by the inversion method and it is a good choice of the instrumental variable for importance sampling of the truncated normal distribution.

After a bit of simplifications, sampling $U\sim\text{Unif}\bigl(G(a),G(b)\bigr)$ and taking $G^{-1}(U)$ is equivalent to take $\tan(U')$ with $U'\sim\text{Unif}\bigl(\arctan(a),\arctan(b)\bigr)$:

a <- 1
b <- 5
nsims <- 10^5
sims <- tan(runif(nsims, atan(a), atan(b)))

Now one has to calculate the weight for each sampled value $x_i$, defined as the ratio $\phi(x)/g(x)$ of the two densities up to normalization, hence we can take $$ w(x) = \exp(-x^2/2)(1+x^2), $$ but it could be safer to take the log-weights:

log_w <- -sims^2/2 + log1p(sims^2)
w <- exp(log_w) # unnormalized weights
w <- w/sum(w)

The weighted sample $(x_i,w(x_i))$ allows to estimate the measure of every interval $[u,v]$ under the target distribution, by summing the weights of each sampled value falling inside the interval:

u <- 2; v<- 4
sum(w[sims>u & sims<v])
## [1] 0.1418

This provides an estimate of the target cumulative function. We can quickly get and plot it with the spatsat package:

F <- spatstat::ewcdf(sims,w)
# estimated F:
curve(F(x), from=a-0.1, to=b+0.1)
# true F:
curve((pnorm(x)-pnorm(a))/(pnorm(b)-pnorm(a)), add=TRUE, col="red")

# approximate probability of u<x<v:
F(v)-F(u)
## [1] 0.1418

Of course, the sample $(x_i)$ is definitely not a sample of the target distribution, but of the instrumental Cauchy distribution, and one gets a sample of the target distribution by performing weighted resampling, for instance using the multinomial sampling:

msample <- rmultinom(1, nsims, w)[,1]
resims <- rep(sims, times=msample)
hist(resims) 

mean(resims>u & resims<v)
## [1] 0.1446

Another method: fast inverse transform sampling

Olver and Townsend developed a sampling method for a broad class of continuous distribution. It is implemented in the chebfun2 library for Matlab as well as the ApproxFun library for Julia. I have recently discovered this library and it sounds very promising (not only for random sampling). Basically this is the inversion method but using powerful approximations of the cdf and the inverse cdf. The input is the target density function up to normalization.

The sample is simply generated by the following code:

using ApproxFun
f = Fun(x -> exp(-x.^2./2), [1,5]);
nsims = 10^5;
x = sample(f,nsims);

As checked below, it yields an estimated measure of the interval $[2,4]$ close to the one previously obtained by importance sampling:

sum((x.>2) & (x.<4))/nsims
## 0.14191