diff --git a/src/Bpp/Numeric/Prob/TruncatedPoissonDistribution.cpp b/src/Bpp/Numeric/Prob/TruncatedPoissonDistribution.cpp new file mode 100644 index 0000000..94ba235 --- /dev/null +++ b/src/Bpp/Numeric/Prob/TruncatedPoissonDistribution.cpp @@ -0,0 +1,143 @@ +// SPDX-FileCopyrightText: The Bio++ Development Group +// +// SPDX-License-Identifier: CECILL-2.1 + +#include "TruncatedPoissonDistribution.h" + +#include +#include + +using namespace bpp; +using namespace std; + +TruncatedPoissonDistribution::TruncatedPoissonDistribution( + double lambda, + size_t maxK +) : + AbstractDiscreteDistribution(maxK, "TruncatedPoisson."), + maxK_(maxK) +{ + if (lambda <= 0.0) + throw Exception("TruncatedPoisson: lambda must be > 0."); + + if (maxK_ == 0) + throw Exception("TruncatedPoisson: maxK must be >= 1."); + + addParameter_(new Parameter("TruncatedPoisson.lambda", lambda, Parameter::R_PLUS_STAR)); + + updateDistribution_(); +} + +TruncatedPoissonDistribution::TruncatedPoissonDistribution( + const TruncatedPoissonDistribution& tpdd +) : + AbstractDiscreteDistribution(tpdd), + maxK_(tpdd.maxK_) +{} + +double TruncatedPoissonDistribution::getLambda() const +{ + return getParameterValue("lambda"); +} + +void TruncatedPoissonDistribution::setLambda(double lambda) +{ + setParameterValue("lambda", lambda); +} + +void TruncatedPoissonDistribution::setMaxK(size_t maxK) +{ + if (maxK == 0) + throw Exception("TruncatedPoisson: maxK must be >= 1."); + + maxK_ = maxK; + updateDistribution_(); +} + +void TruncatedPoissonDistribution::fireParameterChanged(const ParameterList& parameters) +{ + AbstractDiscreteDistribution::fireParameterChanged(parameters); + setLambda(getParameterValue("lambda")); + updateDistribution_(); +} + +void TruncatedPoissonDistribution::updateDistribution_() +{ + distribution_.clear(); + + const double lambda = getLambda(); + + double p = std::exp(-lambda); + distribution_[0.0] = p; + + for (size_t k = 1; k < maxK_; ++k) + { + p *= lambda / static_cast(k); + distribution_[static_cast(k)] = p; + } + + // Compute normalization over truncated support + double sum = 0.0; + for (const auto& kv : distribution_) + sum += kv.second; + + if (sum <= 0.0) + throw Exception("TruncatedPoisson: normalization failed."); + + // Normalize probabilities to sum to 1 + for (auto& kv : distribution_) + kv.second /= sum; +} + +double TruncatedPoissonDistribution::pProb(double x) const +{ + double cdf = 0.0; + for (const auto& kv : distribution_) + { + if (kv.first <= x) + cdf += kv.second; + else + break; + } + return cdf; +} + +double TruncatedPoissonDistribution::qProb(double p) const +{ + if (p < 0.0 || p > 1.0) + throw Exception("TruncatedPoisson::qProb(): p must be in [0,1]."); + + double cdf = 0.0; + for (const auto& kv : distribution_) + { + cdf += kv.second; + if (cdf >= p) + return kv.first; + } + + // Return max category if p ≈ 1 + return distribution_.rbegin()->first; +} + +double TruncatedPoissonDistribution::Expectation(double a) const +{ + double num = 0.0; + double den = 0.0; + + for (const auto& kv : distribution_) + { + const double x = kv.first; + const double p = kv.second; + + if (x < a) + continue; + + num += x * p; + den += p; + } + + if (den == 0.0) + return 0.0; + + return num / den; +} diff --git a/src/Bpp/Numeric/Prob/TruncatedPoissonDistribution.h b/src/Bpp/Numeric/Prob/TruncatedPoissonDistribution.h new file mode 100644 index 0000000..649341f --- /dev/null +++ b/src/Bpp/Numeric/Prob/TruncatedPoissonDistribution.h @@ -0,0 +1,71 @@ +// SPDX-FileCopyrightText: The Bio++ Development Group +// +// SPDX-License-Identifier: CECILL-2.1 + + + +#ifndef BPP_NUMERIC_PROB_TRUNCATEDPOISSONDISTRIBUTION_H +#define BPP_NUMERIC_PROB_TRUNCATEDPOISSONDISTRIBUTION_H + + +#include "AbstractDiscreteDistribution.h" + +namespace bpp +{ + + /** + * @brief Truncated Poisson distribution. + * Usually the Poisson distribution has infinite support, so to make it computationally tractable + * we consider a truncated version of the distribution, where the support is {0...maxK_} + * + * The CDF is computed and normalized such that the sum of probabilities of every category is 1. + * + * @author Roye Tadmor + */ + +class TruncatedPoissonDistribution : + public AbstractDiscreteDistribution +{ +private: + size_t maxK_; // upper truncation + +public: + /** + * @brief Build a new discretized normal distribution. + * @param maxK the number of categories to use. + * @param lambda The lambda parameter. + * + */ + TruncatedPoissonDistribution( + double lambda, + size_t maxK + ); + +TruncatedPoissonDistribution* clone() const { return new TruncatedPoissonDistribution(*this); } + +TruncatedPoissonDistribution(const TruncatedPoissonDistribution&); + +public: + std::string getName() const {return "TruncatedPoisson";} + + double pProb(double x) const; + double qProb(double p) const; + double Expectation(double a) const; + + double getLambda() const; + void setLambda(double lambda); + + size_t getMaxK() const { return maxK_; } + void setMaxK(size_t maxK); + +protected: + void fireParameterChanged(const ParameterList& parameters) override; + +private: + void updateDistribution_(); +}; + +} + + +#endif // BPP_NUMERIC_PROB_TRUNCATEDPOISSONDISTRIBUTION_H