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 |
// $Id$ |
4180 |
09 Jun 22 |
peter |
5 |
|
4180 |
09 Jun 22 |
peter |
6 |
/* |
4180 |
09 Jun 22 |
peter |
Copyright (C) 2022 Peter Johansson |
4180 |
09 Jun 22 |
peter |
8 |
|
4180 |
09 Jun 22 |
peter |
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 |
The yat library is free software; you can redistribute it and/or |
4180 |
09 Jun 22 |
peter |
modify it under the terms of the GNU General Public License as |
4180 |
09 Jun 22 |
peter |
published by the Free Software Foundation; either version 3 of the |
4180 |
09 Jun 22 |
peter |
License, or (at your option) any later version. |
4180 |
09 Jun 22 |
peter |
15 |
|
4180 |
09 Jun 22 |
peter |
The yat library is distributed in the hope that it will be useful, |
4180 |
09 Jun 22 |
peter |
but WITHOUT ANY WARRANTY; without even the implied warranty of |
4180 |
09 Jun 22 |
peter |
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
4180 |
09 Jun 22 |
peter |
General Public License for more details. |
4180 |
09 Jun 22 |
peter |
20 |
|
4180 |
09 Jun 22 |
peter |
You should have received a copy of the GNU General Public License |
4180 |
09 Jun 22 |
peter |
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 |
\brief Infer Poissonian mixed model |
4180 |
09 Jun 22 |
peter |
38 |
|
4180 |
09 Jun 22 |
peter |
Data are described as a mixed model of Poissonians \f$ P(X=k) |
4180 |
09 Jun 22 |
peter |
\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 |
\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 |
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 |
\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 |
Create \a n sub-models and tune their parameters such that the |
4180 |
09 Jun 22 |
peter |
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 |
Initialise the model parameters as \c m and \c tau and |
4180 |
09 Jun 22 |
peter |
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 |
Log likelihood is calculated as |
4180 |
09 Jun 22 |
peter |
\f$ |
4181 |
14 Jun 22 |
peter |
\ln L = \sum_s \ln \left( \sum \tau_i po_i(k_s) \right) = |
4181 |
14 Jun 22 |
peter |
\sum_s \ln \left( \sum \tau_i \frac{m_i^{k_s} e^{-m_i}}{k_s!} \right) = |
4180 |
09 Jun 22 |
peter |
\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 |
\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 |
\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 |
// 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 |
// Find the optimal parameters |
4180 |
09 Jun 22 |
peter |
98 |
void optimize(void); |
4180 |
09 Jun 22 |
peter |
99 |
|
4180 |
09 Jun 22 |
peter |
// 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 |
// 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 |