/src/quantlib/ql/math/distributions/normaldistribution.hpp

Source
/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*
 Copyright (C) 2002, 2003 Ferdinando Ametrano
 Copyright (C) 2000, 2001, 2002, 2003 RiskMap srl
 Copyright (C) 2010 Kakhkhor Abdijalilov

 This file is part of QuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://quantlib.org/

 QuantLib is free software: you can redistribute it and/or modify it
 under the terms of the QuantLib license.  You should have received a
 copy of the license along with this program; if not, please email
 <quantlib-dev@lists.sf.net>. The license is also available online at
 <https://www.quantlib.org/license.shtml>.

 This program is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the license for more details.
*/

/*! \file normaldistribution.hpp
    \brief normal, cumulative and inverse cumulative distributions
*/

#ifndef quantlib_normal_distribution_hpp
#define quantlib_normal_distribution_hpp

#include <ql/math/errorfunction.hpp>
#include <ql/errors.hpp>

namespace QuantLib {

    //! Normal distribution function
    /*! Given x, it returns its probability in a Gaussian normal distribution.
        It provides the first derivative too.

      For average \f$ \mu \f$ and standard deviation \f$ \sigma \f$, the
      density is
      \f[
        f(x) = \frac{1}{\sigma\sqrt{2\pi}}
             \exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right).
      \f]
      The standard deviation must be strictly positive.

        \test the correctness of the returned value is tested by
              checking it against numerical calculations. Cross-checks
              are also performed against the
              CumulativeNormalDistribution and InverseCumulativeNormal
              classes.
    */
    class NormalDistribution {
      public:
        // TODO: Review whether this constructor should remain implicit.
        NormalDistribution(Real average = 0.0,
                           Real sigma = 1.0);
        // function
        Real operator()(Real x) const;
        Real derivative(Real x) const;
      private:
        Real average_, sigma_, normalizationFactor_, denominator_,
            derNormalizationFactor_;
    };

    typedef NormalDistribution GaussianDistribution;


    //! Cumulative normal distribution function
    /*! Given x, it provides an approximation to the cumulative probability
      of a Gaussian normal distribution with average \f$ \mu \f$ and
      standard deviation \f$ \sigma \f$:
      \f[
        F(x) = \frac{1}{\sigma\sqrt{2\pi}}
             \int_{-\infty}^{x}
             \exp\left(-\frac{(t-\mu)^2}{2\sigma^2}\right)dt.
      \f]
      The result is between zero and one, and derivative() returns the
      corresponding normal density. The lower tail uses an asymptotic
      expansion when the direct error-function calculation loses precision.

        For this implementation see M. Abramowitz and I. Stegun,
        Handbook of Mathematical Functions,
        Dover Publications, New York (1972)
    */
    class CumulativeNormalDistribution {
      public:
        // TODO: Review whether this constructor should remain implicit.
        CumulativeNormalDistribution(Real average = 0.0,
                                     Real sigma   = 1.0);
        // function
        Real operator()(Real x) const;
        Real derivative(Real x) const;
      private:
        Real average_, sigma_;
        NormalDistribution gaussian_;
        ErrorFunction errorFunction_;
    };


    //! Inverse cumulative normal distribution function
    /*! Given a probability \f$ p \f$ between zero and one, this class
      provides \f$ y = F^{-1}(p) \f$ such that
      \f[
          F(y) = p,
      \f]
      where \f$ F \f$ is the cumulative normal distribution.

      It uses Acklam's approximation:
      by Peter J. Acklam, University of Oslo, Statistics Division.
      URL: http://home.online.no/~pjacklam/notes/invnorm/index.html

      This class can also be used to generate a gaussian normal
      distribution from a uniform distribution.
      This is especially useful when a gaussian normal distribution
      is generated from a low discrepancy uniform distribution:
      in this case the traditional Box-Muller approach and its
      variants would not preserve the sequence's low-discrepancy.

    */
    class InverseCumulativeNormal {
      public:
        // TODO: Review whether this constructor should remain implicit.
        InverseCumulativeNormal(Real average = 0.0,
                                Real sigma   = 1.0);
        // function
        Real operator()(Real x) const {
            return average_ + sigma_*standard_value(x);
        }
        // value for average=0, sigma=1
        /* Compared to operator(), this method avoids 2 floating point
           operations (we use average=0 and sigma=1 most of the
           time). The speed difference is noticeable.
        */
        static Real standard_value(Real x) {
            Real z;
            if (x < x_low_ || x_high_ < x) {
                z = tail_value(x);
            } else {
                z = x - 0.5;
                Real r = z*z;
                z = (((((a1_*r+a2_)*r+a3_)*r+a4_)*r+a5_)*r+a6_)*z /
                    (((((b1_*r+b2_)*r+b3_)*r+b4_)*r+b5_)*r+1.0);
            }

            // The relative error of the approximation has absolute value less
            // than 1.15e-9.  One iteration of Halley's rational method (third
            // order) gives full machine precision.
            // #define REFINE_TO_FULL_MACHINE_PRECISION_USING_HALLEYS_METHOD
            #ifdef REFINE_TO_FULL_MACHINE_PRECISION_USING_HALLEYS_METHOD
            // error (f_(z) - x) divided by the cumulative's derivative
            const Real r = (f_(z) - x) * M_SQRT2 * M_SQRTPI * exp(0.5 * z*z);
            //  Halley's method
            z -= r/(1+0.5*z*r);
            #endif

            return z;
        }
      private:
        /* Handling tails moved into a separate method, which should
           make the inlining of operator() and standard_value method
           easier. tail_value is called rarely and doesn't need to be
           inlined.
        */
        static Real tail_value(Real x);
        #if defined(QL_PATCH_SOLARIS)
        CumulativeNormalDistribution f_;
        #else
        static const CumulativeNormalDistribution f_;
        #endif
        Real average_, sigma_;
        static const Real a1_;
        static const Real a2_;
        static const Real a3_;
        static const Real a4_;
        static const Real a5_;
        static const Real a6_;
        static const Real b1_;
        static const Real b2_;
        static const Real b3_;
        static const Real b4_;
        static const Real b5_;
        static const Real c1_;
        static const Real c2_;
        static const Real c3_;
        static const Real c4_;
        static const Real c5_;
        static const Real c6_;
        static const Real d1_;
        static const Real d2_;
        static const Real d3_;
        static const Real d4_;
        static const Real x_low_;
        static const Real x_high_;
    };

    // backward compatibility
    typedef InverseCumulativeNormal InvCumulativeNormalDistribution;

    //! Moro Inverse cumulative normal distribution class
    /*! Given x between zero and one as
        the integral value of a gaussian normal distribution
        this class provides the value y such that
        formula here ...

        It uses the Beasley-Springer approximation in the central region,
        with an improved Moro approximation for the tails. See Boris Moro,
        "The Full Monte", 1995, Risk Magazine.

        This class can also be used to generate a gaussian normal
        distribution from a uniform distribution.
        This is especially useful when a gaussian normal distribution
        is generated from a low discrepancy uniform distribution:
        in this case the traditional Box-Muller approach and its
        variants would not preserve the sequence's low-discrepancy.

        Peter J. Acklam's approximation is better and is available
        as QuantLib::InverseCumulativeNormal
    */
    class MoroInverseCumulativeNormal {
      public:
        // TODO: Review whether this constructor should remain implicit.
        MoroInverseCumulativeNormal(Real average = 0.0,
                                    Real sigma   = 1.0);
        // function
        Real operator()(Real x) const;
      private:
        Real average_, sigma_;
        static const Real a0_;
        static const Real a1_;
        static const Real a2_;
        static const Real a3_;
        static const Real b0_;
        static const Real b1_;
        static const Real b2_;
        static const Real b3_;
        static const Real c0_;
        static const Real c1_;
        static const Real c2_;
        static const Real c3_;
        static const Real c4_;
        static const Real c5_;
        static const Real c6_;
        static const Real c7_;
        static const Real c8_;
    };

    //! Maddock's Inverse cumulative normal distribution class
    /*! Given a probability \f$ p \f$ between zero and one, this class
      provides \f$ y = F^{-1}(p) \f$ such that
      \f[
        F(y) = p,
      \f]
      where \f$ F \f$ is the cumulative normal distribution.

        From the boost documentation:
         These functions use a rational approximation devised by
         John Maddock to calculate an initial approximation to the
         result that is accurate to ~10^-19, then only if that has
         insufficient accuracy compared to the epsilon for type double,
         do we clean up the result using Halley iteration.
    */
    class MaddockInverseCumulativeNormal {
      public:
        // TODO: Review whether this constructor should remain implicit.
        MaddockInverseCumulativeNormal(Real average = 0.0,
                                       Real sigma   = 1.0);
        Real operator()(Real x) const;

      private:
        const Real average_, sigma_;
    };

    //! Maddock's cumulative normal distribution class
    class MaddockCumulativeNormal {
      public:
        // TODO: Review whether this constructor should remain implicit.
        MaddockCumulativeNormal(Real average = 0.0,
                                       Real sigma   = 1.0);
        Real operator()(Real x) const;

      private:
        const Real average_, sigma_;
    };


    // inline definitions

    inline NormalDistribution::NormalDistribution(Real average,
                                                  Real sigma)
    : average_(average), sigma_(sigma) {

        QL_REQUIRE(sigma_>0.0,
                   "sigma must be greater than 0.0 ("
                   << sigma_ << " not allowed)");

        normalizationFactor_ = M_SQRT_2*M_1_SQRTPI/sigma_;
        derNormalizationFactor_ = sigma_*sigma_;
        denominator_ = 2.0*derNormalizationFactor_;
    }

    inline Real NormalDistribution::operator()(Real x) const {
        Real deltax = x-average_;
        Real exponent = -(deltax*deltax)/denominator_;
        // debian alpha had some strange problem in the very-low range
        return exponent <= -690.0 ? 0.0 :  // exp(x) < 1.0e-300 anyway
            Real(normalizationFactor_*std::exp(exponent));
    }

    inline Real NormalDistribution::derivative(Real x) const {
        const Real density = (*this)(x);
        return density == 0.0 ? 0.0 :
            (density * (average_ - x)) / derNormalizationFactor_;
    }

    inline CumulativeNormalDistribution::CumulativeNormalDistribution(
                                                 Real average, Real sigma)
    : average_(average), sigma_(sigma) {

        QL_REQUIRE(sigma_>0.0,
                   "sigma must be greater than 0.0 ("
                   << sigma_ << " not allowed)");
    }

    inline Real CumulativeNormalDistribution::derivative(Real x) const {
        Real xn = (x - average_) / sigma_;
        return gaussian_(xn) / sigma_;
    }

    inline InverseCumulativeNormal::InverseCumulativeNormal(
                                                 Real average, Real sigma)
    : average_(average), sigma_(sigma) {

        QL_REQUIRE(sigma_>0.0,
                   "sigma must be greater than 0.0 ("
                   << sigma_ << " not allowed)");
    }

    inline MoroInverseCumulativeNormal::MoroInverseCumulativeNormal(
                                                 Real average, Real sigma)
    : average_(average), sigma_(sigma) {

        QL_REQUIRE(sigma_>0.0,
                   "sigma must be greater than 0.0 ("
                   << sigma_ << " not allowed)");
    }

}


#endif

Coverage Report

Created: 2026-09-28 06:23

Line	Count	Source
1		/* -- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -- */
2
3		/*
4		Copyright (C) 2002, 2003 Ferdinando Ametrano
5		Copyright (C) 2000, 2001, 2002, 2003 RiskMap srl
6		Copyright (C) 2010 Kakhkhor Abdijalilov
7
8		This file is part of QuantLib, a free-software/open-source library
9		for financial quantitative analysts and developers - http://quantlib.org/
10
11		QuantLib is free software: you can redistribute it and/or modify it
12		under the terms of the QuantLib license. You should have received a
13		copy of the license along with this program; if not, please email
14		<quantlib-dev@lists.sf.net>. The license is also available online at
15		<https://www.quantlib.org/license.shtml>.
16
17		This program is distributed in the hope that it will be useful, but WITHOUT
18		ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19		FOR A PARTICULAR PURPOSE. See the license for more details.
20		*/
21
22		/*! \file normaldistribution.hpp
23		\brief normal, cumulative and inverse cumulative distributions
24		*/
25
26		#ifndef quantlib_normal_distribution_hpp
27		#define quantlib_normal_distribution_hpp
28
29		#include <ql/math/errorfunction.hpp>
30		#include <ql/errors.hpp>
31
32		namespace QuantLib {
33
34		//! Normal distribution function
35		/*! Given x, it returns its probability in a Gaussian normal distribution.
36		It provides the first derivative too.
37
38		For average \f$ \mu \f$ and standard deviation \f$ \sigma \f$, the
39		density is
40		\f[
41		f(x) = \frac{1}{\sigma\sqrt{2\pi}}
42		\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right).
43		\f]
44		The standard deviation must be strictly positive.
45
46		\test the correctness of the returned value is tested by
47		checking it against numerical calculations. Cross-checks
48		are also performed against the
49		CumulativeNormalDistribution and InverseCumulativeNormal
50		classes.
51		*/
52		class NormalDistribution {
53		public:
54		// TODO: Review whether this constructor should remain implicit.
55		NormalDistribution(Real average = 0.0,
56		Real sigma = 1.0);
57		// function
58		Real operator()(Real x) const;
59		Real derivative(Real x) const;
60		private:
61		Real average_, sigma_, normalizationFactor_, denominator_,
62		derNormalizationFactor_;
63		};
64
65		typedef NormalDistribution GaussianDistribution;
66
67
68		//! Cumulative normal distribution function
69		/*! Given x, it provides an approximation to the cumulative probability
70		of a Gaussian normal distribution with average \f$ \mu \f$ and
71		standard deviation \f$ \sigma \f$:
72		\f[
73		F(x) = \frac{1}{\sigma\sqrt{2\pi}}
74		\int_{-\infty}^{x}
75		\exp\left(-\frac{(t-\mu)^2}{2\sigma^2}\right)dt.
76		\f]
77		The result is between zero and one, and derivative() returns the
78		corresponding normal density. The lower tail uses an asymptotic
79		expansion when the direct error-function calculation loses precision.
80
81		For this implementation see M. Abramowitz and I. Stegun,
82		Handbook of Mathematical Functions,
83		Dover Publications, New York (1972)
84		*/
85		class CumulativeNormalDistribution {
86		public:
87		// TODO: Review whether this constructor should remain implicit.
88		CumulativeNormalDistribution(Real average = 0.0,
89		Real sigma = 1.0);
90		// function
91		Real operator()(Real x) const;
92		Real derivative(Real x) const;
93		private:
94		Real average_, sigma_;
95		NormalDistribution gaussian_;
96		ErrorFunction errorFunction_;
97		};
98
99
100		//! Inverse cumulative normal distribution function
101		/*! Given a probability \f$ p \f$ between zero and one, this class
102		provides \f$ y = F^{-1}(p) \f$ such that
103		\f[
104		F(y) = p,
105		\f]
106		where \f$ F \f$ is the cumulative normal distribution.
107
108		It uses Acklam's approximation:
109		by Peter J. Acklam, University of Oslo, Statistics Division.
110		URL: http://home.online.no/~pjacklam/notes/invnorm/index.html
111
112		This class can also be used to generate a gaussian normal
113		distribution from a uniform distribution.
114		This is especially useful when a gaussian normal distribution
115		is generated from a low discrepancy uniform distribution:
116		in this case the traditional Box-Muller approach and its
117		variants would not preserve the sequence's low-discrepancy.
118
119		*/
120		class InverseCumulativeNormal {
121		public:
122		// TODO: Review whether this constructor should remain implicit.
123		InverseCumulativeNormal(Real average = 0.0,
124		Real sigma = 1.0);
125		// function
126	2.13M	Real operator()(Real x) const {
127	2.13M	return average_ + sigma_*standard_value(x);
128	2.13M	}
129		// value for average=0, sigma=1
130		/* Compared to operator(), this method avoids 2 floating point
131		operations (we use average=0 and sigma=1 most of the
132		time). The speed difference is noticeable.
133		*/
134	2.13M	static Real standard_value(Real x) {
135	2.13M	Real z;
136	2.13M	if (x < x_low_ \|\| x_high_ < x) {
137	106k	z = tail_value(x);
138	2.02M	} else {
139	2.02M	z = x - 0.5;
140	2.02M	Real r = z*z;
141	2.02M	z = (((((a1_r+a2_)r+a3_)r+a4_)r+a5_)r+a6_)z /
142	2.02M	(((((b1_r+b2_)r+b3_)r+b4_)r+b5_)*r+1.0);
143	2.02M	}
144
145		// The relative error of the approximation has absolute value less
146		// than 1.15e-9. One iteration of Halley's rational method (third
147		// order) gives full machine precision.
148		// #define REFINE_TO_FULL_MACHINE_PRECISION_USING_HALLEYS_METHOD
149		#ifdef REFINE_TO_FULL_MACHINE_PRECISION_USING_HALLEYS_METHOD
150		// error (f_(z) - x) divided by the cumulative's derivative
151		const Real r = (f_(z) - x) * M_SQRT2 * M_SQRTPI * exp(0.5 * z*z);
152		// Halley's method
153		z -= r/(1+0.5zr);
154		#endif
155
156	2.13M	return z;
157	2.13M	}
158		private:
159		/* Handling tails moved into a separate method, which should
160		make the inlining of operator() and standard_value method
161		easier. tail_value is called rarely and doesn't need to be
162		inlined.
163		*/
164		static Real tail_value(Real x);
165		#if defined(QL_PATCH_SOLARIS)
166		CumulativeNormalDistribution f_;
167		#else
168		static const CumulativeNormalDistribution f_;
169		#endif
170		Real average_, sigma_;
171		static const Real a1_;
172		static const Real a2_;
173		static const Real a3_;
174		static const Real a4_;
175		static const Real a5_;
176		static const Real a6_;
177		static const Real b1_;
178		static const Real b2_;
179		static const Real b3_;
180		static const Real b4_;
181		static const Real b5_;
182		static const Real c1_;
183		static const Real c2_;
184		static const Real c3_;
185		static const Real c4_;
186		static const Real c5_;
187		static const Real c6_;
188		static const Real d1_;
189		static const Real d2_;
190		static const Real d3_;
191		static const Real d4_;
192		static const Real x_low_;
193		static const Real x_high_;
194		};
195
196		// backward compatibility
197		typedef InverseCumulativeNormal InvCumulativeNormalDistribution;
198
199		//! Moro Inverse cumulative normal distribution class
200		/*! Given x between zero and one as
201		the integral value of a gaussian normal distribution
202		this class provides the value y such that
203		formula here ...
204
205		It uses the Beasley-Springer approximation in the central region,
206		with an improved Moro approximation for the tails. See Boris Moro,
207		"The Full Monte", 1995, Risk Magazine.
208
209		This class can also be used to generate a gaussian normal
210		distribution from a uniform distribution.
211		This is especially useful when a gaussian normal distribution
212		is generated from a low discrepancy uniform distribution:
213		in this case the traditional Box-Muller approach and its
214		variants would not preserve the sequence's low-discrepancy.
215
216		Peter J. Acklam's approximation is better and is available
217		as QuantLib::InverseCumulativeNormal
218		*/
219		class MoroInverseCumulativeNormal {
220		public:
221		// TODO: Review whether this constructor should remain implicit.
222		MoroInverseCumulativeNormal(Real average = 0.0,
223		Real sigma = 1.0);
224		// function
225		Real operator()(Real x) const;
226		private:
227		Real average_, sigma_;
228		static const Real a0_;
229		static const Real a1_;
230		static const Real a2_;
231		static const Real a3_;
232		static const Real b0_;
233		static const Real b1_;
234		static const Real b2_;
235		static const Real b3_;
236		static const Real c0_;
237		static const Real c1_;
238		static const Real c2_;
239		static const Real c3_;
240		static const Real c4_;
241		static const Real c5_;
242		static const Real c6_;
243		static const Real c7_;
244		static const Real c8_;
245		};
246
247		//! Maddock's Inverse cumulative normal distribution class
248		/*! Given a probability \f$ p \f$ between zero and one, this class
249		provides \f$ y = F^{-1}(p) \f$ such that
250		\f[
251		F(y) = p,
252		\f]
253		where \f$ F \f$ is the cumulative normal distribution.
254
255		From the boost documentation:
256		These functions use a rational approximation devised by
257		John Maddock to calculate an initial approximation to the
258		result that is accurate to ~10^-19, then only if that has
259		insufficient accuracy compared to the epsilon for type double,
260		do we clean up the result using Halley iteration.
261		*/
262		class MaddockInverseCumulativeNormal {
263		public:
264		// TODO: Review whether this constructor should remain implicit.
265		MaddockInverseCumulativeNormal(Real average = 0.0,
266		Real sigma = 1.0);
267		Real operator()(Real x) const;
268
269		private:
270		const Real average_, sigma_;
271		};
272
273		//! Maddock's cumulative normal distribution class
274		class MaddockCumulativeNormal {
275		public:
276		// TODO: Review whether this constructor should remain implicit.
277		MaddockCumulativeNormal(Real average = 0.0,
278		Real sigma = 1.0);
279		Real operator()(Real x) const;
280
281		private:
282		const Real average_, sigma_;
283		};
284
285
286		// inline definitions
287
288		inline NormalDistribution::NormalDistribution(Real average,
289		Real sigma)
290	24.2k	: average_(average), sigma_(sigma) {
291
292	24.2k	QL_REQUIRE(sigma_>0.0,
293	24.2k	"sigma must be greater than 0.0 ("
294	24.2k	<< sigma_ << " not allowed)");
295
296	24.2k	normalizationFactor_ = M_SQRT_2*M_1_SQRTPI/sigma_;
297	24.2k	derNormalizationFactor_ = sigma_*sigma_;
298	24.2k	denominator_ = 2.0*derNormalizationFactor_;
299	24.2k	}
300
301	19.8k	inline Real NormalDistribution::operator()(Real x) const {
302	19.8k	Real deltax = x-average_;
303	19.8k	Real exponent = -(deltax*deltax)/denominator_;
304		// debian alpha had some strange problem in the very-low range
305	19.8k	return exponent <= -690.0 ? 0.0 : // exp(x) < 1.0e-300 anyway
306	19.8k	Real(normalizationFactor_*std::exp(exponent));
307	19.8k	}
308
309	305	inline Real NormalDistribution::derivative(Real x) const {
310	305	const Real density = (*this)(x);
311	305	return density == 0.0 ? 0.0 :
312	305	(density * (average_ - x)) / derNormalizationFactor_;
313	305	}
314
315		inline CumulativeNormalDistribution::CumulativeNormalDistribution(
316		Real average, Real sigma)
317	23.3k	: average_(average), sigma_(sigma) {
318
319	23.3k	QL_REQUIRE(sigma_>0.0,
320	23.3k	"sigma must be greater than 0.0 ("
321	23.3k	<< sigma_ << " not allowed)");
322	23.3k	}
323
324	7.74k	inline Real CumulativeNormalDistribution::derivative(Real x) const {
325	7.74k	Real xn = (x - average_) / sigma_;
326	7.74k	return gaussian_(xn) / sigma_;
327	7.74k	}
328
329		inline InverseCumulativeNormal::InverseCumulativeNormal(
330		Real average, Real sigma)
331	3.84k	: average_(average), sigma_(sigma) {
332
333		QL_REQUIRE(sigma_>0.0,
334	3.84k	"sigma must be greater than 0.0 ("
335	3.84k	<< sigma_ << " not allowed)");
336	3.84k	}
337
338		inline MoroInverseCumulativeNormal::MoroInverseCumulativeNormal(
339		Real average, Real sigma)
340		: average_(average), sigma_(sigma) {
341
342		QL_REQUIRE(sigma_>0.0,
343		"sigma must be greater than 0.0 ("
344		<< sigma_ << " not allowed)");
345		}
346
347		}
348
349
350		#endif