From 9a98a7e61faac9a6cd5c5d69d50a2ad674612efa Mon Sep 17 00:00:00 2001 From: roye Date: Tue, 3 Feb 2026 13:59:33 +0200 Subject: [PATCH 1/3] Add truncated poisson distribution --- .../Prob/TruncatedPoissonDistribution.cpp | 139 ++++++++++++++++++ .../Prob/TruncatedPoissonDistribution.h | 65 ++++++++ 2 files changed, 204 insertions(+) create mode 100644 src/Bpp/Numeric/Prob/TruncatedPoissonDistribution.cpp create mode 100644 src/Bpp/Numeric/Prob/TruncatedPoissonDistribution.h diff --git a/src/Bpp/Numeric/Prob/TruncatedPoissonDistribution.cpp b/src/Bpp/Numeric/Prob/TruncatedPoissonDistribution.cpp new file mode 100644 index 0000000..ad9787d --- /dev/null +++ b/src/Bpp/Numeric/Prob/TruncatedPoissonDistribution.cpp @@ -0,0 +1,139 @@ +#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; +} \ No newline at end of file diff --git a/src/Bpp/Numeric/Prob/TruncatedPoissonDistribution.h b/src/Bpp/Numeric/Prob/TruncatedPoissonDistribution.h new file mode 100644 index 0000000..8c30b28 --- /dev/null +++ b/src/Bpp/Numeric/Prob/TruncatedPoissonDistribution.h @@ -0,0 +1,65 @@ +#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 From a3576bb2c06487aeb402c6fdfc6cde928046ceb5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Laurent=20Gu=C3=A9guen?= Date: Mon, 9 Feb 2026 11:27:59 +0100 Subject: [PATCH 2/3] Update TruncatedPoissonDistribution.cpp --- src/Bpp/Numeric/Prob/TruncatedPoissonDistribution.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/Bpp/Numeric/Prob/TruncatedPoissonDistribution.cpp b/src/Bpp/Numeric/Prob/TruncatedPoissonDistribution.cpp index ad9787d..94ba235 100644 --- a/src/Bpp/Numeric/Prob/TruncatedPoissonDistribution.cpp +++ b/src/Bpp/Numeric/Prob/TruncatedPoissonDistribution.cpp @@ -1,3 +1,7 @@ +// SPDX-FileCopyrightText: The Bio++ Development Group +// +// SPDX-License-Identifier: CECILL-2.1 + #include "TruncatedPoissonDistribution.h" #include @@ -136,4 +140,4 @@ double TruncatedPoissonDistribution::Expectation(double a) const return 0.0; return num / den; -} \ No newline at end of file +} From f678bb5c0cb662112c17b7af8f269ca16d58ec8e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Laurent=20Gu=C3=A9guen?= Date: Mon, 9 Feb 2026 11:28:40 +0100 Subject: [PATCH 3/3] Update TruncatedPoissonDistribution.h --- src/Bpp/Numeric/Prob/TruncatedPoissonDistribution.h | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/Bpp/Numeric/Prob/TruncatedPoissonDistribution.h b/src/Bpp/Numeric/Prob/TruncatedPoissonDistribution.h index 8c30b28..649341f 100644 --- a/src/Bpp/Numeric/Prob/TruncatedPoissonDistribution.h +++ b/src/Bpp/Numeric/Prob/TruncatedPoissonDistribution.h @@ -1,3 +1,9 @@ +// SPDX-FileCopyrightText: The Bio++ Development Group +// +// SPDX-License-Identifier: CECILL-2.1 + + + #ifndef BPP_NUMERIC_PROB_TRUNCATEDPOISSONDISTRIBUTION_H #define BPP_NUMERIC_PROB_TRUNCATEDPOISSONDISTRIBUTION_H