二値データの時系列解析
ベータ分布のベイジアン動的モデル


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



はじめに


時系列の二値データが得られる場合というのは多々あることと思います。連続量だけどデータがどうもよくないので「購入 vs 未購入」に変数変換したという場合などです。 行動実験での時系列データも二値データの場合があるかもしれません。

今回は、二値データ生起確率の時系列推移を推定・定量化してみたいと思います。

例として、「365日のデータで、CVがあった日とCVがなかった日」のデータを想定しました。CVとは、サイトクリックや訪問行動などの目的としたい出来事を指します。

パッケージとサンプルデータの読み込み。

0と1の365個の値、2016年1月1日から12月30日の365個のデータを作成します。うるう年なので、12月30日で365です。

invisible({rm(list=ls());gc();gc()})
#################################################################
# Binaryデータの時系列解析: ベータ分布の動的モデル
#   Author: Mr.Unadon
#################################################################

#packages
library(tidyverse) 
library(rstan) # HMC
library(shinystan)  #収束診断で使用
library(scales) #ggplotで時間をx軸に表示する際に使用

#stan parallel
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())


#Sample Data: 365days CV or nonCV
Dat<-c(0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,
       1,0,0,0,0,0,0,0,0,0,1,1,0,0,1,1,0,1,1,1,
       1,1,1,1,1,1,1,1,1,1,1,0,1,1,1,1,1,1,0,1,
       1,1,1,0,0,1,1,1,0,1,1,0,1,1,1,1,1,1,0,0,
       0,1,0,1,0,0,0,1,1,1,0,1,1,1,1,1,1,0,1,1,
       1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
       1,0,0,0,0,0,0,0,0,0,1,1,0,0,1,1,0,1,1,1,
       1,1,1,1,1,1,1,1,1,1,1,0,1,1,1,1,1,1,0,1,
       0,1,0,1,0,0,0,1,1,1,0,1,1,1,1,1,1,0,1,1,
       1,0,0,0,0,0,0,0,0,0,1,1,0,0,1,1,0,1,1,1,
       0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,
       1,1,1,0,0,1,1,1,0,1,1,0,1,1,1,1,1,1,0,0,
       0,1,0,1,0,0,0,1,1,1,0,1,1,1,1,1,1,0,1,1,
       1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
       1,0,0,0,0,0,0,0,0,0,1,1,0,0,1,1,0,1,1,1,
       1,1,1,1,1,1,1,1,1,1,1,0,1,1,1,1,1,1,0,1,
       0,1,0,1,0,0,0,1,1,1,0,1,1,1,1,1,1,0,1,1,
       1,0,0,0,0,0,0,0,0,0,1,1,0,0,1,1,0,1,1,1,
       1,1,1,1,1)

#as.data.frame and mutate "Date"
df_plot<-data.frame(binaryCV=Dat)%>%
    dplyr::mutate(Date=as.POSIXct(seq(as.Date("2016-01-01"), as.Date("2016-12-30"), by="days")))


データの可視化


ggplot2でデータを見てみましょう。 こういうのが出てきて、「どういうトレンドなの?」と聞かれると僕は嬉しく思います。 普通の時系列データとちがって、分析しないとトレンドが見えないからです。やりがいがあります。

BinaryDat.png

コードはこちらです。

#ggplot
colour1=c( "#EEC900","#8B0A50")
ggplot(df_plot,aes(x=Date,y=binaryCV))+theme_light()+
    scale_color_manual(values=colour1)+
    scale_x_datetime(expand = c(0,0) ,breaks = date_breaks("30 days"))+
    theme(axis.text.x = element_text(angle = 45, hjust = 1))+
    geom_point(aes(colour=as.factor(Dat)),size=3)+ 
    guides(colour=FALSE,fill=FALSE)


ベータ分布の動的モデル


いくつかやり方があると思いますが、 今回は各時点CV生起の有無をベルヌイ分布で表現し(観測方程式に相当?)、 ベルヌイ分布の生起確率をベータ分布で表現します(状態方程式に相当?)。

「生起確率の推移は一時点前の生起確率に依存する」というローカルレベルモデルの発想を採用します。 なお、ここでの状態には、ベータ分布の平均を採用します。 そして、t時点のベータ分布の平均は、一時点前のベータ分布から出てくるものと仮定します。

「あれ、ベータ分布のパラメータってa,bで、平均と分散はどうするの?」 となるとおもいますので、ベータ分布の平均等について整理します。

ベータ分布の平均muはa/(a+b)で定義されます。 また、分散と同じ機能を持つphiはa+bで扱うことができます。(分散は別にちゃんと定義があります)

さて、muとphiをつかってベータ分布のパラメータa,bを定義しましょう。簡単です。

a = mu * phi
b = (1-mu) * phi

なお、今回推定したいのはパラメータaとbではなく、muです。 muは一時点前のベータ分布に従うということにしましたので、 状態遷移の式はこうなります。

for ( i in 1 : t){
  mu[i] ~ beta(a[i-1],b[i-1])
}

Stanコードと実行


では、推定に移ります。データとパラメータを指定して、Stanにパスします。 .stanファイルにモデルコードを書いても、そのまま推定していただいても構いません。


#data to be passed on stan
t=length(Dat)
datastan=list(t=t,Dat=Dat)
parameters=c("a","b","mu","phi")


model<-"
data {
    int t;
    int Dat[t];
}
parameters{
    vector <lower=0, upper=1> [t] mu; 
    real <lower=0> phi; 
}
transformed parameters{
    vector <lower=0> [t] a;
    vector <lower=0> [t] b;
        a = mu*phi;
        b = (1-mu)*phi;
}
model {
    phi ~ gamma(0.01,0.01);
    mu[1] ~ beta(1,1);
        for(i in 2 : t ){
            mu[i] ~ beta(a[i-1],b[i-1]);
        }
    Dat ~ bernoulli(mu);
}

"

#stan fit-------------------------------------------------------------
fit<-stan(model_code=model,data = datastan,pars=parameters,iter=4000,warmup = 1000,chain=3,thin = 10)

#from .stan file
    #fit<-stan(file="BetaDynamic.stan",data=datastan,pars=parameters)
#a priori compiling and read .rds file
    #BetaDynamic<-stan_model("BetaDynamic.stan")
    #saveRDS(BetaDynamic,"BetaDynamic.rds")
    #BetaDynamic<-readRDS("BetaDynamic.rds")
    #fit<-sampling(BetaDynamic,data=datastan,pars=parameters)
#variation bayes
    #fitVB<-vb(BetaDynamic,data=datastan,pars=parameters)


収束の確認と可視化


上手くいったかなぁ。 shinystanで結果を見てみましょう。

launch_shinystan(fit)

trace plotと事後分布をひとつ見てみます。推定自体は成立していそうです。

BetaDynamicShiny.png

 つづいて、shinystanを使ってGelman-Rubin testで収束の確認をします。 Neffが全サンプルの10%以上、MCSEが事後分布SDの10%以下、Rhatが1.1未満でOKだったと認識してます。

お、いい感じですね。全て基準を満たしています。

BetaDynamicGelman.png


二値データ生起確率の推移: 結果の可視化


最後に、MCMCサンプルを取り出して推定結果を確認しましょう。 最初にお見せしたデータに生起確率の推移を重ねて描いていきます。

いやぁ、この瞬間は楽しみですね。

CV確率の推移を上手く可視化できるかな?

0と1の間をさまよう確率の変化で、何かを語ることができるといいな。

さて、こちらがggplotコードです。

#stan results---------------------------------------------------------
#extract samples
mu<-data.frame(rstan::extract(fit)$mu)

#統計量の算出
#このfor文なんとかならんかな・・・
up<-c()
lo<-c()
up50<-c()
lo50<-c()
MAP<-c()
for(i in 1 :t){
    up[i]<-quantile(mu[,i],0.975)
    lo[i]<-quantile(mu[,i],0.025)
    up50[i]<-quantile(mu[,i],0.75)
    lo50[i]<-quantile(mu[,i],0.25)
    MAP[i]<-density(mu[,i])$x[which(density(mu[,i])$y == max(density(mu[,i])$y))]
}

#adding the results to dataframe
df_plot<-df_plot%>%dplyr::mutate(up=up,lo=lo,up50=up50,lo50=lo50,MAP=MAP)


#Result Plot-------------------------------------------------------------
ggplot(df_plot,aes(x=Date))+theme_set(theme_light(base_family = "HiraKakuProN-W3"))+
    scale_color_manual(values=colour1)+
    scale_x_datetime(expand = c(0,0) ,breaks = date_breaks("30 days"))+
    theme(axis.text.x = element_text(angle = 45, hjust = 1))+
    geom_point(aes(y=binaryCV,colour=as.factor(Dat)),size=3)+ 
    geom_hline(yintercept = 0.5,linetype=2,colour="gray20")+
    geom_ribbon(aes(ymin=lo,ymax=up),alpha=0.95)+
    geom_ribbon(aes(ymin=lo50,ymax=up50),alpha=0.99)+
    geom_line(aes(y=MAP))+
    guides(colour=FALSE,fill=FALSE) +
    theme(axis.title = element_text(size = 22), 
    axis.text = element_text(size = 10), 
    plot.title = element_text(size = 18, 
        hjust = 0.5)) +labs(title = "二値データの時系列解析: ベータ分布を用いたCV確率推移の推定", 
    y = "CV確率")

よいしょ!!!

BinaryProb.png

うぉっ、hoxo_m さんでてきた!!!

Enjoy!!

Written on May 7, 2017