yat/statistics/Likelihood.cc

Code
Comments
Other
Rev Date Author Line
4289 02 Feb 23 peter 1 // $Id$
4289 02 Feb 23 peter 2
4289 02 Feb 23 peter 3 /*
4289 02 Feb 23 peter 4   Copyright (C) 2023 Peter Johansson
4289 02 Feb 23 peter 5
4289 02 Feb 23 peter 6   This file is part of the yat library, https://dev.thep.lu.se/trac/yat
4289 02 Feb 23 peter 7
4289 02 Feb 23 peter 8   The yat library is free software; you can redistribute it and/or
4289 02 Feb 23 peter 9   modify it under the terms of the GNU General Public License as
4289 02 Feb 23 peter 10   published by the Free Software Foundation; either version 3 of the
4289 02 Feb 23 peter 11   License, or (at your option) any later version.
4289 02 Feb 23 peter 12
4289 02 Feb 23 peter 13   The yat library is distributed in the hope that it will be useful,
4289 02 Feb 23 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
4289 02 Feb 23 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
4289 02 Feb 23 peter 16   General Public License for more details.
4289 02 Feb 23 peter 17
4289 02 Feb 23 peter 18   You should have received a copy of the GNU General Public License
4289 02 Feb 23 peter 19   along with yat. If not, see <https://www.gnu.org/licenses/>.
4289 02 Feb 23 peter 20 */
4289 02 Feb 23 peter 21
4289 02 Feb 23 peter 22 #include <config.h>
4289 02 Feb 23 peter 23
4289 02 Feb 23 peter 24 #include "Likelihood.h"
4289 02 Feb 23 peter 25
4289 02 Feb 23 peter 26 #include <cassert>
4289 02 Feb 23 peter 27 #include <cmath>
4289 02 Feb 23 peter 28 #include <limits>
4289 02 Feb 23 peter 29
4289 02 Feb 23 peter 30 namespace theplu {
4289 02 Feb 23 peter 31 namespace yat {
4289 02 Feb 23 peter 32 namespace statistics {
4289 02 Feb 23 peter 33
4289 02 Feb 23 peter 34   Likelihood::Likelihood(long double x)
4289 02 Feb 23 peter 35     : ln_x_(std::log(x))
4289 02 Feb 23 peter 36   {
4289 02 Feb 23 peter 37     assert((x==0.0) == zero());
4289 02 Feb 23 peter 38   }
4289 02 Feb 23 peter 39
4289 02 Feb 23 peter 40
4289 02 Feb 23 peter 41   double Likelihood::get(void) const
4289 02 Feb 23 peter 42   {
4289 02 Feb 23 peter 43     return std::exp(ln_x_);
4289 02 Feb 23 peter 44   }
4289 02 Feb 23 peter 45
4289 02 Feb 23 peter 46
4289 02 Feb 23 peter 47   long double& Likelihood::log(void)
4289 02 Feb 23 peter 48   {
4289 02 Feb 23 peter 49     return ln_x_;
4289 02 Feb 23 peter 50   }
4289 02 Feb 23 peter 51
4289 02 Feb 23 peter 52
4289 02 Feb 23 peter 53   const long double& Likelihood::log(void) const
4289 02 Feb 23 peter 54   {
4289 02 Feb 23 peter 55     return ln_x_;
4289 02 Feb 23 peter 56   }
4289 02 Feb 23 peter 57
4289 02 Feb 23 peter 58
4289 02 Feb 23 peter 59   bool Likelihood::zero(void) const
4289 02 Feb 23 peter 60   {
4289 02 Feb 23 peter 61     return std::isinf(ln_x_);
4289 02 Feb 23 peter 62   }
4289 02 Feb 23 peter 63
4289 02 Feb 23 peter 64
4289 02 Feb 23 peter 65   bool Likelihood::operator<(const Likelihood& rhs) const
4289 02 Feb 23 peter 66   {
4289 02 Feb 23 peter 67     return ln_x_ < rhs.ln_x_;
4289 02 Feb 23 peter 68   }
4289 02 Feb 23 peter 69
4289 02 Feb 23 peter 70
4289 02 Feb 23 peter 71   bool Likelihood::operator==(const Likelihood& rhs) const
4289 02 Feb 23 peter 72   {
4289 02 Feb 23 peter 73     return ln_x_ == rhs.ln_x_;
4289 02 Feb 23 peter 74   }
4289 02 Feb 23 peter 75
4289 02 Feb 23 peter 76
4289 02 Feb 23 peter 77   Likelihood& Likelihood::operator+=(const Likelihood& rhs)
4289 02 Feb 23 peter 78   {
4289 02 Feb 23 peter 79     if (rhs.zero())
4289 02 Feb 23 peter 80       return *this;
4289 02 Feb 23 peter 81     if (zero()) {
4289 02 Feb 23 peter 82       *this = rhs;
4289 02 Feb 23 peter 83       return *this;
4289 02 Feb 23 peter 84     }
4289 02 Feb 23 peter 85
4289 02 Feb 23 peter 86     // C = A + B
4289 02 Feb 23 peter 87     // logC = log(A+B) = log(A * (1+B/A)) = logA + log(1+B/A)
4289 02 Feb 23 peter 88     if (ln_x_ > rhs.ln_x_) {
4289 02 Feb 23 peter 89       // ln B/A = lnB - lnA
4289 02 Feb 23 peter 90       // B/A = exp(lnB - lnA)
4289 02 Feb 23 peter 91       double epsilon = std::exp(rhs.ln_x_ - ln_x_);
4289 02 Feb 23 peter 92       assert(epsilon <= 1.0);
4289 02 Feb 23 peter 93       ln_x_ += std::log1p(epsilon);
4289 02 Feb 23 peter 94     }
4289 02 Feb 23 peter 95     else {
4289 02 Feb 23 peter 96       double epsilon = std::exp(ln_x_ - rhs.ln_x_);
4289 02 Feb 23 peter 97       assert(epsilon <= 1.0);
4289 02 Feb 23 peter 98       ln_x_ = rhs.ln_x_ + std::log1p(epsilon);
4289 02 Feb 23 peter 99     }
4289 02 Feb 23 peter 100     return *this;
4289 02 Feb 23 peter 101   }
4289 02 Feb 23 peter 102
4289 02 Feb 23 peter 103
4289 02 Feb 23 peter 104   Likelihood& Likelihood::operator*=(const Likelihood& rhs)
4289 02 Feb 23 peter 105   {
4289 02 Feb 23 peter 106     ln_x_ += rhs.ln_x_;
4289 02 Feb 23 peter 107     return *this;
4289 02 Feb 23 peter 108   }
4289 02 Feb 23 peter 109
4289 02 Feb 23 peter 110
4289 02 Feb 23 peter 111   Likelihood& Likelihood::operator*=(double rhs)
4289 02 Feb 23 peter 112   {
4289 02 Feb 23 peter 113     ln_x_ += std::log(rhs);
4289 02 Feb 23 peter 114     return *this;
4289 02 Feb 23 peter 115   }
4289 02 Feb 23 peter 116
4289 02 Feb 23 peter 117
4335 13 Apr 23 peter 118   Likelihood& Likelihood::operator/=(const Likelihood& rhs)
4335 13 Apr 23 peter 119   {
4335 13 Apr 23 peter 120     ln_x_ -= rhs.ln_x_;
4335 13 Apr 23 peter 121     return *this;
4335 13 Apr 23 peter 122   }
4335 13 Apr 23 peter 123
4335 13 Apr 23 peter 124
4289 02 Feb 23 peter 125   Likelihood& Likelihood::operator/=(double rhs)
4289 02 Feb 23 peter 126   {
4289 02 Feb 23 peter 127     ln_x_ -= std::log(rhs);
4289 02 Feb 23 peter 128     return *this;
4289 02 Feb 23 peter 129   }
4289 02 Feb 23 peter 130
4289 02 Feb 23 peter 131
4289 02 Feb 23 peter 132   Likelihood exp(double x)
4289 02 Feb 23 peter 133   {
4289 02 Feb 23 peter 134     Likelihood L;
4289 02 Feb 23 peter 135     L.log() = x;
4289 02 Feb 23 peter 136     return L;
4289 02 Feb 23 peter 137   }
4289 02 Feb 23 peter 138
4289 02 Feb 23 peter 139   double log(const Likelihood& L)
4289 02 Feb 23 peter 140   {
4289 02 Feb 23 peter 141     return L.log();
4289 02 Feb 23 peter 142   }
4289 02 Feb 23 peter 143
4289 02 Feb 23 peter 144
4289 02 Feb 23 peter 145   Likelihood operator*(double x, const Likelihood& L)
4289 02 Feb 23 peter 146   {
4289 02 Feb 23 peter 147     Likelihood result(L);
4289 02 Feb 23 peter 148     result.log() += std::log(x);
4289 02 Feb 23 peter 149     return result;
4289 02 Feb 23 peter 150   }
4289 02 Feb 23 peter 151
4289 02 Feb 23 peter 152
4289 02 Feb 23 peter 153   Likelihood operator*(const Likelihood& L, double x)
4289 02 Feb 23 peter 154   {
4289 02 Feb 23 peter 155     return x * L;
4289 02 Feb 23 peter 156   }
4289 02 Feb 23 peter 157
4289 02 Feb 23 peter 158
4289 02 Feb 23 peter 159   Likelihood operator/(const Likelihood& L, double x)
4289 02 Feb 23 peter 160   {
4289 02 Feb 23 peter 161     Likelihood result(L);
4289 02 Feb 23 peter 162     result.log() -= std::log(x);
4289 02 Feb 23 peter 163     return result;
4289 02 Feb 23 peter 164   }
4289 02 Feb 23 peter 165
4335 13 Apr 23 peter 166
4335 13 Apr 23 peter 167   Likelihood operator/(const Likelihood& num, const Likelihood& den)
4335 13 Apr 23 peter 168   {
4335 13 Apr 23 peter 169     Likelihood result(num);
4335 13 Apr 23 peter 170     result /= den;
4335 13 Apr 23 peter 171     return result;
4335 13 Apr 23 peter 172   }
4335 13 Apr 23 peter 173
4289 02 Feb 23 peter 174 }}} // of namespace utility, yat, and thep

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