yat/regression/detail/PoissonFitter.cc

Code
Comments
Other
Rev Date Author Line
3658 13 Jul 17 peter 1 // $Id$
3658 13 Jul 17 peter 2
3658 13 Jul 17 peter 3 /*
3658 13 Jul 17 peter 4   Copyright (C) 2017 Peter Johansson
3658 13 Jul 17 peter 5
3658 13 Jul 17 peter 6   This file is part of the yat library, http://dev.thep.lu.se/yat
3658 13 Jul 17 peter 7
3658 13 Jul 17 peter 8   The yat library is free software; you can redistribute it and/or
3658 13 Jul 17 peter 9   modify it under the terms of the GNU General Public License as
3658 13 Jul 17 peter 10   published by the Free Software Foundation; either version 3 of the
3658 13 Jul 17 peter 11   License, or (at your option) any later version.
3658 13 Jul 17 peter 12
3658 13 Jul 17 peter 13   The yat library is distributed in the hope that it will be useful,
3658 13 Jul 17 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
3658 13 Jul 17 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3658 13 Jul 17 peter 16   General Public License for more details.
3658 13 Jul 17 peter 17
3658 13 Jul 17 peter 18   You should have received a copy of the GNU General Public License
3658 13 Jul 17 peter 19   along with yat. If not, see <http://www.gnu.org/licenses/>.
3658 13 Jul 17 peter 20 */
3658 13 Jul 17 peter 21
3658 13 Jul 17 peter 22 #include <config.h>
3658 13 Jul 17 peter 23
3658 13 Jul 17 peter 24 #include "PoissonFitter.h"
3658 13 Jul 17 peter 25
3658 13 Jul 17 peter 26 #include "yat/utility/DiagonalMatrix.h"
3658 13 Jul 17 peter 27
3658 13 Jul 17 peter 28 #include <cassert>
3658 13 Jul 17 peter 29 #include <cmath>
3658 13 Jul 17 peter 30
3658 13 Jul 17 peter 31 namespace theplu {
3658 13 Jul 17 peter 32 namespace yat {
3658 13 Jul 17 peter 33 namespace regression {
3658 13 Jul 17 peter 34 namespace detail {
3658 13 Jul 17 peter 35
3658 13 Jul 17 peter 36   PoissonFitter::PoissonFitter(const utility::Matrix& X,
3658 13 Jul 17 peter 37                                const utility::VectorBase& y,
3658 13 Jul 17 peter 38                                utility::Vector& beta,
3658 13 Jul 17 peter 39                                utility::Matrix& covariance)
3658 13 Jul 17 peter 40     : beta_(beta), covariance_(covariance), X_(X), y_(y)
3658 13 Jul 17 peter 41   {
3658 13 Jul 17 peter 42     init();
3658 13 Jul 17 peter 43     size_t n = y_.size();
3658 13 Jul 17 peter 44     size_t p = X_.columns();
3658 13 Jul 17 peter 45     double olddev = 0;
3658 13 Jul 17 peter 46     double dev = deviate();
3658 13 Jul 17 peter 47     double deltadev = dev - olddev;
3658 13 Jul 17 peter 48     size_t iteration = 0;
3658 13 Jul 17 peter 49     while (std::abs(deltadev) > 1e-6 && iteration < 100) {
3658 13 Jul 17 peter 50       ++iteration;
3662 19 Jul 17 peter 51
3658 13 Jul 17 peter 52       utility::DiagonalMatrix W(mu_);
3658 13 Jul 17 peter 53       // z = eta + (y-mu)/w
3658 13 Jul 17 peter 54       utility::Vector z(eta_);
3658 13 Jul 17 peter 55       for (size_t i=0; i<z.size(); ++i)
3658 13 Jul 17 peter 56         z(i) += (y(i) - mu_(i)) / mu_(i);
3658 13 Jul 17 peter 57
3658 13 Jul 17 peter 58       assert(X_.rows() == n);
3658 13 Jul 17 peter 59       assert(X_.columns() == p);
3658 13 Jul 17 peter 60
3658 13 Jul 17 peter 61       utility::Matrix I = transpose(X_) * (W * X);
3658 13 Jul 17 peter 62       inverse_svd(I, covariance_);
3658 13 Jul 17 peter 63       // avoid matrix * matrix
3658 13 Jul 17 peter 64       beta_ = covariance_ * (transpose(X_) * (W * z));
3658 13 Jul 17 peter 65
3658 13 Jul 17 peter 66       assert(X_.columns() == beta_.size());
3658 13 Jul 17 peter 67       eta_ = X_ * beta_;
3658 13 Jul 17 peter 68       for (size_t i=0; i<eta_.size(); ++i)
3658 13 Jul 17 peter 69         mu_(i) = exp(eta_(i));
3658 13 Jul 17 peter 70
3658 13 Jul 17 peter 71       olddev = dev;
3658 13 Jul 17 peter 72       dev = deviate();
3658 13 Jul 17 peter 73       deltadev = dev - olddev;
3658 13 Jul 17 peter 74     }
3658 13 Jul 17 peter 75   }
3658 13 Jul 17 peter 76
3658 13 Jul 17 peter 77
3658 13 Jul 17 peter 78   double PoissonFitter::deviate(void) const
3658 13 Jul 17 peter 79   {
3658 13 Jul 17 peter 80     double sum = 0;
3658 13 Jul 17 peter 81     for (size_t i=0; i<y_.size(); ++i) {
3662 19 Jul 17 peter 82       if (y_(i) && mu_(i))
3662 19 Jul 17 peter 83         sum += y_(i) * std::log(y_(i) / mu_(i));
3658 13 Jul 17 peter 84     }
3662 19 Jul 17 peter 85     return 2*sum;
3658 13 Jul 17 peter 86   }
3658 13 Jul 17 peter 87
3658 13 Jul 17 peter 88
3658 13 Jul 17 peter 89   void PoissonFitter::init(void)
3658 13 Jul 17 peter 90   {
3658 13 Jul 17 peter 91     mu_ = y_;
3658 13 Jul 17 peter 92     mu_ += sum(y_)/y_.size();
3658 13 Jul 17 peter 93     mu_ *= 0.5;
3658 13 Jul 17 peter 94     eta_.resize(mu_.size());
3658 13 Jul 17 peter 95     for (size_t i=0; i<eta_.size(); ++i)
3658 13 Jul 17 peter 96       eta_(i) = std::log(mu_(i));
3658 13 Jul 17 peter 97
3658 13 Jul 17 peter 98     beta_.resize(X_.columns());
3658 13 Jul 17 peter 99     covariance_.resize(beta_.size(), beta_.size());
3658 13 Jul 17 peter 100   }
3658 13 Jul 17 peter 101
3658 13 Jul 17 peter 102 }}}}

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