yat/statistics/NegativeBinomialExtendedMixture.h

Code
Comments
Other
Rev Date Author Line
4192 11 Aug 22 peter 1 #ifndef theplu_yat_statistics_negative_binomial_extended_mixture_
4192 11 Aug 22 peter 2 #define theplu_yat_statistics_negative_binomial_extended_mixture_
4192 11 Aug 22 peter 3
4192 11 Aug 22 peter 4 // $Id$
4192 11 Aug 22 peter 5
4192 11 Aug 22 peter 6 /*
4192 11 Aug 22 peter 7   Copyright (C) 2022 Peter Johansson
4192 11 Aug 22 peter 8
4192 11 Aug 22 peter 9   This file is part of the yat library, https://dev.thep.lu.se/yat
4192 11 Aug 22 peter 10
4192 11 Aug 22 peter 11   The yat library is free software; you can redistribute it and/or
4192 11 Aug 22 peter 12   modify it under the terms of the GNU General Public License as
4192 11 Aug 22 peter 13   published by the Free Software Foundation; either version 3 of the
4192 11 Aug 22 peter 14   License, or (at your option) any later version.
4192 11 Aug 22 peter 15
4192 11 Aug 22 peter 16   The yat library is distributed in the hope that it will be useful,
4192 11 Aug 22 peter 17   but WITHOUT ANY WARRANTY; without even the implied warranty of
4192 11 Aug 22 peter 18   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
4192 11 Aug 22 peter 19   General Public License for more details.
4192 11 Aug 22 peter 20
4192 11 Aug 22 peter 21   You should have received a copy of the GNU General Public License
4192 11 Aug 22 peter 22   along with yat. If not, see <https://www.gnu.org/licenses/>.
4192 11 Aug 22 peter 23 */
4192 11 Aug 22 peter 24
4192 11 Aug 22 peter 25 #include <yat/utility/MultiMinimizerDerivative.h>
4192 11 Aug 22 peter 26 #include <yat/utility/Vector.h>
4192 11 Aug 22 peter 27
4192 11 Aug 22 peter 28 #include <vector>
4192 11 Aug 22 peter 29
4192 11 Aug 22 peter 30 namespace theplu {
4192 11 Aug 22 peter 31 namespace yat {
4192 11 Aug 22 peter 32   namespace utility {
4192 11 Aug 22 peter 33     class Matrix;
4192 11 Aug 22 peter 34     class MultiMinimizerDerivative;
4192 11 Aug 22 peter 35   }
4192 11 Aug 22 peter 36 namespace statistics {
4192 11 Aug 22 peter 37
4192 11 Aug 22 peter 38   /**
4192 11 Aug 22 peter 39      Data are modeled as
4192 11 Aug 22 peter 40
4192 11 Aug 22 peter 41      \f$ X = y * Z \f$
4192 11 Aug 22 peter 42
4192 11 Aug 22 peter 43      where \f$ y \f$ is a known sample-specific constant and \f$ Z \f$
4192 11 Aug 22 peter 44      is a random variable from a mixture of negative binomial
4192 11 Aug 22 peter 45      distributions.
4192 11 Aug 22 peter 46
4192 11 Aug 22 peter 47      \since New in yat 0.20
4192 11 Aug 22 peter 48    */
4192 11 Aug 22 peter 49   class NegativeBinomialExtendedMixture
4192 11 Aug 22 peter 50   {
4192 11 Aug 22 peter 51   public:
4192 11 Aug 22 peter 52     /**
4192 11 Aug 22 peter 53        \brief Add a data point
4192 11 Aug 22 peter 54
4192 11 Aug 22 peter 55        \param k random variable drawn from a mixture of negative
4192 11 Aug 22 peter 56        binomial distributions.
4192 11 Aug 22 peter 57        \param y sample specific sample constant (see class description)
4192 11 Aug 22 peter 58      */
4192 11 Aug 22 peter 59     void add(unsigned long int k, double y);
4192 11 Aug 22 peter 60
4192 11 Aug 22 peter 61     /**
4192 11 Aug 22 peter 62        The \a alpha values are defined as \f$ alpha = \frac{1}{1-p} \f$
4192 11 Aug 22 peter 63      */
4192 11 Aug 22 peter 64     const utility::Vector& alpha(void) const;
4192 11 Aug 22 peter 65
4192 11 Aug 22 peter 66     /**
4192 11 Aug 22 peter 67        \brief remove all data
4192 11 Aug 22 peter 68      */
4192 11 Aug 22 peter 69     void clear(void);
4192 11 Aug 22 peter 70
4192 11 Aug 22 peter 71     /**
4192 11 Aug 22 peter 72        Fit the data to a model with \a n sub-distributions.
4192 11 Aug 22 peter 73      */
4192 11 Aug 22 peter 74     void fit(size_t n);
4192 11 Aug 22 peter 75
4192 11 Aug 22 peter 76     /**
4192 11 Aug 22 peter 77        Fit the data initiating the search for the optimal model
4192 11 Aug 22 peter 78        parameters at the values passed.
4192 11 Aug 22 peter 79      */
4192 11 Aug 22 peter 80     void fit(const utility::VectorBase& m,
4192 11 Aug 22 peter 81              const utility::VectorBase& alpha,
4192 11 Aug 22 peter 82              const utility::VectorBase& tau);
4192 11 Aug 22 peter 83
4192 11 Aug 22 peter 84     /**
4192 11 Aug 22 peter 85        \return log likelihood of getting data given inferred model
4192 11 Aug 22 peter 86      */
4192 11 Aug 22 peter 87     double logL(void) const;
4192 11 Aug 22 peter 88
4192 11 Aug 22 peter 89     /**
4192 11 Aug 22 peter 90        The mean is calculated as \f$ m = \frac{pr}{1-p} \f$
4192 11 Aug 22 peter 91
4192 11 Aug 22 peter 92        \return the mean of each sub-distribution
4192 11 Aug 22 peter 93      */
4192 11 Aug 22 peter 94     const utility::Vector& m(void) const;
4192 11 Aug 22 peter 95
4192 11 Aug 22 peter 96     /**
4192 11 Aug 22 peter 97        \return tau values in model, describing how much each
4192 11 Aug 22 peter 98        sub-distribution contributes to the total
4192 11 Aug 22 peter 99        distribution.
4192 11 Aug 22 peter 100     */
4192 11 Aug 22 peter 101     const utility::Vector& tau(void) const;
4192 11 Aug 22 peter 102   private:
4192 11 Aug 22 peter 103     std::vector<std::pair<unsigned long int, double>> data_;
4192 11 Aug 22 peter 104
4192 11 Aug 22 peter 105     utility::Vector tau_;
4192 11 Aug 22 peter 106     utility::Vector m_;
4192 11 Aug 22 peter 107     utility::Vector alpha_;
4192 11 Aug 22 peter 108
4192 11 Aug 22 peter 109     double pdf(unsigned long int k, double m, double alpha) const;
4192 11 Aug 22 peter 110
4192 11 Aug 22 peter 111     // p = (alpha-1.0) / alpha = 1.0 - 1.0/alpha
4192 11 Aug 22 peter 112     double p(size_t i) const;
4192 11 Aug 22 peter 113     // r = m*(1-p)/p = m / alpha^2
4192 11 Aug 22 peter 114     double r(size_t i) const;
4192 11 Aug 22 peter 115
4192 11 Aug 22 peter 116     void init_fit(size_t n);
4192 11 Aug 22 peter 117     void optimize(void);
4192 11 Aug 22 peter 118
4192 11 Aug 22 peter 119     bool expectation_step(utility::Matrix& H);
4192 11 Aug 22 peter 120     bool optimize_tau(const utility::Matrix& H);
4192 11 Aug 22 peter 121     bool optimize_param(const utility::Matrix& H,
4192 11 Aug 22 peter 122                         utility::MultiMinimizerDerivative& solver);
4192 11 Aug 22 peter 123
4192 11 Aug 22 peter 124     class Q
4192 11 Aug 22 peter 125     {
4192 11 Aug 22 peter 126     public:
4192 11 Aug 22 peter 127       Q(const std::vector<std::pair<unsigned long int, double>>& data,
4192 11 Aug 22 peter 128         const yat::utility::VectorBase& h);
4192 11 Aug 22 peter 129       double operator()(const yat::utility::VectorBase& x);
4192 11 Aug 22 peter 130       void operator()(const yat::utility::VectorBase& x,
4192 11 Aug 22 peter 131                       yat::utility::VectorMutable& gradient);
4192 11 Aug 22 peter 132
4192 11 Aug 22 peter 133     private:
4192 11 Aug 22 peter 134       const std::vector<std::pair<unsigned long int, double>>& data_;
4192 11 Aug 22 peter 135       const yat::utility::VectorBase& h_;
4192 11 Aug 22 peter 136
4192 11 Aug 22 peter 137       double
4192 11 Aug 22 peter 138       log_L(double m, double alpha, unsigned long int k, double y) const;
4192 11 Aug 22 peter 139
4192 11 Aug 22 peter 140       // return <dlogL/dm, dlogL/dalpha>
4192 11 Aug 22 peter 141       yat::utility::Vector
4192 11 Aug 22 peter 142       calc_gradient(double m, double alpha, unsigned long int k,
4192 11 Aug 22 peter 143                     double y) const;
4192 11 Aug 22 peter 144
4192 11 Aug 22 peter 145       double
4192 11 Aug 22 peter 146       calculate_dp(double p, double r, unsigned long int k) const;
4192 11 Aug 22 peter 147       double
4192 11 Aug 22 peter 148       calculate_dr(double p, double r, unsigned long int k) const;
4192 11 Aug 22 peter 149     };
4192 11 Aug 22 peter 150
4192 11 Aug 22 peter 151   };
4192 11 Aug 22 peter 152
4192 11 Aug 22 peter 153 }}}
4192 11 Aug 22 peter 154 #endif

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