3

Given $X$ a random variable with a normal distribution, what is the best procedure to simulate $X|X\in[a;b]\cup[c;d]$, i.e. we want to simulate the truncated normal distribution only on the intervals $[a;b]$ or $[c;d]$.

I have already seen a question asked on this forum in Truncated normal distribution over a union of intervals, but the answer gives no indication of how to simulate a sample of this law.

KACEFMA.
  • 131
  • 3
    I am astonished by your claim that the answer "gives no indication of how to sample," because it includes working code! Please take another look. However, that code implements the (potentially lethally inefficient) rejection method. Instead, just invert the CDF. Alternatively, sample from a mixture of two Normals, each truncated on one of the intervals. – whuber Dec 20 '21 at 18:40

2 Answers2

6

Assuming wlog that the truncated Normal is a standard $\mathcal N(0,1)$ Normal, the density of the distribution to be simulated writes as $$f(x)\propto e^{-x^2/2}\mathbb I_{(a,b)}(x)+e^{-x^2/2}\mathbb I_{(c,d)}(x)$$ or, with the normalising constant $$f(x)= \dfrac{e^{-x^2/2}\mathbb I_{(a,b)}(x)+e^{-x^2/2}\mathbb I_{(c,d)}(x)}{\sqrt{2\pi}[\Phi(b)-\Phi(a)+\Phi(d)-\Phi(c)]}$$ when assuming that $a<b<c<d$. As noted by whuber, this writes as a mixture of two truncated Normal distributions, $$f(x)= \dfrac{\Phi(b)-\Phi(a)}{\Phi(b)-\Phi(a)+\Phi(d)-\Phi(c)}\dfrac{e^{-x^2/2}\mathbb I_{(a,b)}(x)}{\sqrt{2\pi}[\Phi(b)-\Phi(a)]}\\ \qquad\qquad\qquad+\dfrac{\Phi(d)-\Phi(c)}{\Phi(b)-\Phi(a)+\Phi(d)-\Phi(c)}\dfrac{e^{-x^2/2}\mathbb I_{(c,d)}(x)}{\sqrt{2\pi}[\Phi(d)-\Phi(c)]}$$ Since the probability weights $$\dfrac{\Phi(b)-\Phi(a)}{\Phi(b)-\Phi(a)+\Phi(d)-\Phi(c)}\quad\text{and}\quad\dfrac{\Phi(d)-\Phi(c)}{\Phi(b)-\Phi(a)+\Phi(d)-\Phi(c)}$$ are available, one can pick a component (first or second) according to its probability and then simulate the corresponding truncated Normal (for which I once wrote an efficient accept-reject algorithm that has been implemented in an R package).

Xi'an
  • 105,342
3

Honestly, I don't think you need to do anything more sophisticated than simulating from an unconstrained Normal distribution, and then exclude values outside of your intervals (example in R):

mu = 0
sigma = 1
low1 = -1
high1 = -.5
low2 = 0
high2 = 1
x = rnorm(n, mu, sigma)
in_interval = (x > low1 & x < high1) | (x > low2 & x < high2)
trunc_x = x[in_interval]

There are also more efficient, but more complicated solutions, that let you sample directly from an arbitrary distribution.

Eoin
  • 8,997
  • 2
    This can be incredibly inefficient. Try it with a standard Normal distribution truncated to $(-10,-6)\cup(6,10),$ for instance. – whuber Dec 20 '21 at 18:37