廿TT

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

ガンマ分布の近似的な最尤推定量を用いた PELT アルゴリズムで変化点の検出

今日の川柳

Closed-Form Estimators for the Gamma Distribution Derived From Likelihood Equationsという論文があります。
https://minerva.it.manchester.ac.uk/~saralees/gammapaper.pdf

ガンマ分布のパラメータの新しい推定量を提案しています。

その推定量は閉じた形で求まります。

f:id:abrahamcow:20190606003032p:plain

導出はトリッキーでなかなかおもしろいんだけど、いまどき最尤法でやっても(もっというとベイズでやっても)たいして計算時間かからないので応用先はあんまりないような気がした。

でも尤度を使って変化点の検出とかをやりたい場合、全部の点に対して推定を行わなければいけないので、最尤法の計算コストも馬鹿にならないかもしれない。

なのでPELTアルゴリズム[Rcpp]PELT アルゴリズムで変化点の検出 - 廿TT)とこの推定量を合体させて、Rcppで書いてみました。

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector YZmethod(NumericVector y){
  NumericVector ly = log(y);
  double u1 = mean(y*ly);
  double u2 = mean(ly)*mean(y);
  double bhat = u1 - u2;
  double ahat = mean(y)/bhat;
  NumericVector param(2);
  param[0] = ahat;
  param[1] = bhat;
  return param;
}

double cost_gamma(NumericVector y, int sta, int en) {
  NumericVector suby = y[seq(sta,en)];
  NumericVector param = YZmethod(suby);
  return -sum(dgamma(suby,param[0],param[1],true));
}

// [[Rcpp::export]]
IntegerVector PELTgamma(NumericVector y){
  int n = y.length();
  NumericVector Fv(n+1);
  double pen = 2*log(n);
  Fv[0] = - pen;
  Fv[1] = 0;
  IntegerVector cpts(n);
  cpts[1] = 0;
  IntegerVector R(1);
  IntegerVector cpt_cands;
  int m;
  int tau;
  NumericVector subFv;
  LogicalVector ineq_prune;
  for(int i=1; i<n; i++){
    cpt_cands = R;
    m = cpt_cands.length();
    NumericVector seg_costs(m);
    NumericVector tmp1(m);
    for(int j=0; j<m; j++){
      seg_costs[j] = cost_gamma(y, cpt_cands[j], i);
      tmp1[j] = Fv[cpt_cands[j]] + seg_costs[j] + pen;
    }
    tau = which_min(tmp1);
    Fv[i+1] = tmp1[tau];
    cpts[i] = cpt_cands[tau];
    subFv = Fv[cpt_cands];
    ineq_prune = subFv + seg_costs < Fv[i+1];
    R =  cpt_cands[ineq_prune];
    R.push_back(i);
  }
  int last = cpts[n-1];
  IntegerVector CP(1);
  CP[0]=last;
  while(last > 0){
    last = cpts[last-1];
    CP.push_back(last);
  }
  CP.sort();
  CP.erase(0);
  return CP;
}

乱数でためしてみます。

library(Rcpp)
sourceCpp("~/Documents/PELT.cpp")
set.seed(1)
cps <- PELTgamma(y)
plot.cps <- function(y,cps){
  plot(y,type="l")
  cps <- c(1,cps,length(y))
  means <-sapply(1:(length(cps)-1),function(i){par = YZmethod(y[cps[i]:cps[i+1]]);par[1]*par[2]})
  for(i in 1:(length(cps)-1))segments(cps[i],means[i],cps[i+1],means[i],col="red",lwd=2)
}
set.seed(1)
y0 <- rgamma(50,2,0.1)
cps <-PELTgamma(y0)
plot.cps(y0,cps)
y1 <- c(rgamma(20,2,0.1),rgamma(20,2,0.5))
cps <-PELTgamma(y1)
plot.cps(y,cps)
y2 <- c(rgamma(25,1,1),rgamma(50,0.1,1),rgamma(25,1,0.5))
cps <-PELTgamma(y2)
plot.cps(y2,cps)

変化点0個のとき:
f:id:abrahamcow:20190606003343p:plain

変化点1個のとき:
f:id:abrahamcow:20190606003354p:plain

変化点2個のとき:
f:id:abrahamcow:20190606003405p:plain

なかなかうまく動きそう。

非負の連続量の時系列データ、たとえば風速とか車の時速のデータに使えたりするかもしれません。