二値データのベイズ時系列解析(続編)


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



はじめに


本稿は、ベータ分布を用いた時系列解析例に関する前記事の続編です。

「ベータ分布使うより、ロジットで正規分布をたたんだほうが良いのでは」

<なるほど!>といったお話を頂きましたので、

頂いたアドバイスを元に再度整理をしてみます。

以下は、別記事で同じモデルを取り上げた際に頂いたコメントです。ありがとうございます。

モデリングの発展例なども紹介されています。御覧くださいませ。

ツイートの引用を加筆いたしました(2017/08/15 18:51)

頂いたコメントやアドバイスをもとに、結論から整理すると、次の点がポイントとなりそうです。

  1. ベータ分布のモデルより、今回のモデルの方が分散の解釈がストレート。
  2. 今回のロジットのモデルは、幾つか軽い制約を置くと収束しやすい。
  3. 今回のロジットのモデルの方が、解析時間が断然はやい。
  4. 二値データのトレンドを可視化できるのは、いいですよね

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

さて、時系列の二値データが得られる場合というのは多々あることと思います。

連続量だけどデータがどうもよくないので「購入 vs 未購入」に変数変換したという場合などです。

ということで、前記事と同じデータを使用します。

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

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

下記データから、事象の生起確率の時系列トレンドを推定してみましょう。

invisible({rm(list=ls());gc();gc()})
#################################################################
# 二値データの時系列解析: Gaussian + inv_logitの動的モデル
#   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")),
                  TimePoint=seq(1:length(Dat)))

データの可視化


ggplot2で確認です。

前記事にも同じことを書いていますが、 下のようなデータが出てきて、「どういうトレンドなの?」と聞かれると僕は嬉しく思います。

普通の時系列データとちがって、分析しないとトレンドが見えないからです。やりがいがありますよね。

#ggplot
colour1=c("#00CD66", "#8968CD")
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=2,alpha=0.5)+ 
    guides(colour=FALSE,fill=FALSE)+
    theme(axis.title = element_text(size = 14), 
      axis.text = element_text(size = 10), 
      plot.title = element_text(size = 18, 
                                hjust = 0.5)) +
    labs(title = "二値データの時系列解析2: CV確率推移の推定", y = "CV確率")+
    theme(axis.title = element_text(family = "serif"), 
          axis.text = element_text(family = "serif"), 
          plot.title = element_text(family = "serif")) +
    labs(title = "二値データの時系列解析: 生起確率のトレンド推定(Gaussian_logit)", y = "probability")

BinaryDatGaussianLogitDynamic.png
Figure. サンプル二値データ

Gaussian Logit変換の動的モデル


どういう名前で呼べばいいのかわかりませんが、ここではGaussianLogitモデルと呼ぶことにします。

基本の発想は、最もシンプルな状態空間モデルです。

観測の世界と真の状態の世界を2つに分け、

観測誤差と状態誤差(遷移の幅)を考えるモデルです。

このとき、状態(システム)モデルはマイナス∞から+∞の正規分布の世界で推移していくと仮定します。

観測モデルはベルヌイ分布。その生起確率は、逆ロジットをかませた正規分布で表現されます。

状態の推移は、「一時点前に依存する」「一時点前の状態平均からσの正規分布で移動する」というローカルレベルモデルの発想を採用しました。


Stanコードと実行


では、推定に移ります。まずはStanファイルのコードから。

data {
    int t;
    int Dat[t]; 
}

parameters {
    vector <lower = -5, upper = 5>[t] mu;
    real<lower=0> sigma;
}

transformed parameters {
    vector<lower=0, upper=1>[t] theta;
    theta = inv_logit(mu);
}

model {
    sigma ~ normal(0, 1);
    for (i in 2:t){
        mu[i] ~ normal(mu[i-1], sigma);
    }
    Dat ~ bernoulli(theta);
}

短くていいですね。ここでは、状態モデルの正規分布に上限と下限をあたえて制約を付けています。 こうしないと、sigmaパラメータのNeffが怪しくなってきます。他にもいい方法があるかも。未検証です。

続いて、R側のコードです。 収束は事前に確認しております。

コードを回して、推定を実行します。

#data to be passed on stan
t=length(Dat)
datastan=list(t=t,Dat=Dat)

#
# stan model
#
#model<-stan_model("GaussianLogitDynamic.stan")
#saveRDS(model,"GaussianLogitDynamic.rds")
GLD<-readRDS("GaussianLogitDynamic.rds")

#
# stan fit
#
fit = sampling(GLD, data = datastan, chains = 4, iter = 4000, warmup = 2000, thin=4, seed=1234)

#Gelman-Rubin 
launch_shinystan(fit)

(  ´Д`)ハァハァ


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


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

#Sample Extract and Summarize
df_result<-rstan::extract(fit)$theta%>%
    as.data.frame()%>%
    gather()%>%
    mutate(TimePoint=rep(1:length(Dat),each=2000))%>%
    group_by(TimePoint)%>%
    summarize(EAP=mean(value),
              up=quantile(value,0.975),lo=quantile(value,0.025),
              up50=quantile(value,0.75),lo50=quantile(value,0.25))%>%
    ungroup()%>%
    left_join(df_plot,by="TimePoint")

#Result Plot-------------------------------------------------------------
ggplot(df_result,aes(x=Date))+
    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=2,alpha=0.5)+ 
    geom_hline(yintercept = 0.5,linetype=2,colour="gray20")+
    geom_ribbon(aes(ymin=lo,ymax=up),alpha=0.2)+
    geom_ribbon(aes(ymin=lo50,ymax=up50),alpha=0.5)+
    geom_line(aes(y=EAP))+
    guides(colour=FALSE,fill=FALSE) +
    theme(axis.title = element_text(size = 14), 
          axis.text = element_text(size = 10), 
          plot.title = element_text(size = 18, 
                                    hjust = 0.5)) +
    labs(title = "二値データの時系列解析2: CV確率推移の推定", y = "CV確率")+
    theme(axis.title = element_text(family = "serif"), 
    axis.text = element_text(family = "serif"), 
    plot.title = element_text(family = "serif")) +labs(title = "二値データの時系列解析: 生起確率のトレンド推定(Gaussian_Logit)", 
    y = "probability")

よいしょ!!!

GaussianLogitDynamic.png
Figure Gaussian Logitモデルの推定結果(薄い帯: 95%CI, 濃い帯: 50%CI, 実線:事後平均)

前記事(Beta分布のモデル)の結果はこちら

BetaDynamicProb.png
Figure Beta分布の動的モデル推定結果

結論は冒頭に述べたとおりです

Gaussian Logitのモデルのほうが好き。

快きアドバイスを頂きました、

@berobero11様、@ibaibabaibai様、@simizu706様に、

重ねて御礼申し上げます!

Enjoy!!

Written on August 15, 2017