Untitled


An Improved Algorithm for Generation of Truncated 
Normal Distributed Random Numbers

Vahe V. Sahakyan

Institute for Informatics and Automation Problems of NAS RA

e-mail: vahe@spir.me

Abstract

In this paper we discuss the computational problems of random numbers generation
distributed by truncated normal distribution. It is shown that the standard methods
and libraries have a limit for truncation point caused by the limit on the smallest
number representable by double precision format. Theoretically the problems arise
starting from the truncation point ≈ 40, but in practical calculations the limit is lower,
starting from ≈ 8.5. An improved method is represented, based on the combination
of two approximation algorithms, which with the represented coefficients has 4.5 times
more coverage interval than the standard methods.

Keywords: Random numbers generation, Truncated normal distribution, Error
function, Inverse error function.

1. Introduction

The need of generation of numbers distributed by truncated normal distribution arises across
many statistical calculations and modelling problems. The family of normal and truncated
normal distributions is defined using the error function, and their calculation is strictly
related to the calculation of the error function and the inverse error function. Unfortunately,
this functions are hard to compute numerically. Particularly, in case of the error function
difficulties are visible in tail regions: it convergences to 0 or 1 values very fast and the floating
point representation of the resulting values lose their accuracy. Obviously, the inverse case
will convergence to +∞ or −∞ values very fast for arguments near 1 or 0. Due to these
issues the generation of numbers distributed by truncated normal distribution requires a
different approach. In this paper we are going to address possible ways to overcome those
difficulties and generate numbers distributed by such distribution.

2. Definitions

The probability density function (PDF) of standard normal distribution is defined over the
interval (−∞, +∞). That function is usually denoted as Φ(0, 1, x) and is represented by the
following formula:

73

Mathematical Problems of Computer Science 42, 73–80, 2014.



74 An Improved Algorithm for Generation of Truncated Normal Distributied Random Numbers

f0HxL

f1HxL

0 1 2 3 4
0.0

0.5

1.0

1.5

2.0
F0HxL

F1HxL

0 1 2 3 4 5
0.0

0.2

0.4

0.6

0.8

1.0

1.2

Φ(0, 1, x) =
1√
2π

e−
x
2

2 .

In this article we will observe a distribution defined only on the positive range [0,∞),
often called a half-normal standard distribution (N+). The normal distribution is symmetric
about the origin, hence, it is enough to examine the half-normal case and reconstruct the
normal distribution from it. The same is true for the truncated case.

PDF of half-normal standard distribution on interval [0, +∞) is defined as follows:

f(x) =

{

2 Φ(0, 1, x) x ≥ 0
0 x < 0

=







√

2
π
e−

x
2

2 x ≥ 0
0 x < 0

.

We will say, that a is distributed by half-normal standard distribution truncated by z
(N+

z
), if a follows N+ and satisfies a ≥ z, where z ≥ 0. The PDF of this distribution will

be:

fz(x) =















√
2

π
e
−

x
2

2

erfc

(

z
√

2

) x > z

0 x ≤ z
, (1)

where erfc is a complementary error function, which is defined as

erfc(x) =
2√
π

∫ ∞

x

e−t
2

dt.

The cumulative density function (CDF) of N+
z

can be easily found by integrating (1) on
the interval (−∞, x], and it will be

Fz(x) =
∫

x

−∞
fz(t) dt =















1−
erfc

(

x
√

2

)

erfc

(

z
√

2

) x > z

0 x ≤ z
. (2)

Fig.  1.    PDF and CDF of half-normal standard distribution and truncated by 1 half-normal
standard distributions.



V. Sahakyan 75

calculated

expected

2 4 6 8 10
0.0

0.2

0.4

0.6

0.8

1.0

1
erfc−1(erfc(x))

using [5] for erfc−1

and [6] for erfc.

3. Numbers Generation

There are two common approaches to generate numbers distributed by an arbitrary distri-
bution. The first one is the Von Neumann’s rejection technique [10]. On truncated normal
distribution this technique has very a low efficiency, because the probability that the gen-
erated number satisfies x > z, where x ∼ N+ is equal to f(z). In case, if z = 5 the
f(5) ≈ 3 ∗ 10−6, so approximately 99.9997% of probes will be rejected. Obviously, in appli-
cations this method is not acceptable.

The second approach is the inverse transform sampling method [2]. It claims that if u
is uniformly distributed on the interval (0, 1), then F −1

z
(u) follows the distribution in which

CDF is Fz. We will take u = 1−u to make the equation simple, and in this case the inverse
of (2) will be

F −1
z

(u) =
√
2 erfc−1

(

erfc

(

z√
2

)

u

)

. (3)

Note, that during the inverse function calculation we are not considering the case when
x < z, because in that case Fz(x) = 0, therefore the probability of that is 0, so it goes
beyond our interests.

As we can see, computation of F −1
z

requires computation of erfc and erfc−1. The argument
range of erfc is [0,∞), therefore, the values are ranged in [0, 1]. The u is a random number
from the range (0, 1), hence, the argument range for erfc−1 is [0, 1].

From a theoretical point of view the expression F −1
z

is correct, however, from com-
putational perspective it could imply numerical errors. For example, if z = 40, then
erfc

(

40√
2

)

≈ 7.3 × 10−350, which is not representable in double precision format [3] and
resulting in 0, therefore also erfc−1(0) = ∞. Figure 2 shows the actual numerical calculation
results. It is clearly visible, that starting from values near 6 the result is ∞, even though
theoretically it is becoming not computable starting from values higher than 28.28. The
break point corresponds to the value z = 6

√
2 ≈ 8.49.

Conclusion from here is that independently from precision and accuracy of calculation,
when erfc and erfc−1 are computed separately, starting from z ≥ 40 the final result will
always be ∞ in double precision representation. Obviously, their calculation should be done
combined.

Fig.  2.   This graphic shows calculated  and  expected  values  of



76 An Improved Algorithm for Generation of Truncated Normal Distributied Random Numbers

4. Numerical Approximations

The usual approach to the computation of erfc, is to divide the axis into segments and use
different polynomials and asymptotic approximations for those segments ([1], [6], [7], [8] and
[9]). Let us have a look at approach discussed in [1]. At first erfc can be approximated in
the following form:

erfc(x) = t exp
(

−x2 + P(t)
)

x > 0, (4)

where

t =
2

2 + x
,

and P(t) is polynomial (0 ≤ t ≤ 1) which is found by Chebyshev polynomial approximation.
The sample coefficients for 28 degree polynomial are shown in Table 1.

Now let us have a look at the calculation of erfc−1. As suggested in [1], we apply the
Halley’s method on g(x) = erfc(x) − y and solve it for g(x) = 0. In that case the iteration
variable will be

xn+1 = xn −
2g(xn)g

′(xn)

2[g′(xn)]
2 −g(xn)g′′(xn)

= xn +
erfc(xn)−y

2√
π
exp(−xn2)−xn(erfc(xn)−y)

.

Here, the following form of approximation is used for the initial value of xn:

erfc−1(x) ≈ R(v), 0 < x < 1, (5)
where R(v) is a rational polynomial and

v =

√

−2 log
(

x

2

)

. (6)

R has the following form and the coefficients are found by rational polynomial approxi-
mation [1]

R(v) =
0.70711v3 + 15.6585v2 + 11.5099v −36.4132

v2 + 22.1444v + 22.3164
.

Table 1. Coefficients of polynomial P(t). Degree column corresponds to degree of variable t in
polynomial P .

The represented form of both of these approximations are not the best ones from accuracy
and performance perspective, but they have very simple and convenient form for further
transformations. Accuracy can also be improved by using higher degree polynomials.

5. Improved Algorithm

As was shown in Section 3, with standard approach and using standard libraries the com-
putation of (3) is impossible for values more than 6

√
2. Here we will describe a combined

method from two algorithms described in the previous section. Taking into consideration
the argument of erfc−1 in (3) and (4), (6) will get



V. Sahakyan 77

degree value degree value

0 −1.265512123484646 14 8.81583262585187∗101
1 9.99999999999995∗10−1 15 −2.449786109546412∗102
2 3.750000000013982∗10−1 16 5.592802242564082∗102
3 8.33333331763132∗10−2 17 −1.037357191852551∗103
4 −8.59374904138571∗10−2 18 1.541879449019999∗103
5 −1.437503619675783∗10−1 19 −1.824712058124052∗103
6 −9.17876834732462∗10−2 20 1.714187964988546∗103
7 2.940960653857621∗10−2 21 −1.272415091307403∗103
8 1.351072101958271∗10−1 22 7.388485218111491∗102
9 9.90519562132564∗10−2 23 −3.293364188005826∗102
10 1.566341639700399∗10−1 24 1.090336950910116∗102
11 −1.389503257698621 25 −2.530009680666059∗101
12 5.910981818651237 26 3.677189095748008
13 −2.552672826806604∗101 27 −2.522015791327478∗10−1

v =

√

√

√

√

√−2 log




t u exp
(

−z2
2
+ P(t)

)

2



 =
√

z2 −2P(t)−2 log(t u) + log(4),

where

t =
4

4 +
√
2z

.

By this simple insertion we combine both computation algorithms and as a result exempt
from exponential part in (4). Based on this computation the algorithm will be defined as
follows:

function InverseTruncatedCDF(z, u)
y ← erfc(z)∗u
t ← 4/(4 + sqrt(2)∗z)
v ← sqrt(z ∗z −2∗P(t)−2∗ log(t∗u) + log(4))
x ← R(v)
for i ← 1, 4 do

a ← erfc(x)−y
b ← 2/sqrt(π)∗ exp(−x∗x)−x∗a
if b = 0 then

break

end if

x ← x + a/b
end for

return x
end function

Table 1: Coefficients of polynomial P(t). Degree column corresponds to degree of variable t
in polynomial P .



78 An Improved Algorithm for Generation of Truncated Normal Distributied Random Numbers

PDF

histogram

9.8 10.0 10.2 10.4 10.6 10.8
0

2

4

6

8

10
PDF

histogram

19.9 20.0 20.1 20.2 20.3 20.4
0

5

10

15

20

6

+
10 and N+20, respectively.

The algorithm consists of maximum 4 loops of Halley’s method. This amount was taken
by applying a series of numerical experiments and the results show that the increase of loops
doesn’t improve the result.

Drawback of this algorithm is that with the increase of z the initial x is also increasing.
With coefficients represented in paper for z = 39 we will have x ≈ 28.2833 and the value of
exp(−x ∗ x) ≈ 3.87 ∗ 10−348, which is already not representable in double precision format.
Because of that the Halley’s method improvement steps are not possible, so accuracy is lost.
Nevertheless, with current coefficients it is covering range [0, 38] with maximum error smaller
than 10−14. The result can be improved by taking better approximation for initial x, which
means using a higher degree polynomial for R.

InverseTruncatedCDF can be used for random number generation following distri-
bution N+

z
. The following algorithm can be used for that purpose:

function RandTruncatedNormal(z)
u ← rand() ⊲ uniformly distributed random number from interval (0, 1)
r ← InverseTruncatedCDF(z, u)
return

√
2∗ r

end function

In Figure 3 the result of numbers generation are shown using the above algorithm.

6. Conclusion and Remarks

The problems related to limits of double precision number representation are making im-
possible the calculation of such formulas in tail regions. The suggested way is one of the
approaches that is improving the coverage range. Another way would be the use of Cheby-
shev multi-variate approximation, but as a rule a better coverage means higher degree of
polynomial and consequently a slow computational time. The choice of algorithm is strictly
dependent on the initial problem and here we represented a relatively light-weight method,
also with a possibility to be improved by increasing degrees of approximation polynomials
P and R.

samples calculated by
RandTruncatedNormal function for N
Fig.  3.   Probability  density  function  (PDF)  and  histogram   of     10



V. Sahakyan 7 9

Refer ences

[1 ] W . H . P r e s s , S . Te u ko ls ky, W . T. V e t t e r lin g , a n d B . P . Fla n n e r y,
Numerical R ecipes: The Art of Scienti¯c Computing, Th ir d E d it io n , N e w Y o r k, 2 0 0 7 .

[2 ] L . D e vr o ye , Non-Uniform R andom Variate Generation, S p r in g e r -V e r la g , N e w Y o r k,
1 9 8 6 .

[3 ] A m e r ic a n N a t io n a l S t a n d a r d s In s t it u t e a n d In s t it u t e o f E le c t r ic a l a n d E le c t r o n ic E n g i-
n e e r s , Standard for B inary F loating-P oint Arithmetic, A N S I/ IE E E S t a n d a r d 7 5 4 -1 9 8 5 ,
1 9 8 5 .

[4 ] M. A b r a m o wit z a n d I. A . S t e g u n , Handbook of M athematical F unctions with F ormulas,
Graphs, and M athematical Tables, N a t io n a l B u r e a u o f S t a n d a r d s A p p lie d Ma t h e m a t ic s
S e r ie s - 5 5 , N e w Y o r k, 1 9 6 4 .

[5 ] F. W . J. Olve r , D . W . L o z ie r , R . F. B o is ve r t a n d C. W . Cla r k, \ E r r o r Fu n c t io n s ,
D a ws o n s a n d Fr e s n e l In t e g r a ls " NIST Handbook of M athematical F unctions, Ca m b r id g e
U n ive r s it y P r e s s , p . 1 6 6 , N e w Y o r k, 2 0 1 0

[6 ] S . L . Mo s h ie r , Ce p h e s Ma t h L ib r a r y. [On lin e ]. A va ila b le : h t t p :/ / www.n e t lib .o r g / c e p h e s
[7 ] J. F. H a r t , e t a l., Computer Approximations, Th e S ia m S e r ie s in A p p lie d Ma t h e m a t ic s ,

N e w Y o r k, 1 9 6 8
[8 ] W . J. Co d y, \ R a t io n a l Ch e b ys h e v a p p r o xim a t io n s fo r t h e e r r o r fu n c t io n " M ath. Comp.,

vo l. 2 3 , p p . 6 3 1 -6 3 7 , 1 9 6 9
[9 ] GN U P r o je c t , GN U S c ie n t i¯ c L ib r a r y. [On lin e ]. A va ila b le :

h t t p s :/ / www.g n u .o r g / s o ft wa r e / g s l/
[1 0 ] J. vo n N e u m a n n , \ V a r io u s t e c h n iqu e s u s e d in c o n n e c t io n wit h r a n d o m d ig it s . Mo n t e

Ca r lo m e t h o d s " , Nat. B ureau Standards, vo l. 1 2 , p p . 3 6 -3 8 , 1 9 5 1 .

Submitted 05.09.2014, accepted 28.11.2014.

Îïñí³Í ëï³Ý¹³ñï ÝáñÙ³É µ³ßËí³Í Ãí»ñÇ ·»Ý»ñ³óÙ³Ý
µ³ñ»É³íí³Í ³É·áñÇÃÙ

ì. ê³Ñ³ÏÛ³Ý

²Ù÷á÷áõÙ

Ðá¹í³ÍáõÙùÝݳñÏí³Í»ÝÏïñí³Íëï³Ý¹³ñïÝáñٳɵ³ßËí³ÍÃí»ñÇ·»Ý»ñ³óÙ³Ý
·áñÍÁÝóóáõÙ ³é³ç³óáÕ Ñ³ßíáÕ³Ï³Ý ËݹÇñÝ»ñÁ: ²å³óáõÛóí»É ¿, áñ ëï³Ý¹³ñï
»Õ³Ý³ÏÝ»ñÁ ¨ ·ñ³¹³ñ³ÝÝ»ñÁ áõÝ»Ý ÏïñÙ³Ý Ï»ïÇ ÁÝïñáõÃÛ³Ý ë³Ñٳݳ÷³ÏáõÙ, áñÁ
å³Ûٳݳíáñí³Í ¿ ÏñÏݳÏÇ ×ßïáõÃÛ³Ùµ Ãí»ñÇ Ý»ñϳ۳óÙ³Ý ë³Ñٳݳ÷³ÏáõÙÝ»ñáí:
î»ë³Ï³Ýáñ»Ý ËݹÇñÝ»ñÝ ³é³ç³ÝáõÙ »Ý ¼ 4 0 -Çó Ù»Í ÏïñÙ³Ý Ï»ïÇ Ñ³Ù³ñ, ë³Ï³ÛÝ
·áñÍݳϳÝáõÙ ³Û¹ ë³ÑÙ³ÝÁ ß³ï ³í»ÉÇ ÷áùñ ¿` ¼ 8 :5 -Çó ëÏë³Í: ²é³ç³ñÏí³Í ¿
Ýáñ »Õ³Ý³Ï, áñÇ ÑÇÙùáõÙ »ñÏáõ Ùáï³ñÏÙ³Ý ³É·áñÇÃÙÝ»ñÇ ÙdzíáñáõÙÝ ¿: ÎÇñ³é»Éáí
Ñá¹í³ÍáõÙÝ»ñϳ۳óí³Í·áñͳÏÇóÝ»ñÁ` ÏïñÙ³ÝÏ»ïÇѳٳñëï³óíáõÙ ¿ 4.5³Ý·³Ù
³í»ÉÇ Ù»Í Ñ³ßí³ñÏ»ÉÇ ïÇñáõÛÃ, ù³Ý ëï³Ý¹³ñï »Õ³Ý³ÏÝ»ñÇ ¹»åùáõÙ:



8 0 An Improved Algorithm for Generation of Truncated Normal Distributied Random Numbers

Óëó÷øåííûé àëãîðèòì ãåíåðàöèè ñëó÷àéíûõ ÷èñåë ñ
óñå÷åííûì íîðìàëüíûì ðàñïðåäåëåíèåì

Â. Ñààêÿí

Àííîòàöèÿ

 ýòîé ðàáîòå ðàññìîòðåíû âû÷èñëèòåëüíûå çàäà÷è ãåíåðàöèè ñëó÷àéíûõ
÷èñåë ñ óñå÷åííûì íîðìàëüíûì ðàñïðåäåëåíèåì.  ðàáîòå ïîêàçàíî, ÷òî
ñòàíäàðòíûå ìåòîäû è áèáëèîòåêè èìåþò îãðàíè÷åíèÿ òî÷êè óñå÷åíèÿ. Ýòè
îãðàíè÷åíèÿ îáóñëîâëåíû ëèìèòîì íà íàèìåíüøåå ÷èñëî, ïðåäñòàâèìîå ñ
äâîéíîé òî÷íîñòüþ.Òåîðåòè÷åñêè, ñëîæíîñòè âîçíèêàþò íà÷èíàÿ ñ óñå÷åíèÿ â
òî÷êå ¼ 4 0 , íî â ïðàêòè÷åñêèõ ðàñ÷åòàõ ïðåäåë íàìíîãî íèæå, íà÷èíàÿ ñ ¼ 8 :5 .
Ïðåäñòàâëåí óñîâåðøåíñòâîâàííûé ìåòîä, ñîçäàííûé íà îñíîâå êîìáèíàöèè
äâóõ àïïðîêñèìàöèîííûõ àëãîðèòìîâ. Â ñðàâíåíèè ñî ñòàíäàðòíûìè ìåòîäàìè,
ïðè ïîêàçàííûõ êîýôôèöèåíòàõ, äàííûé ìåòîä îáåñïå÷èâàåò â 4.5 ðàç áîëüøèé
èíòåðâàë ïîêðûòèÿ.


	7_Vahe_Sahakyan.pdf (p.1-6)
	1_Vahe.pdf (p.1-6)