yat/statistics/PoissonMixture.h

Code
Comments
Other
Rev Date Author Line
4180 09 Jun 22 peter 1 #ifndef theplu_yat_statistics_poisson_mixture_
4180 09 Jun 22 peter 2 #define theplu_yat_statistics_poisson_mixture_
4180 09 Jun 22 peter 3
4180 09 Jun 22 peter 4 // $Id$
4180 09 Jun 22 peter 5
4180 09 Jun 22 peter 6 /*
4180 09 Jun 22 peter 7   Copyright (C) 2022 Peter Johansson
4180 09 Jun 22 peter 8
4180 09 Jun 22 peter 9   This file is part of the yat library, https://dev.thep.lu.se/yat
4180 09 Jun 22 peter 10
4180 09 Jun 22 peter 11   The yat library is free software; you can redistribute it and/or
4180 09 Jun 22 peter 12   modify it under the terms of the GNU General Public License as
4180 09 Jun 22 peter 13   published by the Free Software Foundation; either version 3 of the
4180 09 Jun 22 peter 14   License, or (at your option) any later version.
4180 09 Jun 22 peter 15
4180 09 Jun 22 peter 16   The yat library is distributed in the hope that it will be useful,
4180 09 Jun 22 peter 17   but WITHOUT ANY WARRANTY; without even the implied warranty of
4180 09 Jun 22 peter 18   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
4180 09 Jun 22 peter 19   General Public License for more details.
4180 09 Jun 22 peter 20
4180 09 Jun 22 peter 21   You should have received a copy of the GNU General Public License
4180 09 Jun 22 peter 22   along with yat. If not, see <https://www.gnu.org/licenses/>.
4180 09 Jun 22 peter 23 */
4180 09 Jun 22 peter 24
4180 09 Jun 22 peter 25 #include <yat/utility/Vector.h>
4180 09 Jun 22 peter 26
4180 09 Jun 22 peter 27 #include <map>
4180 09 Jun 22 peter 28
4180 09 Jun 22 peter 29 namespace theplu {
4180 09 Jun 22 peter 30 namespace yat {
4180 09 Jun 22 peter 31   namespace utility {
4180 09 Jun 22 peter 32     class Matrix;
4180 09 Jun 22 peter 33   }
4180 09 Jun 22 peter 34 namespace statistics {
4180 09 Jun 22 peter 35
4180 09 Jun 22 peter 36   /**
4180 09 Jun 22 peter 37      \brief Infer Poissonian mixed model
4180 09 Jun 22 peter 38
4180 09 Jun 22 peter 39      Data are described as a mixed model of Poissonians \f$ P(X=k)
4180 09 Jun 22 peter 40      \sum \tau_i po_i(k) \f$ where \f$ po_i(k) = \frac{m^k e^{-m}}{k!} \f$
4180 09 Jun 22 peter 41
4180 09 Jun 22 peter 42      \since New in yat 0.20
4180 09 Jun 22 peter 43    */
4180 09 Jun 22 peter 44   class PoissonMixture
4180 09 Jun 22 peter 45   {
4180 09 Jun 22 peter 46   public:
4180 09 Jun 22 peter 47     /**
4180 09 Jun 22 peter 48        Add \c n data point(s)
4180 09 Jun 22 peter 49      */
4180 09 Jun 22 peter 50     void add(unsigned long int k, unsigned long int n=1);
4180 09 Jun 22 peter 51
4180 09 Jun 22 peter 52     /**
4180 09 Jun 22 peter 53        \brief Remove all data points
4180 09 Jun 22 peter 54      */
4180 09 Jun 22 peter 55     void clear(void);
4180 09 Jun 22 peter 56
4180 09 Jun 22 peter 57     /**
4180 09 Jun 22 peter 58        Create \a n sub-models and tune their parameters such that the
4180 09 Jun 22 peter 59        logL is maximized.
4180 09 Jun 22 peter 60      */
4180 09 Jun 22 peter 61     void fit(size_t n);
4180 09 Jun 22 peter 62
4180 09 Jun 22 peter 63     /**
4180 09 Jun 22 peter 64        Initialise the model parameters as \c m and \c tau and
4180 09 Jun 22 peter 65        fit the data by maximizing logL.
4180 09 Jun 22 peter 66      */
4180 09 Jun 22 peter 67     void fit(const yat::utility::VectorBase& m,
4180 09 Jun 22 peter 68              const yat::utility::VectorBase& tau);
4180 09 Jun 22 peter 69
4180 09 Jun 22 peter 70     /**
4180 09 Jun 22 peter 71        Log likelihood is calculated as
4180 09 Jun 22 peter 72        \f$
4181 14 Jun 22 peter 73        \ln L = \sum_s \ln \left( \sum \tau_i po_i(k_s) \right) =
4181 14 Jun 22 peter 74        \sum_s \ln \left( \sum \tau_i \frac{m_i^{k_s} e^{-m_i}}{k_s!} \right) =
4180 09 Jun 22 peter 75        \f$
4180 09 Jun 22 peter 76      */
4180 09 Jun 22 peter 77     double logL(void);
4180 09 Jun 22 peter 78
4180 09 Jun 22 peter 79     /**
4180 09 Jun 22 peter 80        \return mean of model \c i
4180 09 Jun 22 peter 81      */
4180 09 Jun 22 peter 82     double mean(size_t i) const;
4180 09 Jun 22 peter 83
4180 09 Jun 22 peter 84     /**
4180 09 Jun 22 peter 85        \return tau corresponding to model \c i
4180 09 Jun 22 peter 86      */
4180 09 Jun 22 peter 87     double tau(size_t i) const;
4180 09 Jun 22 peter 88   private:
4180 09 Jun 22 peter 89     std::map<unsigned long int, unsigned long int> count_;
4180 09 Jun 22 peter 90     yat::utility::Vector tau_;
4180 09 Jun 22 peter 91     yat::utility::Vector m_;
4180 09 Jun 22 peter 92
4180 09 Jun 22 peter 93     // guess the parameters using a Gaussian mixture model
4180 09 Jun 22 peter 94     void init_fit(size_t n);
4180 09 Jun 22 peter 95     bool expectation_step(yat::utility::Matrix& H) const;
4180 09 Jun 22 peter 96
4180 09 Jun 22 peter 97     // Find the optimal parameters
4180 09 Jun 22 peter 98     void optimize(void);
4180 09 Jun 22 peter 99
4180 09 Jun 22 peter 100     // return true if sum squared difference is > 1e-9
4180 09 Jun 22 peter 101     bool optimize_tau(const yat::utility::Matrix& H);
4180 09 Jun 22 peter 102     // return true if sum squared difference is > 1e-9
4180 09 Jun 22 peter 103     bool optimize_m(const yat::utility::Matrix& H);
4180 09 Jun 22 peter 104   };
4180 09 Jun 22 peter 105
4180 09 Jun 22 peter 106 }}}
4180 09 Jun 22 peter 107 #endif