yat/statistics/ln_pdf.cc

Code
Comments
Other
Rev Date Author Line
4336 14 Apr 23 peter 1 // $Id$
4336 14 Apr 23 peter 2
4336 14 Apr 23 peter 3 /*
4336 14 Apr 23 peter 4   Copyright (C) 2023 Peter Johansson
4336 14 Apr 23 peter 5
4336 14 Apr 23 peter 6   This file is part of the yat library, https://dev.thep.lu.se/trac/yat
4336 14 Apr 23 peter 7
4336 14 Apr 23 peter 8   The yat library is free software; you can redistribute it and/or
4336 14 Apr 23 peter 9   modify it under the terms of the GNU General Public License as
4336 14 Apr 23 peter 10   published by the Free Software Foundation; either version 3 of the
4336 14 Apr 23 peter 11   License, or (at your option) any later version.
4336 14 Apr 23 peter 12
4336 14 Apr 23 peter 13   The yat library is distributed in the hope that it will be useful,
4336 14 Apr 23 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
4336 14 Apr 23 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
4336 14 Apr 23 peter 16   General Public License for more details.
4336 14 Apr 23 peter 17
4336 14 Apr 23 peter 18   You should have received a copy of the GNU General Public License
4336 14 Apr 23 peter 19   along with yat. If not, see <https://www.gnu.org/licenses/>.
4336 14 Apr 23 peter 20 */
4336 14 Apr 23 peter 21
4336 14 Apr 23 peter 22 #include <config.h>
4336 14 Apr 23 peter 23
4336 14 Apr 23 peter 24 #include "ln_pdf.h"
4336 14 Apr 23 peter 25
4336 14 Apr 23 peter 26 #include <gsl/gsl_math.h>
4336 14 Apr 23 peter 27 #include <gsl/gsl_sf_gamma.h>
4336 14 Apr 23 peter 28
4336 14 Apr 23 peter 29 #include <cassert>
4336 14 Apr 23 peter 30 #include <cmath>
4336 14 Apr 23 peter 31
4336 14 Apr 23 peter 32 namespace theplu {
4336 14 Apr 23 peter 33 namespace yat {
4336 14 Apr 23 peter 34 namespace statistics {
4336 14 Apr 23 peter 35
4336 14 Apr 23 peter 36   double bernoulli_ln_pdf(unsigned int k, double p)
4336 14 Apr 23 peter 37   {
4336 14 Apr 23 peter 38     assert(k == 0 || k == 1);
4336 14 Apr 23 peter 39     assert(p > 0.0 && p < 1.0);
4336 14 Apr 23 peter 40     // pdf = p^k * (1-p)^(1-k)
4336 14 Apr 23 peter 41     if (k)
4336 14 Apr 23 peter 42       return std::log(p);
4336 14 Apr 23 peter 43     return std::log(1.0 - p);
4336 14 Apr 23 peter 44   }
4336 14 Apr 23 peter 45
4336 14 Apr 23 peter 46
4336 14 Apr 23 peter 47   double binomial_ln_pdf(unsigned int k, double p, unsigned int n)
4336 14 Apr 23 peter 48   {
4336 14 Apr 23 peter 49     assert(p > 0.0 && p < 1.0);
4336 14 Apr 23 peter 50     assert(k <= n);
4336 14 Apr 23 peter 51
4336 14 Apr 23 peter 52     // pdf = choose(n, k) * p^k * (1-p)^(n-k)
4336 14 Apr 23 peter 53     return gsl_sf_lnchoose(n,k) + k * std::log(p) + (n-k) * std::log(1.0-p);
4336 14 Apr 23 peter 54   }
4336 14 Apr 23 peter 55
4336 14 Apr 23 peter 56
4336 14 Apr 23 peter 57   double exponential_ln_pdf(double x, double m)
4336 14 Apr 23 peter 58   {
4336 14 Apr 23 peter 59     assert(m > 0);
4336 14 Apr 23 peter 60     // pdf = 1.0/m * exp(-x/m)
4336 14 Apr 23 peter 61
4336 14 Apr 23 peter 62     return -std::log(m) - x/m;
4336 14 Apr 23 peter 63   }
4336 14 Apr 23 peter 64
4336 14 Apr 23 peter 65
4336 14 Apr 23 peter 66   double gaussian_ln_pdf(double x, double s)
4336 14 Apr 23 peter 67   {
4336 14 Apr 23 peter 68     assert(s > 0);
4336 14 Apr 23 peter 69     // pdf = 1/(s sqrt(2pi)) * exp(-1/2 * (x/s)^2)
4336 14 Apr 23 peter 70
4337 14 Apr 23 peter 71     return -std::log(s) - 0.5 * (M_LN2 + M_LNPI + std::pow(x/s, 2));
4336 14 Apr 23 peter 72   }
4336 14 Apr 23 peter 73
4336 14 Apr 23 peter 74
4336 14 Apr 23 peter 75   double geometric_ln_pdf(unsigned int k, double p)
4336 14 Apr 23 peter 76   {
4336 14 Apr 23 peter 77     assert(p > 0 && p < 1.0);
4336 14 Apr 23 peter 78     assert(k>0);
4336 14 Apr 23 peter 79     // pdf = p*(1.0-p)^k-1
4336 14 Apr 23 peter 80     return std::log(p) + (k-1)*std::log(1.0-p);
4336 14 Apr 23 peter 81   }
4336 14 Apr 23 peter 82
4336 14 Apr 23 peter 83
4336 14 Apr 23 peter 84   double hypergeometric_ln_pdf(unsigned int k, unsigned int n1,
4336 14 Apr 23 peter 85                                unsigned int n2, unsigned int t)
4336 14 Apr 23 peter 86   {
4336 14 Apr 23 peter 87     assert(k <= n1);
4336 14 Apr 23 peter 88     assert(k <= t);
4336 14 Apr 23 peter 89     assert(t-k <= n2);
4336 14 Apr 23 peter 90     assert(t <= n1+n2);
4336 14 Apr 23 peter 91     // pdf = choose(n1,k)*choose(n2,t-k)/choose(n1+n2,t)
4336 14 Apr 23 peter 92
4336 14 Apr 23 peter 93     return gsl_sf_lnchoose(n1, k) + gsl_sf_lnchoose(n2, t-k) -
4336 14 Apr 23 peter 94       gsl_sf_lnchoose(n1+n2, t);
4336 14 Apr 23 peter 95   }
4336 14 Apr 23 peter 96
4336 14 Apr 23 peter 97
4336 14 Apr 23 peter 98   double negative_binomial_ln_pdf(unsigned int k, double p, double n)
4336 14 Apr 23 peter 99   {
4336 14 Apr 23 peter 100     assert(p > 0 && p < 1.0);
4336 14 Apr 23 peter 101     // pdf = Sigma(n+k) / Sigma(k+1)Sigma(n) * p^n * (1-p)^k
4336 14 Apr 23 peter 102     return gsl_sf_lngamma(n+k) - gsl_sf_lngamma(k+1) - gsl_sf_lngamma(n) +
4336 14 Apr 23 peter 103       n * std::log(p) + k * std::log(1.0-p);
4336 14 Apr 23 peter 104   }
4336 14 Apr 23 peter 105
4336 14 Apr 23 peter 106
4336 14 Apr 23 peter 107   double negative_hypergeometric_ln_pdf(unsigned int k, unsigned int n1,
4336 14 Apr 23 peter 108                                         unsigned int n2, unsigned int t)
4336 14 Apr 23 peter 109   {
4336 14 Apr 23 peter 110     return gsl_sf_lnchoose(k+t+1, k) + gsl_sf_lnchoose(n1+n2-t-k, n1-k) -
4336 14 Apr 23 peter 111       gsl_sf_lnchoose(n1+n2, n1);
4336 14 Apr 23 peter 112
4336 14 Apr 23 peter 113   }
4336 14 Apr 23 peter 114
4336 14 Apr 23 peter 115
4336 14 Apr 23 peter 116   double poisson_ln_pdf(unsigned int k, double m)
4336 14 Apr 23 peter 117   {
4336 14 Apr 23 peter 118     assert(m > 0);
4336 14 Apr 23 peter 119     // pdf = m^k exp(-m) / k!
4336 14 Apr 23 peter 120     return k*std::log(m) - m - gsl_sf_lnfact(k);
4336 14 Apr 23 peter 121   }
4336 14 Apr 23 peter 122
4336 14 Apr 23 peter 123
4336 14 Apr 23 peter 124   double ugaussian_ln_pdf(double x)
4336 14 Apr 23 peter 125   {
4336 14 Apr 23 peter 126     // pdf = 1/(sqrt(2pi)) * exp(-1/2 * x^2)
4336 14 Apr 23 peter 127
4337 14 Apr 23 peter 128     return - 0.5 * (M_LN2 + M_LNPI + x*x);
4336 14 Apr 23 peter 129   }
4336 14 Apr 23 peter 130
4336 14 Apr 23 peter 131
4336 14 Apr 23 peter 132 }}} // of namespace utility, yat, and thep

Generated on Sat Oct 14 02:24:05 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)