yat/regression/detail/NegativeBinomialFitter.cc

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

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