yat/statistics/NegativeBinomialExtendedMixture.cc

Code
Comments
Other
Rev Date Author Line
4192 11 Aug 22 peter 1 // $Id$
4192 11 Aug 22 peter 2
4192 11 Aug 22 peter 3 /*
4192 11 Aug 22 peter 4   Copyright (C) 2022 Peter Johansson
4192 11 Aug 22 peter 5
4192 11 Aug 22 peter 6   This file is part of the yat library, https://dev.thep.lu.se/yat
4192 11 Aug 22 peter 7
4192 11 Aug 22 peter 8   The yat library is free software; you can redistribute it and/or
4192 11 Aug 22 peter 9   modify it under the terms of the GNU General Public License as
4192 11 Aug 22 peter 10   published by the Free Software Foundation; either version 3 of the
4192 11 Aug 22 peter 11   License, or (at your option) any later version.
4192 11 Aug 22 peter 12
4192 11 Aug 22 peter 13   The yat library is distributed in the hope that it will be useful,
4192 11 Aug 22 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
4192 11 Aug 22 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
4192 11 Aug 22 peter 16   General Public License for more details.
4192 11 Aug 22 peter 17
4192 11 Aug 22 peter 18   You should have received a copy of the GNU General Public License
4192 11 Aug 22 peter 19   along with yat. If not, see <https://www.gnu.org/licenses/>.
4192 11 Aug 22 peter 20 */
4192 11 Aug 22 peter 21
4192 11 Aug 22 peter 22 #include <config.h>
4192 11 Aug 22 peter 23
4192 11 Aug 22 peter 24 #include "NegativeBinomialExtendedMixture.h"
4192 11 Aug 22 peter 25
4192 11 Aug 22 peter 26 #include <yat/statistics/GaussianMixture.h>
4192 11 Aug 22 peter 27 #include <yat/utility/BFGS2.h>
4252 18 Nov 22 peter 28 #include <yat/utility/Exception.h>
4192 11 Aug 22 peter 29 #include <yat/utility/Matrix.h>
4192 11 Aug 22 peter 30
4192 11 Aug 22 peter 31 #include <gsl/gsl_randist.h>
4192 11 Aug 22 peter 32 #include <gsl/gsl_sf_gamma.h>
4192 11 Aug 22 peter 33 #include <gsl/gsl_sf_psi.h>
4192 11 Aug 22 peter 34
4192 11 Aug 22 peter 35 #include <cassert>
4192 11 Aug 22 peter 36 #include <cmath>
4252 18 Nov 22 peter 37 #include <sstream>
4192 11 Aug 22 peter 38 #include <utility>
4192 11 Aug 22 peter 39 #include <vector>
4192 11 Aug 22 peter 40
4192 11 Aug 22 peter 41 namespace theplu {
4192 11 Aug 22 peter 42 namespace yat {
4192 11 Aug 22 peter 43 namespace statistics {
4192 11 Aug 22 peter 44
4192 11 Aug 22 peter 45   void NegativeBinomialExtendedMixture::add(unsigned long int k, double y)
4192 11 Aug 22 peter 46   {
4192 11 Aug 22 peter 47     data_.push_back(std::make_pair(k, y));
4192 11 Aug 22 peter 48   }
4192 11 Aug 22 peter 49
4192 11 Aug 22 peter 50
4192 11 Aug 22 peter 51   const utility::Vector& NegativeBinomialExtendedMixture::alpha(void) const
4192 11 Aug 22 peter 52   {
4192 11 Aug 22 peter 53     return alpha_;
4192 11 Aug 22 peter 54   }
4192 11 Aug 22 peter 55
4192 11 Aug 22 peter 56
4192 11 Aug 22 peter 57   void NegativeBinomialExtendedMixture::clear(void)
4192 11 Aug 22 peter 58   {
4192 11 Aug 22 peter 59     data_.clear();
4192 11 Aug 22 peter 60   }
4192 11 Aug 22 peter 61
4192 11 Aug 22 peter 62
4192 11 Aug 22 peter 63   // Calculate P(Z_i=j) given current param_
4192 11 Aug 22 peter 64   bool NegativeBinomialExtendedMixture::expectation_step(utility::Matrix& H)
4192 11 Aug 22 peter 65   {
4192 11 Aug 22 peter 66     // loop over values
4192 11 Aug 22 peter 67     for (size_t j=0; j<data_.size(); ++j) {
4192 11 Aug 22 peter 68       unsigned long int k = data_[j].first;
4192 11 Aug 22 peter 69       double y = data_[j].second;
4192 11 Aug 22 peter 70       utility::VectorView column = H.column_view(j);
4192 11 Aug 22 peter 71       // loop over models
4192 11 Aug 22 peter 72       for (size_t i=0; i<column.size(); ++i)
4192 11 Aug 22 peter 73         column(i) = tau_(i) * pdf(k, y*m_(i), alpha_(i));
4192 11 Aug 22 peter 74       column *= 1.0 / sum(column);
4192 11 Aug 22 peter 75     }
4192 11 Aug 22 peter 76     return false;
4192 11 Aug 22 peter 77   }
4192 11 Aug 22 peter 78
4192 11 Aug 22 peter 79
4192 11 Aug 22 peter 80   void NegativeBinomialExtendedMixture::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     m_ = m;
4192 11 Aug 22 peter 85     alpha_ = alpha;
4192 11 Aug 22 peter 86     tau_ = tau;
4192 11 Aug 22 peter 87     optimize();
4192 11 Aug 22 peter 88   }
4192 11 Aug 22 peter 89
4192 11 Aug 22 peter 90
4192 11 Aug 22 peter 91   void NegativeBinomialExtendedMixture::fit(size_t n)
4192 11 Aug 22 peter 92   {
4192 11 Aug 22 peter 93     if (data_.empty())
4192 11 Aug 22 peter 94       return;
4192 11 Aug 22 peter 95     init_fit(n);
4192 11 Aug 22 peter 96     optimize();
4192 11 Aug 22 peter 97   }
4192 11 Aug 22 peter 98
4192 11 Aug 22 peter 99
4192 11 Aug 22 peter 100   void NegativeBinomialExtendedMixture::init_fit(size_t n)
4192 11 Aug 22 peter 101   {
4192 11 Aug 22 peter 102     m_.resize(n);
4192 11 Aug 22 peter 103     alpha_.resize(n, 1.1);
4192 11 Aug 22 peter 104     tau_.resize(n, 1.0/n);
4192 11 Aug 22 peter 105
4192 11 Aug 22 peter 106     statistics::GaussianMixture gm;
4192 11 Aug 22 peter 107     for (const auto& x : data_)
4192 11 Aug 22 peter 108       gm.add(x.first / x.second);
4192 11 Aug 22 peter 109     gm.fit(n);
4192 11 Aug 22 peter 110
4192 11 Aug 22 peter 111     for (size_t i=0; i<n; ++i)
4192 11 Aug 22 peter 112       m_(i) = gm.mean(i);
4192 11 Aug 22 peter 113   }
4192 11 Aug 22 peter 114
4192 11 Aug 22 peter 115
4192 11 Aug 22 peter 116   double NegativeBinomialExtendedMixture::logL(void) const
4192 11 Aug 22 peter 117   {
4192 11 Aug 22 peter 118     double sum = 0;
4192 11 Aug 22 peter 119     for (size_t i=0; i<data_.size(); ++i) {
4192 11 Aug 22 peter 120       double L = 0;
4192 11 Aug 22 peter 121       for (size_t j=0; j<tau_.size(); ++j)
4192 11 Aug 22 peter 122         L += tau_(j) * pdf(data_[i].first, data_[i].second * m_(j), alpha_(j));
4192 11 Aug 22 peter 123       if (L)
4192 11 Aug 22 peter 124         sum += std::log(L);
4192 11 Aug 22 peter 125     }
4192 11 Aug 22 peter 126     return sum;
4192 11 Aug 22 peter 127   }
4192 11 Aug 22 peter 128
4192 11 Aug 22 peter 129
4192 11 Aug 22 peter 130   const yat::utility::Vector& NegativeBinomialExtendedMixture::m(void) const
4192 11 Aug 22 peter 131   {
4192 11 Aug 22 peter 132     return m_;
4192 11 Aug 22 peter 133   }
4192 11 Aug 22 peter 134
4192 11 Aug 22 peter 135
4192 11 Aug 22 peter 136   void NegativeBinomialExtendedMixture::optimize(void)
4192 11 Aug 22 peter 137   {
4192 11 Aug 22 peter 138     assert(data_.size());
4192 11 Aug 22 peter 139     utility::BFGS2 solver(tau_.size());
4192 11 Aug 22 peter 140     yat::utility::Matrix H(tau_.size(), data_.size());
4192 11 Aug 22 peter 141     bool changed = true;
4192 11 Aug 22 peter 142     for (size_t i=0; i<100 && changed; ++i) {
4192 11 Aug 22 peter 143       changed = expectation_step(H);
4192 11 Aug 22 peter 144       if (optimize_tau(H))
4192 11 Aug 22 peter 145         changed = true;
4192 11 Aug 22 peter 146       if (optimize_param(H, solver))
4192 11 Aug 22 peter 147         changed=true;
4192 11 Aug 22 peter 148     }
4192 11 Aug 22 peter 149   }
4192 11 Aug 22 peter 150
4192 11 Aug 22 peter 151
4192 11 Aug 22 peter 152   bool NegativeBinomialExtendedMixture::optimize_tau(const utility::Matrix& H)
4192 11 Aug 22 peter 153   {
4192 11 Aug 22 peter 154     yat::utility::Vector prev_tau(tau_);
4192 11 Aug 22 peter 155     // tau is updated as the average
4192 11 Aug 22 peter 156     for (size_t i=0; i<tau_.size(); ++i) {
4192 11 Aug 22 peter 157       tau_(i) = sum(H.row_const_view(i));
4192 11 Aug 22 peter 158     }
4192 11 Aug 22 peter 159     tau_ *= 1.0 / sum(tau_);
4192 11 Aug 22 peter 160     double norm2 = norm2_squared(yat::utility::Vector(tau_-prev_tau));
4192 11 Aug 22 peter 161     return norm2 > 1e-9;
4192 11 Aug 22 peter 162   }
4192 11 Aug 22 peter 163
4192 11 Aug 22 peter 164
4192 11 Aug 22 peter 165   bool NegativeBinomialExtendedMixture::
4192 11 Aug 22 peter 166   optimize_param(const utility::Matrix& H,
4192 11 Aug 22 peter 167                  utility::MultiMinimizerDerivative& solver)
4192 11 Aug 22 peter 168   {
4192 11 Aug 22 peter 169     /*
4192 11 Aug 22 peter 170       Maximize Q given current H matrix
4192 11 Aug 22 peter 171       Q = sum_i sum_j P(Z_j=i | K_j=k_j; params) ln(L(param_i; k_j, i)
4192 11 Aug 22 peter 172         = sum_i sum_j H_ij ln(L(param_i; k_j, i)
4192 11 Aug 22 peter 173       where i runs over models and j over samples
4192 11 Aug 22 peter 174
4192 11 Aug 22 peter 175       L = tau*Gamma(r+k)/(Gamma(k+1)*Gamma(r)) * (1-p)^r * p^k
4192 11 Aug 22 peter 176       ln L = ln(tau) + ln(Gamma(r+k)) - ln(Gamma(k+1)) - ln(Gamma(r)) +
4192 11 Aug 22 peter 177       r*ln(1-p) + k*ln(p)
4192 11 Aug 22 peter 178     */
4192 11 Aug 22 peter 179
4192 11 Aug 22 peter 180     yat::utility::Vector old_m = m_;
4192 11 Aug 22 peter 181     yat::utility::Vector old_alpha = alpha_;
4192 11 Aug 22 peter 182     yat::utility::Vector param(2);
4192 11 Aug 22 peter 183     for (size_t i=0; i<tau_.size(); ++i) {
4192 11 Aug 22 peter 184       yat::utility::VectorConstView h = H.row_const_view(i);
4192 11 Aug 22 peter 185       param(0) = m_(i);
4192 11 Aug 22 peter 186       param(1) = alpha_(i);
4192 11 Aug 22 peter 187       Q func(data_, h);
4252 18 Nov 22 peter 188       size_t m_epochs = 1000;
4252 18 Nov 22 peter 189       solver(param, func, utility::MultiMinimizerDerivative::Gradient(1e-3),
4252 18 Nov 22 peter 190              m_epochs);
4252 18 Nov 22 peter 191       if (solver.epochs() == m_epochs) {
4252 18 Nov 22 peter 192         std::ostringstream msg;
4252 18 Nov 22 peter 193         msg << "NegativeBinomialExtendedMixture: maximal epochs reached";
4252 18 Nov 22 peter 194         throw utility::runtime_error(msg.str());
4252 18 Nov 22 peter 195       }
4192 11 Aug 22 peter 196       m_(i) = param(0);
4192 11 Aug 22 peter 197       alpha_(i) = param(1);
4192 11 Aug 22 peter 198     }
4192 11 Aug 22 peter 199
4192 11 Aug 22 peter 200     double norm2 = norm2_squared(yat::utility::Vector(m_ - old_m)) +
4192 11 Aug 22 peter 201       norm2_squared(yat::utility::Vector(alpha_ - old_alpha));
4192 11 Aug 22 peter 202     return norm2 > 1e-9;
4192 11 Aug 22 peter 203   }
4192 11 Aug 22 peter 204
4192 11 Aug 22 peter 205
4192 11 Aug 22 peter 206   double NegativeBinomialExtendedMixture::p(size_t i) const
4192 11 Aug 22 peter 207   {
4192 11 Aug 22 peter 208     return 1.0 - 1.0/alpha_(i);
4192 11 Aug 22 peter 209   }
4192 11 Aug 22 peter 210
4192 11 Aug 22 peter 211
4192 11 Aug 22 peter 212   double NegativeBinomialExtendedMixture::pdf(unsigned long int k, double m,
4192 11 Aug 22 peter 213                                               double alpha) const
4192 11 Aug 22 peter 214   {
4192 11 Aug 22 peter 215     if (m <= 0) {
4192 11 Aug 22 peter 216       if (m == 0 && k==0)
4192 11 Aug 22 peter 217         return 1.0;
4192 11 Aug 22 peter 218       return 0.0;
4192 11 Aug 22 peter 219     }
4192 11 Aug 22 peter 220
4192 11 Aug 22 peter 221     if (alpha <= 1.0)
4192 11 Aug 22 peter 222       return gsl_ran_poisson_pdf(k, m);
4192 11 Aug 22 peter 223     double p = 1 - 1.0 / alpha;
4192 11 Aug 22 peter 224     double r = m / (alpha - 1);
4192 11 Aug 22 peter 225     return gsl_ran_negative_binomial_pdf(k, 1.0-p, r);
4192 11 Aug 22 peter 226   }
4192 11 Aug 22 peter 227
4192 11 Aug 22 peter 228
4192 11 Aug 22 peter 229   double NegativeBinomialExtendedMixture::r(size_t i) const
4192 11 Aug 22 peter 230   {
4192 11 Aug 22 peter 231     return m_(i) / (alpha_(i)-1.0);
4192 11 Aug 22 peter 232   }
4192 11 Aug 22 peter 233
4192 11 Aug 22 peter 234
4192 11 Aug 22 peter 235   const yat::utility::Vector& NegativeBinomialExtendedMixture::tau(void) const
4192 11 Aug 22 peter 236   {
4192 11 Aug 22 peter 237     return tau_;
4192 11 Aug 22 peter 238   }
4192 11 Aug 22 peter 239
4192 11 Aug 22 peter 240
4192 11 Aug 22 peter 241   NegativeBinomialExtendedMixture::
4192 11 Aug 22 peter 242   Q::Q(const std::vector<std::pair<unsigned long int, double>>& data,
4192 11 Aug 22 peter 243        const yat::utility::VectorBase& h)
4192 11 Aug 22 peter 244     : data_(data), h_(h)
4192 11 Aug 22 peter 245   {
4192 11 Aug 22 peter 246   }
4192 11 Aug 22 peter 247
4192 11 Aug 22 peter 248
4192 11 Aug 22 peter 249   double NegativeBinomialExtendedMixture::Q::operator()
4192 11 Aug 22 peter 250     (const yat::utility::VectorBase& x)
4192 11 Aug 22 peter 251   {
4192 11 Aug 22 peter 252     assert(x.size()==2);
4192 11 Aug 22 peter 253     double m = x(0);
4192 11 Aug 22 peter 254     double alpha = x(1);
4192 11 Aug 22 peter 255     double sum = 0;
4192 11 Aug 22 peter 256
4192 11 Aug 22 peter 257     for (size_t i=0; i<data_.size(); ++i) {
4192 11 Aug 22 peter 258       unsigned long int k = data_[i].first;
4192 11 Aug 22 peter 259       double y = data_[i].second;
4192 11 Aug 22 peter 260       sum -= h_(i) * log_L(m, alpha, k, y);
4192 11 Aug 22 peter 261     }
4192 11 Aug 22 peter 262     return sum;
4192 11 Aug 22 peter 263   }
4192 11 Aug 22 peter 264
4192 11 Aug 22 peter 265
4192 11 Aug 22 peter 266   void NegativeBinomialExtendedMixture::Q::operator()
4192 11 Aug 22 peter 267     (const utility::VectorBase& x, utility::VectorMutable& gradient)
4192 11 Aug 22 peter 268   {
4192 11 Aug 22 peter 269     assert(gradient.size() == 2);
4192 11 Aug 22 peter 270     assert(data_.size() == h_.size());
4192 11 Aug 22 peter 271     gradient.all(0.0);
4192 11 Aug 22 peter 272
4192 11 Aug 22 peter 273     assert(x.size()==2);
4192 11 Aug 22 peter 274     double m = x(0);
4192 11 Aug 22 peter 275     double alpha = x(1);
4192 11 Aug 22 peter 276
4192 11 Aug 22 peter 277     for (size_t i=0; i<data_.size(); ++i) {
4192 11 Aug 22 peter 278       unsigned long int k = data_[i].first;
4192 11 Aug 22 peter 279       double y = data_[i].second;
4192 11 Aug 22 peter 280       gradient -= h_(i) * calc_gradient(m, alpha, k, y);
4192 11 Aug 22 peter 281     }
4192 11 Aug 22 peter 282   }
4192 11 Aug 22 peter 283
4192 11 Aug 22 peter 284
4192 11 Aug 22 peter 285   // return dlnL/dp
4192 11 Aug 22 peter 286   double
4192 11 Aug 22 peter 287   NegativeBinomialExtendedMixture::Q::calculate_dp(double p, double r,
4192 11 Aug 22 peter 288                                                    unsigned long int k) const
4192 11 Aug 22 peter 289   {
4192 11 Aug 22 peter 290     /*
4192 11 Aug 22 peter 291       logL = log(Gamma(yr+k)) - log(Gamma(k+1) - log(Gamma(yr)) -
4192 11 Aug 22 peter 292       yrlog(1-p) + klog(p)
4192 11 Aug 22 peter 293
4192 11 Aug 22 peter 294       dlnL/dp = d yr*ln(1-p)/dp + d k*ln(p)/dp =
4192 11 Aug 22 peter 295               = -yr/(1-p) + k/p
4192 11 Aug 22 peter 296     */
4192 11 Aug 22 peter 297     return k/p - r/(1-p);
4192 11 Aug 22 peter 298   }
4192 11 Aug 22 peter 299
4192 11 Aug 22 peter 300
4192 11 Aug 22 peter 301   double
4192 11 Aug 22 peter 302   NegativeBinomialExtendedMixture::Q::calculate_dr(double p, double r,
4192 11 Aug 22 peter 303                                                    unsigned long int k) const
4192 11 Aug 22 peter 304   {
4192 11 Aug 22 peter 305     /*
4192 11 Aug 22 peter 306       logL = log(Gamma(yr+k)) - log(Gamma(k+1) - log(Gamma(yr)) -
4192 11 Aug 22 peter 307       yrlog(1-p) + klog(p)
4192 11 Aug 22 peter 308
4192 11 Aug 22 peter 309       dlnL/dr = y*d ln(Gamma(y*r+k))/dr - y*d ln(Gamma(y*r))/dr + ln(1-p)
4192 11 Aug 22 peter 310     */
4192 11 Aug 22 peter 311     return gsl_sf_psi(r+k) - gsl_sf_psi(r) + std::log(1.0-p);
4192 11 Aug 22 peter 312   }
4192 11 Aug 22 peter 313
4192 11 Aug 22 peter 314
4192 11 Aug 22 peter 315   /*
4192 11 Aug 22 peter 316     return log L =
4192 11 Aug 22 peter 317     log(Gamma(yr+k)) - log(Gamma(k+1) - log(Gamma(yr)) - yrlog(1-p) + klog(p)
4192 11 Aug 22 peter 318   */
4192 11 Aug 22 peter 319   double
4192 11 Aug 22 peter 320   NegativeBinomialExtendedMixture::Q::log_L(double m, double alpha,
4192 11 Aug 22 peter 321                                             unsigned long int k,
4192 11 Aug 22 peter 322                                             double y) const
4192 11 Aug 22 peter 323   {
4192 11 Aug 22 peter 324     if (alpha > 1.001) {
4192 11 Aug 22 peter 325       double r = y * m / (alpha-1.0);
4192 11 Aug 22 peter 326       double p = 1.0 - 1.0/alpha;
4192 11 Aug 22 peter 327       return gsl_sf_lngamma(r+k) - gsl_sf_lngamma(k+1) - gsl_sf_lngamma(r) +
4192 11 Aug 22 peter 328         r*log(1.0-p) + k*log(p);
4192 11 Aug 22 peter 329     }
4192 11 Aug 22 peter 330
4192 11 Aug 22 peter 331     // alpha = 1+epsilon
4192 11 Aug 22 peter 332     // r = m / (alpha-1) = m / epsilon
4192 11 Aug 22 peter 333     // p = (alpha-1) / alpha = epsilon / (1+epsilon)
4192 11 Aug 22 peter 334     // 1-p = (1+epsilon-epsilon)/(1+epsilon) = 1/(1+epsilon)
4192 11 Aug 22 peter 335
4192 11 Aug 22 peter 336     // Taylor expand around alpha=1, f(1) + f'(1) * (alpha-1)
4192 11 Aug 22 peter 337     // where f(1) is the Poisson logL
4192 11 Aug 22 peter 338     // logL = k*log(m) - m - log(Gamma(k+1)) +
4192 11 Aug 22 peter 339     //    + [(k-1)*k/(2*m) + m/2 - k] (alpha-1)
4192 11 Aug 22 peter 340
4192 11 Aug 22 peter 341     double res = k*std::log(y*m) - y*m - gsl_sf_lngamma(k+1);
4192 11 Aug 22 peter 342     if (alpha > 1.0)
4192 11 Aug 22 peter 343       res += (alpha-1.0) * 0.5*((k-1)*k/(y*m) + y*m - 2*k);
4192 11 Aug 22 peter 344     else // penalizing alpha < 1
4192 11 Aug 22 peter 345       res -= 1000*std::pow(alpha-1.0, 2);
4192 11 Aug 22 peter 346
4192 11 Aug 22 peter 347     return res;
4192 11 Aug 22 peter 348   }
4192 11 Aug 22 peter 349
4192 11 Aug 22 peter 350
4192 11 Aug 22 peter 351   utility::Vector
4192 11 Aug 22 peter 352   NegativeBinomialExtendedMixture::Q::calc_gradient(double m, double alpha,
4192 11 Aug 22 peter 353                                                     unsigned long int k,
4192 11 Aug 22 peter 354                                                     double y) const
4192 11 Aug 22 peter 355   {
4192 11 Aug 22 peter 356     utility::Vector res(2);
4192 11 Aug 22 peter 357     if (alpha > 1.001) {
4192 11 Aug 22 peter 358       double r = m / (alpha-1.0);
4192 11 Aug 22 peter 359       double p = 1.0 - 1.0/alpha;
4192 11 Aug 22 peter 360       double dp = calculate_dp(p, y*r, k);
4192 11 Aug 22 peter 361       double dr = y*calculate_dr(p, y*r, k);
4192 11 Aug 22 peter 362       // df/dm = df/dr * dr/dm + df/dp * dp/dm
4192 11 Aug 22 peter 363       //        = df/dr * 1/(alpha-1)
4192 11 Aug 22 peter 364       res(0) = dr / (alpha-1.0);
4192 11 Aug 22 peter 365
4192 11 Aug 22 peter 366       // df/dalpha = df/dp dp/dalpha + df/dr dr/dalpha =
4192 11 Aug 22 peter 367       //           = df/dp 1/alpha^2 - df/dr m/(alpha-1)^2
4192 11 Aug 22 peter 368       res(1) = dp / std::pow(alpha,2) - dr*m/std::pow(alpha-1, 2);
4192 11 Aug 22 peter 369     }
4192 11 Aug 22 peter 370     else {
4192 11 Aug 22 peter 371       // Taylor expansion at alpha=1:
4192 11 Aug 22 peter 372       // logL = k*log(m*y) - m*y - log(Gamma(k+1)) +
4192 11 Aug 22 peter 373       //    + [(k-1)*k/(2*m*y) + m*y/2 - k] (alpha-1)
4192 11 Aug 22 peter 374
4192 11 Aug 22 peter 375       // df/dm
4192 11 Aug 22 peter 376       res(0) = k/(m*y) - 1;
4192 11 Aug 22 peter 377       if (alpha > 1.0)
4192 11 Aug 22 peter 378         res(0) += (alpha-1.0) * 0.5 * (1.0 - (k-1)*k/(m*m*y));
4192 11 Aug 22 peter 379
4192 11 Aug 22 peter 380       //df/dalpha
4192 11 Aug 22 peter 381       res(1) = (k-1)*k/(2*m*y) + (m*y)/2 - k;
4192 11 Aug 22 peter 382       if (alpha < 1) // penalizing alpha < 1
4192 11 Aug 22 peter 383         res(1) -= 2*1000*(alpha-1);
4192 11 Aug 22 peter 384     }
4192 11 Aug 22 peter 385     return res;
4192 11 Aug 22 peter 386   }
4192 11 Aug 22 peter 387 }}}

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