1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2003 Ferdinando Ametrano
5
6 This file is part of QuantLib, a free-software/open-source library
7 for financial quantitative analysts and developers - http://quantlib.org/
8
9 QuantLib is free software: you can redistribute it and/or modify it
10 under the terms of the QuantLib license. You should have received a
11 copy of the license along with this program; if not, please email
12 <quantlib-dev@lists.sf.net>. The license is also available online at
13 <http://quantlib.org/license.shtml>.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the license for more details.
18*/
19
20/*! \file binomialdistribution.hpp
21 \brief Binomial distribution
22*/
23
24#ifndef quantlib_binomial_distribution_h
25#define quantlib_binomial_distribution_h
26
27#include <ql/math/factorial.hpp>
28#include <ql/math/beta.hpp>
29
30namespace QuantLib {
31
32 inline Real binomialCoefficientLn(BigNatural n, BigNatural k) {
33
34 QL_REQUIRE(n>=k, "n<k not allowed");
35
36 return Factorial::ln(n)-Factorial::ln(n: k)-Factorial::ln(n: n-k);
37
38 }
39
40 inline Real binomialCoefficient(BigNatural n, BigNatural k) {
41
42 return std::floor(x: 0.5+std::exp(x: binomialCoefficientLn(n, k)));
43
44 }
45
46 //! Binomial probability distribution function
47 /*! formula here ...
48 Given an integer k it returns its probability in a Binomial
49 distribution with parameters p and n.
50 */
51 class BinomialDistribution {
52 public:
53 /*! \deprecated Use `auto` or `decltype` instead.
54 Deprecated in version 1.29.
55 */
56 QL_DEPRECATED
57 typedef Real argument_type;
58
59 /*! \deprecated Use `auto` or `decltype` instead.
60 Deprecated in version 1.29.
61 */
62 QL_DEPRECATED
63 typedef Real result_type;
64
65 BinomialDistribution(Real p, BigNatural n);
66 // function
67 Real operator()(BigNatural k) const;
68 private:
69 BigNatural n_;
70 Real logP_, logOneMinusP_;
71 };
72
73 //! Cumulative binomial distribution function
74 /*! Given an integer k it provides the cumulative probability
75 of observing kk<=k:
76 formula here ...
77
78 */
79 class CumulativeBinomialDistribution {
80 public:
81 /*! \deprecated Use `auto` or `decltype` instead.
82 Deprecated in version 1.29.
83 */
84 QL_DEPRECATED
85 typedef Real argument_type;
86
87 /*! \deprecated Use `auto` or `decltype` instead.
88 Deprecated in version 1.29.
89 */
90 QL_DEPRECATED
91 typedef Real result_type;
92
93 CumulativeBinomialDistribution(Real p, BigNatural n);
94 // function
95 Real operator()(BigNatural k) const {
96 if (k >= n_)
97 return 1.0;
98 else
99 return 1.0 - incompleteBetaFunction(a: k+1, b: n_-k, x: p_);
100 }
101 private:
102 BigNatural n_;
103 Real p_;
104 };
105
106
107 inline BinomialDistribution::BinomialDistribution(Real p,
108 BigNatural n)
109 : n_(n) {
110
111 if (p==0.0) {
112 logP_ = -QL_MAX_REAL;
113 logOneMinusP_ = 0.0;
114 } else if (p==1.0) {
115 logP_ = 0.0;
116 logOneMinusP_ = -QL_MAX_REAL;
117 } else {
118 QL_REQUIRE(p>0, "negative p not allowed");
119 QL_REQUIRE(p<1.0, "p>1.0 not allowed");
120
121 logP_ = std::log(x: p);
122 logOneMinusP_ = std::log(x: 1.0-p);
123 }
124 }
125
126
127 inline
128 CumulativeBinomialDistribution::CumulativeBinomialDistribution(
129 Real p, BigNatural n)
130 : n_(n), p_(p) {
131
132 QL_REQUIRE(p>=0, "negative p not allowed");
133 QL_REQUIRE(p<=1.0, "p>1.0 not allowed");
134
135 }
136
137 inline Real BinomialDistribution::operator()(BigNatural k) const {
138
139 if (k > n_) return 0.0;
140
141 // p==1.0
142 if (logP_==0.0)
143 return (k==n_ ? 1.0 : 0.0);
144 // p==0.0
145 else if (logOneMinusP_==0.0)
146 return (k==0 ? 1.0 : 0.0);
147 else
148 return std::exp(x: binomialCoefficientLn(n: n_, k) +
149 k * logP_ + (n_-k) * logOneMinusP_);
150 }
151
152
153
154 /*! Given an odd integer n and a real number z it returns p such that:
155 1 - CumulativeBinomialDistribution((n-1)/2, n, p) =
156 CumulativeNormalDistribution(z)
157
158 \pre n must be odd
159 */
160 inline Real PeizerPrattMethod2Inversion(Real z, BigNatural n) {
161
162 QL_REQUIRE(n%2==1,
163 "n must be an odd number: " << n << " not allowed");
164
165 Real result = (z/(n+1.0/3.0+0.1/(n+1.0)));
166 result *= result;
167 result = std::exp(x: -result*(n+1.0/6.0));
168 result = 0.5 + (z>0 ? 1 : -1) * std::sqrt(x: (0.25 * (1.0-result)));
169 return result;
170 }
171
172}
173
174
175#endif
176

source code of quantlib/ql/math/distributions/binomialdistribution.hpp