廿TT

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

R: シニングによる非定常ポアソン過程のシミュレーション

ポアソン過程は適当に点を間引いてやっても(非定常)ポアソン過程になる。

この性質を使って非定常ポアソン過程をシミュレートすることができて、それをシニング(thinning)と呼ぶ。

詳しくは An Introduction to Statistical Computing: A Simulation-based Approach (Wiley Series in Computational Statistics) などに出ている。

やっていることは棄却サンプリング(R による棄却サンプリング法(おれが)入門 - 廿TT)と同じである。

シミュレーション

さて、シニングのアルゴリズムを愚直に実装するとこんな感じになるだろう。

get_NHPP <-function(iter,lambda,lambda2){
  Ts = c()
  Tast = 0
  for(i in 1:iter){
    Tast = Tast + rexp(1,lambda2)
    if(runif(1) < lambda(Tast)/lambda2){
      Ts =c(Ts,Tast)
    }
  }
  return(Ts)
}

ここで lambda2 は非定常ポアソン過程の強度関数(intensity function) lambda より大きい値になるように選ぶ。

で、R では実は for 文で回す必要なくて、

get_NHPP <-function(iter,lambda,lambda2){
    Tast = cumsum(rexp(iter,lambda2))
    Ts = ifelse(runif(iter) < lambda(Tast)/lambda2,Tast,NA)
    Ts[!is.na(Ts)]
}

と書いても同じ結果が得られる。

GO モデル

ソフトウェアのバグの発見数のモデルとして有名なものに Goel-Okumoto モデル(G-O モデル)がある。

累積強度関数の形は指数分布の分布関数を定数倍したのと同じであり、以下の微分方程式の解から出てくる。

d\Lambda(t)/dt=a(m-\Lambda(t))

バグの総数を表すパラメータが m で、a は比例定数。

つまり、バグがいっぱい残っているときはバグが発見されやすく、少なくなってくると発見されにくくなる、というメカニズムを反映している。

今回はこれをシミュレートしてみる。

#累積強度関数
GO <- function(pars){
  function(t){
    m=pars[1]
    a=pars[2]
    m*(1-exp(-a*t))
  }}

#強度関数
go <- function(pars){
  function(t){
    m=pars[1]
    a=pars[2]
    m*a*exp(-a*t)
  }}

y <-get_NHPP(iter =10000, go(c(1000,0.2)),200)

強度関数のノンパラメトリックな推定量は以下の lambdahat のようにして求める。

xv <-seq(1,ceiling(max(y)),by=1)
lambdahat <-sapply(xv,function(t)sum(t-1 < y & y<t))

詳しくは http://ebsa.ism.ac.jp/ebooks/sites/default/files/ebook/1868/pdf/vol2_ch3.pdf などを参照。

ノンパラメトリックな推定量をグラフにして、G-O モデルの強度関数を重ねてみる。

plot(xv,lambdahat)
curve(go(c(1000,0.2))(x),add=TRUE,col="firebrick",lwd=2)

f:id:abrahamcow:20170523032736p:plain


また、累積強度関数についても同じようにプロットしてみる。

plot(y,1:length(y),type="s",xlab="",ylab="")
curve(GO(c(1000,0.2))(x),add=TRUE,col="firebrick",lwd=2)

f:id:abrahamcow:20170523032732p:plain

ただしくシミュレーションできているらしいことがわかる。

最尤推定

非定常ポアソン過程の尤度関数は、In All Likelihood: Statistical Modelling and Inference Using Likelihood などに記述がある。

尤度関数は以下に比例する。

\prod_{i=1}^{n} \lambda ( t_{i}) \exp\left(-\Lambda(T)\right)

t_i はイベントの発生した時刻、T は観測期間の長さ、\lambda(x) は強度関数、\Lambda(x) は累積強度関数で、\lambda(x)積分したもの。

これを optim に突っ込めばパラメータの最尤推定ができる。

こんなふうだ。

loglik <- function(par){
  maxT <- max(y)
  lambda <- log(go(exp(par))(y))
  sum(lambda) - GO(exp(par))(maxT)
}
opt1 <-optim(c(7,0),loglik,control = list(fnscale=-1))
exp(opt1$par)

以上です。

In All Likelihood: Statistical Modelling and Inference Using Likelihood

In All Likelihood: Statistical Modelling and Inference Using Likelihood