廿TT

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

RcppNumerical の optim_lbfgsで最尤推定がしたい(ポアソン回帰)

今日の川柳

RcppNumerical パッケージには例題的にロジスティック回帰(RcppNumerical/fastLR.cpp at master · yixuan/RcppNumerical · GitHub)が実装されているので、ちょっとだけかえてポアソン回帰のコードを書いてみました。


C++ のコード:

// [[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 PoissonReg: public MFuncGrad
{
private:
  const MapMat X;
  const MapVec Y;
  const int n;
  Eigen::VectorXd xbeta;
  Eigen::VectorXd intensity;
public:
  PoissonReg(const MapMat x_, const MapVec y_) :
  X(x_),
  Y(y_),
  n(X.rows()),
  xbeta(n),
  intensity(n)
  {}
  
  double f_grad(Constvec& beta, Refvec grad)
  {
    // Negative log likelihood
    xbeta.noalias() = X * beta;
    const double yxbeta = Y.dot(xbeta);
    for(int i = 0; i < n; i++)
      intensity[i] = std::exp(xbeta[i]);
    const double f = intensity.sum() - yxbeta;
    
    // Gradient
    grad.noalias() = X.transpose() * (intensity - Y);
    
    return f;
  }
  
  Eigen::VectorXd current_xb() const { return xbeta; }
  Eigen::VectorXd current_p()  const { return intensity; }
};

// [[Rcpp::export]]
Rcpp::List optPoisson(Rcpp::NumericMatrix x, Rcpp::NumericVector y,
                   Rcpp::NumericVector start,
                   double eps_f, double eps_g, int maxit)
{
  const MapMat xx = Rcpp::as<MapMat>(x);
  const MapVec yy = Rcpp::as<MapVec>(y);
  // Negative log likelihood
  PoissonReg nll(xx, yy);
  // Initial guess
  Rcpp::NumericVector b = Rcpp::clone(start);
  MapVec beta(b.begin(), b.length());
  
  double fopt;
  int status = optim_lbfgs(nll, beta, fopt, maxit, eps_f, eps_g);
  if(status < 0)
    Rcpp::warning("algorithm did not converge");
  
  return Rcpp::List::create(
    Rcpp::Named("coefficients")      = b,
    Rcpp::Named("fitted.values")     = nll.current_p(),
    Rcpp::Named("linear.predictors") = nll.current_xb(),
    Rcpp::Named("loglikelihood")     = -fopt,
    Rcpp::Named("converged")         = (status >= 0)
  );
}

コンパイルすると #pragma clang diagnostic pop みたいなエラーメッセージが大量にでるがよくわからないのでとりあえず無視。


R のコード:

library("Rcpp")
library("rbenchmark")
sourceCpp("Poisreg.cpp")
counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)
d.AD <- data.frame(treatment, outcome, counts)
glm.D93 <- glm(counts ~ outcome + treatment, family = poisson())

X <-model.matrix(counts ~ outcome + treatment,data=d.AD)
beta <- rep(0,5) #initial value
out <- optPoisson(X,counts,beta,1e-8,1e-8,100)

glm とだいたい同じ値が求まる。

> print(round(out$coefficients,4))
[1]  3.0445 -0.4543 -0.2930  0.0000  0.0000
> print(round(glm.D93$coefficients,4))
(Intercept)    outcome2    outcome3  treatment2  treatment3 
     3.0445     -0.4543     -0.2930      0.0000      0.0000 
> predict(glm.D93)
       1        2        3        4        5        6        7        8        9 
3.044522 2.590267 2.751535 3.044522 2.590267 2.751535 3.044522 2.590267 2.751535 
> out$linear.predictors
[1] 3.044522 2.590268 2.751535 3.044522 2.590268 2.751535 3.044522 2.590268
[9] 2.751535
> fitted(glm.D93)
       1        2        3        4        5        6        7        8        9 
21.00000 13.33333 15.66667 21.00000 13.33333 15.66667 21.00000 13.33333 15.66667 
> out$fitted.values
[1] 21.00000 13.33334 15.66666 21.00000 13.33334 15.66666 21.00000 13.33334
[9] 15.66666

はやいです。

> benchmark(glm.D93 = glm(counts ~ outcome + treatment, family = poisson()),
+           out = optPoisson(X,counts,beta,1e-8,1e-8,100))
     test replications elapsed relative user.self sys.self user.child sys.child
1 glm.D93          100   0.160      160     0.156    0.001          0         0
2     out          100   0.001        1     0.001    0.000          0         0

abrahamcow.hatenablog.com
前にこんなの書いてたのか。完全に忘れてた。

github.com