廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

RcppNumerical でワイブル分布のパラメータの最尤推定

GitHub - yixuan/RcppNumerical: Rcpp Integration for Numerical Computing Libraries を参考に、optim_lbfgs を使ってワイブル分布のパラメータの最尤推定を試した.

変数 f に対数尤度,


\displaystyle l=\sum\{ \log m -\log \eta  + (m-1)( \log y_i -\log \eta ) - (y/\eta)^m\}

g1, g2 に対数対数の一回微分,


\displaystyle \frac{\partial l}{\partial m}=\sum\{ 1/m-\log(\eta) + ( \log y_i -\log \eta ) - (y_i/\eta)^m \log(y_i/\eta) \}


\displaystyle \frac{\partial l}{\partial \eta}=\sum\{-1/\eta - (m-1)/\eta  - m(y_i^m)/(\eta^{m+1}) \}

を与えています.

// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppNumerical)]]

#include <RcppNumerical.h>

using namespace Numer;

typedef Eigen::Map<Eigen::MatrixXd> MapMat;
typedef Eigen::Map<Eigen::VectorXd> MapVec;

class WeibullLL: public MFuncGrad
{
private:
  const MapVec Y;
public:
  WeibullLL(const MapVec y_) : Y(y_) {}
  
  double f_grad(Constvec& par, Refvec grad){

    double f = 0;
    double g1 = 0;
    double g2 = 0;
    const int n = Y.size();
    for(int i=0;i<n;i++){
      f = f + log(par[0]/par[1])+(par[0]-1)*log(Y[i]/par[1]) - pow((Y[i]/par[1]),par[0]);
      g1 = g1 + 1/par[0] + log(Y[i]/par[1]) - pow(Y[i]/par[1],par[0])*log(Y[i]/par[1]);
      g2 = g2 + -1/par[1] - (par[0]-1)/par[1] + par[0]*(pow(Y[i],par[0])/pow(par[1],par[0]+1));
    }
    grad[0] = -g1;
    grad[1] = -g2;
    return -f;
  }
};

// [[Rcpp::export]]
Rcpp::List weibull_ML(Rcpp::NumericVector y, Rcpp::NumericVector init){
  const MapVec yy = Rcpp::as<MapVec>(y);
  Eigen::VectorXd par(2);
  par[0] = init[0];
  par[1] = init[1];
  double fopt;
  WeibullLL f(yy);
  int res = optim_lbfgs(f, par, fopt);
  return Rcpp::List::create(
    Rcpp::Named("par") = par,
    Rcpp::Named("value") = -fopt,
    Rcpp::Named("status") = res
  );
}

R からはこんな感じに使う

library(Rcpp)
sourceCpp("weibull.cpp")
set.seed(1)
y <- rweibull(10,2,20)

res1 <-MASS::fitdistr(y,"weibull")
res2 <- weibull_ML(y,c(1,1))

結果は R の optim でやったのと一致する.

> res1$estimate
    shape     scale 
 1.922598 18.163144 
> res2$par
[1]  1.922792 18.164729

abrahamcow.hatenablog.com