1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2003 Ferdinando Ametrano
5 Copyright (C) 2004 Walter Penschke
6
7 This file is part of QuantLib, a free-software/open-source library
8 for financial quantitative analysts and developers - http://quantlib.org/
9
10 QuantLib is free software: you can redistribute it and/or modify it
11 under the terms of the QuantLib license. You should have received a
12 copy of the license along with this program; if not, please email
13 <quantlib-dev@lists.sf.net>. The license is also available online at
14 <http://quantlib.org/license.shtml>.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the license for more details.
19*/
20
21/*! \file poissondistribution.hpp
22 \brief Poisson distribution
23*/
24
25#ifndef quantlib_poisson_distribution_hpp
26#define quantlib_poisson_distribution_hpp
27
28#include <ql/math/factorial.hpp>
29#include <ql/math/incompletegamma.hpp>
30
31namespace QuantLib {
32
33 //! Poisson distribution function
34 /*! Given an integer \f$ k \f$, it returns its probability
35 in a Poisson distribution.
36
37 \test the correctness of the returned value is tested by
38 checking it against known good results.
39 */
40 class PoissonDistribution {
41 public:
42 /*! \deprecated Use `auto` or `decltype` instead.
43 Deprecated in version 1.29.
44 */
45 QL_DEPRECATED
46 typedef Real argument_type;
47
48 /*! \deprecated Use `auto` or `decltype` instead.
49 Deprecated in version 1.29.
50 */
51 QL_DEPRECATED
52 typedef Real result_type;
53
54 PoissonDistribution(Real mu);
55 // function
56 Real operator()(BigNatural k) const;
57 private:
58 Real mu_, logMu_;
59 };
60
61
62 //! Cumulative Poisson distribution function
63 /*! This function provides an approximation of the
64 integral of the Poisson distribution.
65
66 For this implementation see
67 "Numerical Recipes in C", 2nd edition,
68 Press, Teukolsky, Vetterling, Flannery, chapter 6
69
70 \test the correctness of the returned value is tested by
71 checking it against known good results.
72 */
73 class CumulativePoissonDistribution {
74 public:
75 /*! \deprecated Use `auto` or `decltype` instead.
76 Deprecated in version 1.29.
77 */
78 QL_DEPRECATED
79 typedef Real argument_type;
80
81 /*! \deprecated Use `auto` or `decltype` instead.
82 Deprecated in version 1.29.
83 */
84 QL_DEPRECATED
85 typedef Real result_type;
86
87 CumulativePoissonDistribution(Real mu) : mu_(mu) {}
88 Real operator()(BigNatural k) const {
89 return 1.0 - incompleteGammaFunction(a: k+1, x: mu_);
90 }
91 private:
92 Real mu_;
93 };
94
95
96 //! Inverse cumulative Poisson distribution function
97 /*! \test the correctness of the returned value is tested by
98 checking it against known good results.
99 */
100 class InverseCumulativePoisson {
101 public:
102 /*! \deprecated Use `auto` or `decltype` instead.
103 Deprecated in version 1.29.
104 */
105 QL_DEPRECATED
106 typedef Real argument_type;
107
108 /*! \deprecated Use `auto` or `decltype` instead.
109 Deprecated in version 1.29.
110 */
111 QL_DEPRECATED
112 typedef Real result_type;
113
114 InverseCumulativePoisson(Real lambda = 1.0);
115 Real operator()(Real x) const;
116 private:
117 Real lambda_;
118 Real calcSummand(BigNatural index) const;
119 };
120
121
122
123 // inline definitions
124
125 inline PoissonDistribution::PoissonDistribution(Real mu)
126 : mu_(mu) {
127
128 QL_REQUIRE(mu_>=0.0,
129 "mu must be non negative (" << mu_ << " not allowed)");
130
131 if (mu_!=0.0) logMu_ = std::log(x: mu_);
132 }
133
134 inline Real PoissonDistribution::operator()(BigNatural k) const {
135 if (mu_==0.0) {
136 if (k==0) return 1.0;
137 else return 0.0;
138 }
139 Real logFactorial = Factorial::ln(n: k);
140 return std::exp(x: k*std::log(x: mu_) - logFactorial - mu_);
141 }
142
143
144 inline InverseCumulativePoisson::InverseCumulativePoisson(Real lambda)
145 : lambda_(lambda) {
146 QL_REQUIRE(lambda_ > 0.0, "lambda must be positive");
147 }
148
149 inline Real InverseCumulativePoisson::operator()(Real x) const {
150 QL_REQUIRE(x >= 0.0 && x <= 1.0,
151 "Inverse cumulative Poisson distribution is "
152 "only defined on the interval [0,1]");
153
154 if (x == 1.0)
155 return QL_MAX_REAL;
156
157 Real sum = 0.0;
158 BigNatural index = 0;
159 while (x > sum) {
160 sum += calcSummand(index);
161 index++;
162 }
163
164 return Real(index-1);
165 }
166
167 inline Real InverseCumulativePoisson::calcSummand(BigNatural index) const {
168 return std::exp(x: -lambda_) * std::pow(x: lambda_, y: Integer(index)) /
169 Factorial::get(n: index);
170 }
171
172}
173
174
175#endif
176

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