stanでゼロ過剰ポアソン&指数・双極関数
(W)AICと(W)BICによるモデル比較

author: Unadon (見習い飯炊き兵) 動作環境:Mac OS Sierra 10.12.1; R version3.3.1; rstan 2.10.1

ゼロ過剰ポアソン、指数・双極関数フィット、モデル比較

ExpHypSim.jpeg InformationCriteria.jpeg

Figure. 上左: ShinyAppによる指数関数と双極関数のシミュレーション、
上右: (W)AICと(W)BICによるモデル比較、
パネル下: WinBugsの検索人気度の推移(ゼロ過剰ポアソンを用いた指数関数と双極関数のフィッティング)

はじめに

時系列データの各時点のデータに対して、ゼロ過剰ポアソン分布を適用。

ポアソン分布の部分の平均パラメータlambdaに対して、指数関数と双極関数をフィットさせてモデル比較しました。

これが妥当かどうかはさておき(時系列データに線を当てはめるっていいの?各時点の値は独立じゃないのに?)、

「こういうことができるんだ」という参考になればとおもいます。

テーマはGibbsSamplingのベイズソフトウェア「Winbugsの未来」。

GoogleTrendsの検索人気度データを使用します。

本稿で実施すること

  1. GoogleTrendsからWinbugsの検索人気度を取得

  2. ローデータの可視化(Line, Point, Histogram)

  3. Zero-Inflated Poissonと指数関数のフィット

  4. Zero-Inflated Poissonと双極関数のフィット

  5. モデル比較(AIC, BIC, WAIC, WBIC)

なお、コードが膨大なので、一部省略しながらご紹介させていただきます。

1. GoogleTrendsからWinbugsの検索人気度を取得

早速データ取得と可視化まで行きます。googleアカウントのメールアドレスとパスワードを入力し、実行してください。

############################################################
#  Zero-Inflated Poisson Exponential & Hyperbolic
#    Google Trends Analysis "Winbugs Future"
#       Author: Mr.Unadon
############################################################
invisible({rm(list=ls());gc();gc()})
#グーグルアカウントの設定
usr <- "" #mail adress
psw <- "" # PassWord

#キーワード
keyWords<-c("winbugs")

#データ取得の開始日と終了日(3ヶ月以上だと週1のサンプリング)
startDate<-c("2004-11-01")
endDate<-c("2016-11-01")

#read package
library(gtrendsR) # for Scraping Google Trends
library(ggplot2) # Graph
library(plotly)# Graph
library(ggraptR)# Graphic GUI
library(stringi) # Change Encoding
library(stringr) # Change Encoding
library(loo)  # WAIC and Leave-One-Out
library(rstan)
library(formattable)

#Stan Parallel processing 
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

#############################################################
#1. グーグルトレンドからデータを取得(Winbugsの検索人気度)
#############################################################

#login
gconnect(usr, psw)

#Scraping with gtrends
GetTrend <- gtrends(query = stri_encode(keyWords, "", "UTF-8"),  # Set UTF-8 or shift-jis
                    geo = "JP", # get data of "Japan Regions: JP"
                    start_date = startDate, 
                    end_date = endDate
)

2.ローデータの可視化(Line, Point, Histogram)

データの取得ができました。

geom_line,goem_point, geom_histogramの3種類のグラフを出力します。

plotlyを使ってみましょう。ggplot2で作成した図を、ggplotly()で出力するだけで、インタラクティブなグラフに早変わりです。

# Geometric Line Plot
Dat<-as.data.frame(GetTrend$trend)# ここでデータフレームにする
gline<-ggplot(Dat, aes(x = start, y = hits)) +
  geom_line(size = 1.0, colour = "turquoise4") +theme_set(theme_light(base_family="HiraKakuProN-W3"))+
  labs(colour = "検索クエリ", x = "日付", y = "検索人気度", title = GetTrend$meta) 
ggplotly(gline)

# Geometric scatter Plot
Dat<-as.data.frame(GetTrend$trend)# ここでデータフレームにする
g<-ggplot(Dat, aes(x = start, y = hits)) +
  geom_point(size = 1.0, colour = "turquoise4") +theme_set(theme_light(base_family="HiraKakuProN-W3"))+
  labs(colour = "検索クエリ", x = "日付", y = "検索人気度", title = GetTrend$meta) 
ggplotly(g)


# Geometric Histogram
ghist<-ggplot(Dat, aes(x = hits)) +
  geom_histogram(alpha = 0.5, colour = "turquoise4",fill="turquoise4") +theme_set(theme_light(base_family="HiraKakuProN-W3"))+
  labs(colour = "検索クエリ", x = "日付", y = "検索人気度", title = GetTrend$meta) 
ggplotly(ghist)

lineでプロットするとなんとなく右肩下がりなトレンドが確認できます。

Pointで描くと、ゼロが多く、また散らばり具合もわりと広範であることが。ヒストグラムにすると、ゼロの多さがよくわかります。

今回、検索した単語は”winbugs”でした。ベイズ統計の近傍にいる人くらいしか検索しないワードでしょう。 ゼロがおおくて当然といえば当然です。

ですので今回は、”検索があった時の人気度の推移”をモデルで表現してみたいと思います。

つまるところ、散布図の”ゼロじゃない上半分のトレンド”を数式で表現したい。

という願いを叶えるための試みです。 今回使用する数式は、指数関数(パラメータはgamma)と双極関数(パラメータk, 今回はHypGammaと表記)。

指数関数は”なだらかに”しかし”一気にゼロまで降下”していくトレンド。双極関数は”急速に”しかし”なかなかゼロには近づかない”トレンドを表現します。

ExpHypSim.jpeg

こちらで少し、遊んでみてください。

ShinyAppを用いた指数・双極関数のシミュレーション

3.Zero-Inflated Poissonと指数関数のフィット

さて、GoogleTrendのデータをベイズしていきます。 ここで、指数関数を”ExtinctionModel”と呼称することにします。指数関数はゼロに極限しやすいためです。

GoogleTrendsのデータは相対性を表した値なので、値が0でも”検索がゼロ”なわけではありません。 したがって、本来的にExtinctionという単語は不適切なところではあります。すみません。でも変更がとても大変なのでこのまま行きます

推定の実行

#######################################################################
#2. Exponential "Extinction model" fitting  using Zero-Inflated Poisson
#   ゼロ過剰ポアソンのポアソンの部分に指数関数をフィット
#######################################################################

#data to be passed on stan
rate<-(Dat$hits)
N<-length(rate)
date<-c(1:N)
date<-date
dataStan=list(date=date,rate=rate,N=N)

#parameters kick back from stan
parameters=c("gamma","lambda","theta","log_lik")


#Zero-Inflated Poisson Exponential "Extinction model" Fitting-------------------------
  #       model<-stan_model("ZIPExp.stan")
  #       saveRDS(model,"ZIPExp.rds")
#Stan Fit
ZIPExp <-readRDS("ZIPExp.rds") 
fitZP<-sampling(ZIPExp,pars=parameters,data=dataStan,iter=1000,warmup = 500,chain=3)

.stanファイルをC++言語にコンパイルした後、読み込んでMCMCを行っています。

方法は、stanモデルを書いた後(model<-““という部分を除いたRファイル、以下のコードのように)、

  1. そのファイルを.stanの拡張子で保存

  2. model<-stan_model(“ZIPExp.stan”)で先にコンパイル

  3. saveRDS(model,”ZIPExp.rds”)で拡張子rdsにして保存

  4. ZIPExp <-readRDS(“ZIPExp.rds”) で読み込み

  5. sampling(ZIPExp,pars=parameters・・・)でサンプリングを実行(stan()ではなくsampling()にしておく)

.stanファイルの中身はこちらです。ゼロ過剰ポアソン分布のポアソンの部分。その平均パラメータlambdaを定義する部分で、指数関数を使っています。

lambda[i] = 100*(gamma^date[i])

の部分です。100がかかっているのは、データのスケール0から100だからです。

.stanファイル内のStanモデル

data {
  int N;
  int date[N];
  int <lower=0,upper=100>rate[N];
}

parameters{
  real <lower=0,upper=1>gamma;
  real <lower=0, upper=1>theta;
}

transformed parameters {
  real <lower=0, upper=100>lambda[N];
  for(i in 1:N){
    lambda[i] = 100*(gamma^date[i]);
  }
}

model {
  for(i in 1:N){
    if (rate[i] == 0) {
      target+=log(1- theta)+theta*exp(-lambda[i]);
    } else {
      target+=bernoulli_lpmf(1|theta)+ poisson_lpmf(rate[i]|lambda[i]);
    }
  }
}

generated quantities {
  vector [N]log_lik;
  for (i in 1:N){
    log_lik[i] = (rate[i] > 0) ?
    bernoulli_lpmf(1 | theta) + poisson_lpmf(rate[i] | lambda[i]) :
      log_sum_exp(bernoulli_lpmf(0 | theta),
                  poisson_lpmf(0 | lambda[i]));
  }
}

推定結果を表で確認しておきます。この表の出し方については、前記事を御覧ください。

ZIP & Exponential Model: Summary

続いて、ローデータ(散布図)にフィットした数式を重ね書きしていきます(コードは最下部)。

ゼロの多いデータですが、ゼロではない部分のトレンドを上手く表現できている感はありますね。ゼロに到達しそうな勢いです。

4.Zero-Inflated Poissonと双極関数のフィット

指数関数を”ExtinctionModel”としたので、双極関数を使ったモデルは”Survive Model”としましょう。双極関数が、なかなかゼロに到達しないトレンドを持つためです。

.stanファイルのモデルは、指数関数の部分と以下だけが異なります。

lambda[i] = 100/(1+date[i]*gamma)

推定の実行

#####################################################################
#3. Hyperbolic  "Survival model" fitting using Zero-Inflated Poisson
####################################################################

#data to be passed on stan
rate<-(Dat$hits)
N<-length(rate)
date<-c(1:N)
date<-date
dataStan=list(date=date,rate=rate,N=N)

# parameters kick back from stan
parameters=c("gamma","lambda","theta","log_lik")

#zero-inflated poisson & Hyperbolic "Survival model" Fitting-------------------------
    #   model<-stan_model("ZIPHyp.stan")
    #   saveRDS(model,"ZIPHyp.rds")
ZIPHyp <-readRDS("ZIPHyp.rds") 
fitZP2<-sampling(ZIPHyp,data=dataStan,pars=parameters,iter=1000,warmup = 500,chain=3)
Summary<-summary(fitZP2)$summary[,c(1,3,4,8,9,10)]

結果の表をまとめておきます。

ZIP & Hyperbolic Model: Summary

指数関数を使ったフィッティングの結果に、さらに双極関数を上書きします(コードは最下部)

うーん、どっちが良いのでしょうか。見た目には判断できない感じがします。

では、情報量基準を使って比較してみましょう。メインはWAICで行いますが、AICとBIC、ついでにWBICも算出してみます。

5.モデル比較(AIC, BIC, WAIC, WBIC)

今回使用する情報量規準。AICやBICは耳にされたことが有る方も多いと思います。

AIC系は、”予測の良さ”を表現するもので、BIC系は”真のモデルとの距離”を表現するものです。

W(Widely)がイニシャルにつくと、事後分布が正規分布で近似できない場合であっても使えるものになります。

今回はゼロ過剰ポアソンなるものを使用していますので、本当はAICとBICは不適当です。

また、今回は検索人気度のトレンドを推定して予測すること(未来の時系列データの予測ではありません。既に手元にあるデータの時点でもう一度データが仮に取れた時のデータに対する予測です。週に一回のデータの間引き方が変わっても予測が良いモデルの選択)を念頭においたもので、”検索人気度を生成した真のモデル”を突き止めたいわけではありませんので、WBICというのもちょっと違います。

ですが一応、勉強とおもって全部算出してみることにしました。

WAICとWBICの算出方法はこちらを参考にさせていただきました。

なお、WBICは事後分布の逆温度なるものを1/log(N)に固定して推定を実行する必要があるため、stanモデルを別で書く必要があります。

そのコードを掲載しておきます。

modelブロックの

target+=1/log(N)log(1- theta)+thetaexp(-lambda[i]);

の部分だけが異なります。

data {
  int N;
  int date[N];
  int <lower=0,upper=100>rate[N];
}

parameters{
  real <lower=0,upper=1>gamma;
  real <lower=0, upper=1>theta;
  
}

transformed parameters {
  real <lower=0, upper=100>lambda[N];
  for(i in 1:N){
    lambda[i] = 100*(gamma^date[i]);
  }
}

model {
  for(i in 1:N){
    if (rate[i] == 0) {
      target+=1/log(N)*log(1- theta)+theta*exp(-lambda[i]);
    } else {
      target+=1/log(N)*bernoulli_lpmf(1|theta)+ poisson_lpmf(rate[i]|lambda[i]);
    }
  }
}

generated quantities {
  vector [N]log_lik;
  for (i in 1:N){
    log_lik[i] = (rate[i] > 0) ?
    bernoulli_lpmf(1 | theta) + poisson_lpmf(rate[i] | lambda[i]) :
      log_sum_exp(bernoulli_lpmf(0 | theta),
                  poisson_lpmf(0 | lambda[i]));
  }
}

さて、指数関数を使ったモデルと双極関数を使ったモデルで、 AIC,BIC,WAIC,WBICを求めました。

結果はパッケージ{ggraptR}の、関数ggraptR()でサッと書きました。

InformationCriteria.jpeg

いずれの指標も、Extinction Model、すなわち指数関数をトレンドにフィットさせたものを支持しました。

説明不足な点が多々ありましたが、紙面の都合上ご容赦ください。

最後に、今回使用したコードをまとめて掲載しておきます。

それでは、

Enjoy

############################################################
#  Zero-Inflated Poisson Exponential & Hyperbolic
#    Google Trends Analysis "Winbugs Future"
#       Author: Mr.Unadon
############################################################
invisible({rm(list=ls());gc();gc()})

#グーグルアカウントの設定
usr <- "" #mail adress
psw <- "" # PassWord

#キーワード
keyWords<-c("winbugs")

#データ取得の開始日と終了日(3ヶ月以上だと週1のサンプリング)
startDate<-c("2004-11-01")
endDate<-c("2016-11-01")

#read package
library(gtrendsR) # for Scraping Google Trends
library(ggplot2) # Graph
library(plotly)# Graph
library(ggraptR)# Graphic GUI
library(stringi) # Change Encoding
library(stringr) # Change Encoding
library(loo)  # WAIC and Leave-One-Out
library(rstan)
library(formattable)

#Stan Parallel processing 
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())



#############################################################
#1. グーグルトレンドからデータを取得(Winbugsの検索人気度)
#############################################################

#login
gconnect(usr, psw)

#Scraping with gtrends
GetTrend <- gtrends(query = stri_encode(keyWords, "", "UTF-8"),  # Set UTF-8 or shift-jis
                    geo = "JP", # get data of "Japan Regions: JP"
                    start_date = startDate, 
                    end_date = endDate
)

#1.1 Geometoric point plot with ggplot2-----------------------------------------

# Geometric Line Plot
Dat<-as.data.frame(GetTrend$trend)
gline<-ggplot(Dat, aes(x = start, y = hits)) +
  geom_line(size = 1.0, colour = "turquoise4") +theme_set(theme_light(base_family="HiraKakuProN-W3"))+
  labs(colour = "検索クエリ", x = "日付", y = "検索人気度", title = GetTrend$meta) 
ggplotly(gline)

# Geometric scatter Plot
Dat<-as.data.frame(GetTrend$trend)
g<-ggplot(Dat, aes(x = start, y = hits)) +
  geom_point(size = 1.0, colour = "turquoise4") +theme_set(theme_light(base_family="HiraKakuProN-W3"))+
  labs(colour = "検索クエリ", x = "日付", y = "検索人気度", title = GetTrend$meta) 
ggplotly(g)


# Geometric Histogram
ghist<-ggplot(Dat, aes(x = hits)) +
  geom_histogram(alpha = 0.5, colour = "turquoise4",fill="turquoise4") +theme_set(theme_light(base_family="HiraKakuProN-W3"))+
  labs(colour = "検索クエリ", x = "日付", y = "検索人気度", title = GetTrend$meta) 
ggplotly(ghist)




#######################################################################
#2. Exponential "Extinction model" fitting  using Zero-Inflated Poisson
#   ゼロ過剰ポアソンのポアソンの部分に指数関数をフィット
#######################################################################

#data to be passed on stan
rate<-(Dat$hits)
N<-length(rate)
date<-c(1:N)
date<-date
dataStan=list(date=date,rate=rate,N=N)

#parameters kick back from stan
parameters=c("gamma","lambda","theta","log_lik")


#Zero-Inflated Poisson Exponential "Extinction model" Fitting-------------------------
  #       model<-stan_model("ZIPExp.stan")
  #       saveRDS(model,"ZIPExp.rds")
#Stan Fit
ZIPExp <-readRDS("ZIPExp.rds") 
fitZP<-sampling(ZIPExp,pars=parameters,data=dataStan,iter=1000,warmup = 500,chain=3)


#formattableで結果を整理-----------------------------
Summary<-summary(fitZP)$summary[,c(1,3,4,8,9,10)]

#結果をデータフレームに
df<-data.frame(Summary)
#RhatのたまにあるNaNは数値認識されないのでゼロで置き換え
df$Rhat[is.nan(df$Rhat)]<-0
#小数点第3ケタまでにする
for(i in 1 : length(df)){df[,i]<-digits(df[,i],3)}
#対数尤度の項lp___を取り除く
df<-df[-length(df$mean),]

formattable::formattable(df, list(
    mean = formattable::formatter("span",style = x ~ formattable::style(color=ifelse(x <0, "steelblue","tomato"))),
    sd = formattable::color_bar("orange"),
    X2.5. = formattable::formatter("span",style = x ~ ifelse(x <0, formattable::style(color = "steelblue", font.weight = "bold"),NA)),
    X97.5. = formattable::formatter("span",style = x ~ ifelse(x <0, formattable::style(color = "steelblue", font.weight = "bold"),NA)),
    n_eff = formattable::formatter("span",style = x ~ formattable::style(color = ifelse(rank(-x) <= c(length(df$n_eff)-1), "gray", "tomato")),
                      x ~ sprintf("%.2f (rank: %g)", x, rank(-x))),
    Rhat = formattable::formatter("span", x ~digits(x,2),
                     style = x ~ formattable::style(color = ifelse(x>=1.1, "tomato", "green")),
                     x ~ icontext(ifelse(x<1.1, "ok", "remove"), ifelse(x<1.1, "Yes", "No")))
))
#-------------------------------------------------------------------------------------


#MCMC samples Extracting
extZP<-rstan::extract(fitZP)
gammaZP<-extZP$gamma

#MAP estimate
gammaMAPZP  = c()
gammaMAPZP<- density(gammaZP)$x[which(density(gammaZP)$y == max(density(gammaZP)$y))]
gammaMAPZP
LL<-extZP$log_lik

#for Model Prediction Line plot
predZP<-c(rep(0,N))
for(i in 1:N){
  predZP[i]<-gammaMAPZP^(i)
}

#Rescaling
predZP<-predZP*100


#MAP (as MLE) Log_likelihood for AIC and BIC
LLmle  = c()
for ( i in 1:length(LL[1,])){
  LLmle[i]<- density(LL[,i])$x[which(density(LL[,i])$y == max(density(LL[,i])$y))]
}


#2.1 Information Criteria: AIC, BIC, WAIC and WBIC-------------------------------------

#AIC
  AIC<- -2*sum(LLmle)+2*2
  AIC

#BIC
  BIC<- -2*sum(LLmle)+log(length(LL[1,]))*2
  BIC

#WAIC
  library(loo)
  LL <- extract_log_lik(fitZP)
  WAIC <- waic(LL)
  print(WAIC, digits = 4)

#WBIC
  ZIPExpWBIC <-readRDS("ZIPExpWBIC.rds") 
    #       model<-stan_model("ZIPExpWBIC.stan")
    #       saveRDS(model,"ZIPExpWBIC.rds")

parameters=c("gamma","lambda","theta","log_lik")
fitZPwbic<-sampling(ZIPExpWBIC,pars=parameters,data=dataStan,iter=1000,warmup = 500,chain=3)


#formattableで結果を整理---------------------------
Summary<-summary(fitZPwbic)$summary[,c(1,3,4,8,9,10)]

#結果をデータフレームに
df<-data.frame(Summary)
#RhatのたまにあるNaNは数値認識されないのでゼロで置き換え
df$Rhat[is.nan(df$Rhat)]<-0
#小数点第3ケタまでにする
for(i in 1 : length(df)){df[,i]<-digits(df[,i],3)}
#対数尤度の項lp___を取り除く
df<-df[-length(df$mean),]

formattable::formattable(df, list(
    mean = formattable::formatter("span",style = x ~ formattable::style(color=ifelse(x <0, "steelblue","tomato"))),
    sd = formattable::color_bar("orange"),
    X2.5. = formattable::formatter("span",style = x ~ ifelse(x <0, formattable::style(color = "steelblue", font.weight = "bold"),NA)),
    X97.5. = formattable::formatter("span",style = x ~ ifelse(x <0, formattable::style(color = "steelblue", font.weight = "bold"),NA)),
    n_eff = formattable::formatter("span",style = x ~ formattable::style(color = ifelse(rank(-x) <= c(length(df$n_eff)-1), "gray", "tomato")),
                                   x ~ sprintf("%.2f (rank: %g)", x, rank(-x))),
    Rhat = formattable::formatter("span", x ~digits(x,2),
                                  style = x ~ formattable::style(color = ifelse(x>=1.1, "tomato", "green")),
                                  x ~ icontext(ifelse(x<1.1, "ok", "remove"), ifelse(x<1.1, "Yes", "No")))
))


LLwbic <- extract_log_lik(fitZPwbic)
WBIC <- - mean(rowSums(LLwbic))
WBIC


#put Fitting results into data.frame
modelName<-rep(c("Extinction"),4)
ICs<-c("AIC","BIC","WAIC","WBIC")
value<-c(AIC,BIC,WAIC$waic,WBIC)

dfZIPExp<-data.frame(list(Model=modelName,ICs=ICs,value=value))


#ggplot
fitDat<-data.frame(list(predZP=predZP,date=date,rate=rate,x=Dat$start))

g2<-ggplot(fitDat,aes(x=x,y=rate))+
  theme_set(theme_bw(base_family="HiraKakuProN-W3"))+
  geom_point(size = 1.0, colour = "turquoise4") +
  ylim(0,100)+
  geom_line(aes(x,predZP),colour="darkblue",size=1)+
  labs(colour = "検索クエリ", x = "日付", y = "検索人気度", title = GetTrend$meta) 
ggplotly(g2)


#####################################################################
#3. Hyperbolic  "Survival model" fitting using Zero-Inflated Poisson
####################################################################

#data to be passed on stan
rate<-(Dat$hits)
N<-length(rate)
date<-c(1:N)
date<-date
dataStan=list(date=date,rate=rate,N=N)

# parameters kick back from stan
parameters=c("gamma","lambda","theta","log_lik")

#zero-inflated poisson & Hyperbolic "Survival model" Fitting-------------------------
    #   model<-stan_model("ZIPHyp.stan")
    #   saveRDS(model,"ZIPHyp.rds")
ZIPHyp <-readRDS("ZIPHyp.rds") 
fitZP2<-sampling(ZIPHyp,data=dataStan,pars=parameters,iter=1000,warmup = 500,chain=3)
Summary<-summary(fitZP2)$summary[,c(1,3,4,8,9,10)]


#formattableで結果を整理------------------------
#結果をデータフレームに
df<-data.frame(Summary)
#RhatのたまにあるNaNは数値認識されないのでゼロで置き換え
df$Rhat[is.nan(df$Rhat)]<-0
#小数点第3ケタまでにする
for(i in 1 : length(df)){df[,i]<-digits(df[,i],3)}
#対数尤度の項lp___を取り除く
df<-df[-length(df$mean),]

formattable::formattable(df, list(
    mean = formattable::formatter("span",style = x ~ formattable::style(color=ifelse(x <0, "steelblue","tomato"))),
    sd = formattable::color_bar("orange"),
    X2.5. = formattable::formatter("span",style = x ~ ifelse(x <0, formattable::style(color = "steelblue", font.weight = "bold"),NA)),
    X97.5. = formattable::formatter("span",style = x ~ ifelse(x <0, formattable::style(color = "steelblue", font.weight = "bold"),NA)),
    n_eff = formattable::formatter("span",style = x ~ formattable::style(color = ifelse(rank(-x) <= c(length(df$n_eff)-1), "gray", "tomato")),
                                   x ~ sprintf("%.2f (rank: %g)", x, rank(-x))),
    Rhat = formattable::formatter("span", x ~digits(x,2),
                                  style = x ~ formattable::style(color = ifelse(x>=1.1, "tomato", "green")),
                                  x ~ icontext(ifelse(x<1.1, "ok", "remove"), ifelse(x<1.1, "Yes", "No")))
))


#-------------------------------------------------------------------------------------


# MCMC Sampling Extraction
extZP2<-rstan::extract(fitZP2)
LL2<-extZP2$log_lik
gammaZP2<-extZP2$gamma

#MapEstimate
gammaMAPZP2  = c()
gammaMAPZP2<- density(gammaZP2)$x[which(density(gammaZP2)$y == max(density(gammaZP2)$y))]
gammaMAPZP2


#for Prediction Line plot
predZP2<-c(rep(0,N))
for(i in 1:N){
  predZP2[i]<-1/(1+(gammaMAPZP2*(i)))
}
#rescaling
predZP2<-predZP2*100


# MAP (as MLE) Log_likelihood 
LLmle2  = c()
for ( i in 1:length(LL2[1,])){
  LLmle2[i]<- density(LL2[,i])$x[which(density(LL2[,i])$y == max(density(LL2[,i])$y))]
}


# Results as dataframe
fDat<-data.frame(list(predZP=predZP,predZP2=predZP2,date=date,rate=rate,x=Dat$start))

#plot
g3<-ggplot(fDat,aes(x=x,y=rate))+
  theme_set(theme_bw(base_family="HiraKakuProN-W3"))+
  geom_point(size = 1.0, colour = "turquoise4") +
  ylim(0,100)+
  #geom_line(aes(x,predZP),colour="darkblue",size=1)+   # Extinction model line
  geom_line(aes(x,predZP2),colour="springgreen4",size=1)+    # Survival Model
  labs(colour = "検索クエリ", x = "日付", y = "検索人気度", title = GetTrend$meta) 
ggplotly(g3)

#   3.1 Information Criteria: AIC, BIC, WAIC and WBIC---------------------------------

library(loo)
LL2 <- extract_log_lik(fitZP2)
WAIC2 <- waic(LL2)
print(WAIC2, digits = 4)

AIC2<- -2*sum(LLmle2)+2*2
AIC2

BIC2<- -2*sum(LLmle2)+log(length(LL2[1,]))*2
BIC2

#WBIC
ZIPHypWBIC <-readRDS("ZIPHypWBIC.rds") 
  #       model<-stan_model("ZIPHypWBIC.stan")
  #       saveRDS(model,"ZIPHypWBIC.rds")
parameters=c("gamma","lambda","theta","log_lik")
fitZP2wbic<-sampling(ZIPHypWBIC,pars=parameters,data=dataStan,iter=1000,warmup = 500,chain=3)
Summary<-summary(fitZPwbic)$summary[,c(1,3,4,8,9,10)]

#formattableで結果を整理----------------------
#結果をデータフレームに
df<-data.frame(Summary)
#RhatのたまにあるNaNは数値認識されないのでゼロで置き換え
df$Rhat[is.nan(df$Rhat)]<-0
#小数点第3ケタまでにする
for(i in 1 : length(df)){df[,i]<-digits(df[,i],3)}
#対数尤度の項lp___を取り除く
df<-df[-length(df$mean),]

formattable::formattable(df, list(
    mean = formattable::formatter("span",style = x ~ formattable::style(color=ifelse(x <0, "steelblue","tomato"))),
    sd = formattable::color_bar("orange"),
    X2.5. = formattable::formatter("span",style = x ~ ifelse(x <0, formattable::style(color = "steelblue", font.weight = "bold"),NA)),
    X97.5. = formattable::formatter("span",style = x ~ ifelse(x <0, formattable::style(color = "steelblue", font.weight = "bold"),NA)),
    n_eff = formattable::formatter("span",style = x ~ formattable::style(color = ifelse(rank(-x) <= c(length(df$n_eff)-1), "gray", "tomato")),
                                   x ~ sprintf("%.2f (rank: %g)", x, rank(-x))),
    Rhat = formattable::formatter("span", x ~digits(x,2),
                                  style = x ~ formattable::style(color = ifelse(x>=1.1, "tomato", "green")),
                                  x ~ icontext(ifelse(x<1.1, "ok", "remove"), ifelse(x<1.1, "Yes", "No")))
))

#WBIC
LL2wbic <- extract_log_lik(fitZP2wbic)
WBIC2 <- - mean(rowSums(LL2wbic))
WBIC2


#put Fitting results into data.frame
modelName<-rep(c("Survival"),4)
ICs<-c("AIC","BIC","WAIC","WBIC")
value<-c(AIC2,BIC2,WAIC2$waic,WBIC2)

dfZIPHyp<-data.frame(list(Model=modelName,ICs=ICs,value=value))


##########################################################################
#4. Fitting Summary Graphics
##########################################################################

#plot
g3<-ggplot(fDat,aes(x=x,y=rate))+
    theme_set(theme_bw(base_family="HiraKakuProN-W3"))+
    geom_point(size = 1.0, colour = "turquoise4") +
    ylim(0,100)+
    geom_line(aes(x,predZP),colour="darkblue",size=1)+   # Extinction model line
    geom_line(aes(x,predZP2),colour="springgreen4",size=1)+    # Survival Model
    labs(colour = "検索クエリ", x = "日付", y = "検索人気度", title = GetTrend$meta) 
ggplotly(g3)


# Summary Information Criteria Bar plot
dfSum<-rbind(dfZIPExp,dfZIPHyp)

ggraptR()



Written on February 6, 2017