廿TT

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

窓打ち切りデータからのワイブル再生過程のパラメータ推定(RStan版)

詳細は 窓打ち切りデータからのワイブル再生過程のパラメータの最尤推定 - 廿TT を参照してください。

パラメータを推定するだけなら最尤法の枠組みのなかで十分なんですが、RStan でもできるかな、という練習です。

#シミュレーションでデータを生成
simu1 <- function(){
  #set <- numeric(0)
  set <-vector("list",length(pics))
  for(j in 1:pics){
    x <- rweibull(5000, shape=3, scale=10)
    t2 <- cumsum(x)
    o = runif(1,100,1000)
    sta <- which.max(t2 > o)
    en <- which.max(t2 > o+win)
    Z <- x[(sta):(en)]  
    N <- length(Z)
    A <- numeric(N)
    A[1] <- 1
    D <- rep(1,N)
    D[N] <- 0
    if(N==1){Z=win}else{
      Z[1] <- t2[sta] - o
      Z[N] <- x[en] - (t2[en] - (o+win))
    }    
    set[[j]] <-data.frame(Z,D,A)
  }
#  return(set)
  return(do.call("rbind",set))
}
pics = 100
win = 10
dat <- simu1()
J <- dim(dat)[1]
status <- numeric(J)
status <- ifelse(dat$D == 0 & dat$A ==0,1,status)
status <- ifelse(dat$D == 1 & dat$A ==1,2,status)
status <- ifelse(dat$D == 0 & dat$A ==1,3,status)
dat4stan <-list(J = J, tt = dat$Z, status = status)
library(rstan)
stanmodel <- stan_model("weibull_renewal.stan")

Stan のコードはこうです。

data {
  int<lower=1> J;
  int status[J];
  real<lower=0> tt[J];
}
parameters {
  real<lower=0> m; 
  real<lower=0> eta;
}
model {
  for(j in 1:J){
    if(status[j] == 0){
      increment_log_prob(weibull_log(tt[j], m, eta));
    }
    if(status[j] == 1){
      increment_log_prob(weibull_ccdf_log(tt[j], m, eta));
    }
     if(status[j] == 2){
     increment_log_prob(weibull_ccdf_log(tt[j], m, eta) - log(eta) - lgamma(1+1/m));
    }
   if(status[j] == 3){
     increment_log_prob(
     log1m(
     gamma_cdf( pow(tt[j]/eta, m),1+1/m, 1)
      + tt[j]*(1-weibull_cdf(tt[j], m, eta))/(eta*tgamma(1+1/m))
     )
     );
  }
  }
}

 \log(1-x) を log1m(x) で表すところと、べき乗を pow で表すところがつまづきポイントでした。

> traceplot(fit)
> ex <- extract(fit)
> summary(ex$m)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  2.028   2.964   3.209   3.232   3.481   5.220 
> summary(ex$eta)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  8.307   9.537   9.817   9.842  10.130  12.180 

f:id:abrahamcow:20150901045708p:plain

基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門

基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門