/usr/include/ql/math/distributions/poissondistribution.hpp is in libquantlib0-dev 1.1-2build1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 | /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
/*
Copyright (C) 2003 Ferdinando Ametrano
Copyright (C) 2004 Walter Penschke
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
<http://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 poissondistribution.hpp
\brief Poisson distribution
*/
#ifndef quantlib_poisson_distribution_hpp
#define quantlib_poisson_distribution_hpp
#include <ql/math/factorial.hpp>
#include <ql/math/incompletegamma.hpp>
namespace QuantLib {
//! Poisson distribution function
/*! Given an integer \f$ k \f$, it returns its probability
in a Poisson distribution.
\test the correctness of the returned value is tested by
checking it against known good results.
*/
class PoissonDistribution : public std::unary_function<Real,Real> {
public:
PoissonDistribution(Real mu);
// function
Real operator()(BigNatural k) const;
private:
Real mu_, logMu_;
};
//! Cumulative Poisson distribution function
/*! This function provides an approximation of the
integral of the Poisson distribution.
For this implementation see
"Numerical Recipes in C", 2nd edition,
Press, Teukolsky, Vetterling, Flannery, chapter 6
\test the correctness of the returned value is tested by
checking it against known good results.
*/
class CumulativePoissonDistribution
: public std::unary_function<Real,Real> {
public:
CumulativePoissonDistribution(Real mu) : mu_(mu) {}
Real operator()(BigNatural k) const {
return 1.0 - incompleteGammaFunction(k+1, mu_);
}
private:
Real mu_;
};
//! Inverse cumulative Poisson distribution function
/*! \test the correctness of the returned value is tested by
checking it against known good results.
*/
class InverseCumulativePoisson : public std::unary_function<Real,Real> {
public:
InverseCumulativePoisson(Real lambda = 1.0);
Real operator()(Real x) const;
private:
Real lambda_;
Real calcSummand(BigNatural index) const;
};
// inline definitions
inline PoissonDistribution::PoissonDistribution(Real mu)
: mu_(mu) {
QL_REQUIRE(mu_>=0.0,
"mu must be non negative (" << mu_ << " not allowed)");
if (mu_!=0.0) logMu_ = std::log(mu_);
}
inline Real PoissonDistribution::operator()(BigNatural k) const {
if (mu_==0.0) {
if (k==0) return 1.0;
else return 0.0;
}
Real logFactorial = Factorial::ln(k);
return std::exp(k*std::log(mu_) - logFactorial - mu_);
}
inline InverseCumulativePoisson::InverseCumulativePoisson(Real lambda)
: lambda_(lambda) {
QL_REQUIRE(lambda_ > 0.0, "lambda must be positive");
}
inline Real InverseCumulativePoisson::operator()(Real x) const {
QL_REQUIRE(x >= 0.0 && x <= 1.0,
"Inverse cumulative Poisson distribution is "
"only defined on the interval [0,1]");
if (x == 1.0)
return QL_MAX_REAL;
Real sum = 0.0;
BigNatural index = 0;
while (x > sum) {
sum += calcSummand(index);
index++;
}
return Real(index-1);
}
inline Real InverseCumulativePoisson::calcSummand(BigNatural index) const {
return std::exp(-lambda_) * std::pow(lambda_, Integer(index)) /
Factorial::get(index);
}
}
#endif
|