12/22(日)にBUGS/Stan勉強会#2がドリコム株式会社にて催されました。そこで2つ発表をしました。そのうちの1つ「『予測にいかす統計モデリングの基本』の売上データの分析をトレースしてみた」に関する詳細＆補足＆苦労話をここで書きたいと思います。RStanというパッケージでRからStanというMCMCサンプリングソフトを使っています。
最初に発表内容のスライドは以下になります。ざっと見るにはこれで十分です。
以降ではスライドごとに簡単に補足していきます。
予測にいかす統計モデリングの基本―ベイズ統計入門から応用まで (KS理工学専門書)
- 作者:樋口 知之
- 発売日: 2011/04/07
- メディア: 単行本（ソフトカバー）
まずは元となった書籍の紹介です。時系列解析の第一人者による分かりやすく丁寧に書かれた良書です。入門的な事項から実際にデータ分析する際のノウハウまで幅広く書かれています。パーティクルフィルタ（粒子フィルタ）に知りたい人もこの本が今のところ一番と思います。僕はBUGSやStanでパラメータ（の分布）を推定してしまうので、書籍の中でよく出てくる最尤推定や最適化などはあまり深く読んではいません。
我々が統計モデリングを勉強＆実践する際のはじめの難所は手頃なデータが得られるかどうかです。手頃なデータが公開されることは僕にとっては論文よりも価値があります。今回のデータは問い合わせをしましたところ、データ提供者からデータ配布の許諾を得ていないとのことで公開は難しいとのことでした。そこで、分析結果がなるべく保持されるようにデータをゼロから作成しました。かなり苦労しました。以下のようなデータです。公開します（data_20131226.txt）。発表資料と異なり、各分解項の真の値も入れてあります。自由に使ってください（売上の単位が[千円]なのに小数点4桁以上あってごめんなさい…）。
こういったデータを渡されてまずやることは可視化です。経験的に売上は曜日や天気と関係してるのかな？とかイベントとは正の相関があるだろーとかをプロットしていく過程で色々見ておきたいことがでてきますよね。いつもおおよそ10枚ぐらい図を作ります。次のステップとしては相関係数求めてみたり、重回帰であてはまり具合を見ますかね。しかしながら重回帰はガチガチに制限された条件が前提となっています。詳しくは久保先生の緑本の1章を読んでください。通常の手法ではこれ以上なかなかよい結果や考察ができません。せっかくの業務知識を生かした“もっと気の利いた重回帰”をしたいとは思いませんか？統計モデリングはその問いに対する一つの解答です。
書籍の第8章では最尤推定でパラメータの値を推定しているようです。最尤推定で求めた値を信じて痛い目あったことある人いませんかね？僕だけですか、そうですか。BUGSやStanを使ってモデリングすることでパラメータの“分布”を求めることができます。色々なメリットがありますが、以下の久保先生のツイートがとてもよかったので参考までに引用しておきます。
library(lme4) がRから撤退するかも…とのこと．無くてもいいんじゃない…などとlme4にはまったく愛着のない私などは考えるのであった．だいたい，ややこしい統計モデルの最尤推定なんてアブないだけのような気がする．
— 久保拓弥 (@KuboBook) 2013, 12月 10
統計モデリングは和で分解が基本です。パラメータ同士の積とかあまり好ましくないです。もちろん仕方ない時はありますが。また書籍では行列形式（DLM; Dynamic Linear Model）で説明を展開していきますが、とてもスパースな行列ですし、BUGSやStanでは効果ごとにモデル式を立てた方が実行速度も速く、分かりやすいかと思います。
二次トレンドモデルは、「≈（にょろイコール）の左側の差分と右側の差分は似ているだろう（だが同じではない）」という経験からきています。真のトレンドは書籍の分析結果とほぼ同一になるように値をとっています。
この週効果の入れ方は統計数理研究所の赤池流です。微妙な週効果の変動にフィットできる反面、パラメータが多くて収束しにくくなるデメリットもあります。速度や収束しやすさを重視するならば、7日間の合計がピッタリ0となるような制約も考えられます。休日効果にかんしては、2012年・2013年のカレンダーを選んだせいでD2=1は3日しかありませんでした。最近は祝日は土日につなげるのが流行みたいです。祝日のデータは検索するといくつかヒットします。さらに会社が休みであろう三が日やGWまわりは祝日扱いとして追加しました。ちなみにデータを生成した時に使った真のパラメータの値は書籍にあわせて（b0, b1, b2）=（1, 0, 0.4）です。
天気のデータはこちらの出現確率からstochasticに生成しました。ほんとは天気のつながりとかあるのでそうは簡単ではないですが、ランダムよりマシだろうということで。雨のうち確率的に10%程度を大雨にしました。真のc_rainの値は書籍（Web）の分析結果にあわせて11.7です。
イベントに関するデータは全くよりどころがなかったので直感でテキトーに作りました。イベントのあった日数だけは書籍におおよそあわせてあります。Event_valに変換するところは勉強になりました。立ち上がるところの参加人数（8000人）と立ち上がる係数（0.00012）は書籍（Web）の分析結果にあわせてあります。はじめはこれらの値すらもStanで推定させようと気張りましたが死にました。パラメータが崖のようになっていて、値が少し変わると結果が大幅に変わるような場合はMCMCサンプリングによる推定は難しいと理解しています。またc_eventはイベントの影響の時間的変動を表す部分です。データがないかわりに、かなりイジワルな変動を作って与えました。
AR項は短期的（ここでは月のオーダー）な振動や変化を表現するのに便利です。この点は船渡川先生の本を読むと分かりやすかったです。
Stanコードは以下になります。
data { int<lower=1> N; int<lower=1> N_obs; int<lower=0> N_na; int<lower=1> Idx_obs[N_obs]; int<lower=1> Idx_na[N_na]; int<lower=0, upper=1> D1[N]; int<lower=0, upper=1> D2[N]; int<lower=0, upper=6> Wday[N]; real Rain_val[N]; real Event_val[N]; real<lower=0> Sale_obs[N_obs]; } parameters { real trend[N]; real s[N]; real<lower=0, upper=1> b1; real<lower=0, upper=1> b2; real<lower=0, upper=1> b3; real c_rain; real c_event[N]; real ar[N]; real c_ar[2]; real<lower=0> sale_na[N_na]; real<lower=0, upper=100> s_trend; real<lower=0, upper=100> s_s; real<lower=0, upper=100> s_event; real<lower=0, upper=100> s_ar; real<lower=0, upper=100> s_r; } model { real week[N]; real rain[N]; real event[N]; real sale_mu[N]; for(i in 3:N) trend[i] ~ normal(2*trend[i-1] - trend[i-2], s_trend); for(i in 7:N) s[i] ~ normal(-s[i-1]-s[i-2]-s[i-3]-s[i-4]-s[i-5]-s[i-6], s_s); for(i in 1:7) week[i] <- s[i] + D1[i]*b1*(s[i-Wday[i]]-s[i]); for(i in 8:N) week[i] <- s[i] + D1[i]*b1*(s[i-Wday[i]]-s[i]) + D2[i]*(b2*(s[i-Wday[i]-2]-s[i]) + b3*(s[i-Wday[i]-1]-s[i])); for(i in 1:N) rain[i] <- c_rain * Rain_val[i]; for(i in 2:N) c_event[i] ~ normal(c_event[i-1], s_event); for(i in 1:N) event[i] <- c_event[i] * Event_val[i]; for(i in 3:N) ar[i] ~ normal(c_ar[1]*ar[i-1] + c_ar[2]*ar[i-2], s_ar); for(i in 1:N) sale_mu[i] <- trend[i] + week[i] + rain[i] + event[i] + ar[i]; for(i in 1:N_obs) Sale_obs[i] ~ normal(sale_mu[Idx_obs[i]], s_r); for(i in 1:N_na) sale_na[i] ~ normal(sale_mu[Idx_na[i]], s_r); }
難しいところはあまりないと思います。Stanは欠損値のデータは別ループで回す必要があります。そのため3-6行目と58-61行目のような書き方が必要になります。観測データと欠損値のインデックスはRからデータとして渡してやります。あとb0,b1,b2のところは書籍の条件に忠実に実装しようとすれば以下のようになります。しかしながらデータ数が少ないために収束が若干悪くなったり、別にその制限はかけなくてもいいんじゃないかなという気もしたのではずしました。
parameters { ... real<lower=0, upper=1> b1; real<lower=0, upper=1> b2; real<lower=0, upper=1> sum_b2_b3; ... } transformed parameters { real<lower=0, upper=1> b3; b3 <- sum_b2_b3 - b2; }
昔はそう思っていましたが、Stanにある機能を使った方が高速かつ後で手間が少ない場合が多いです。
kickするRコードは以下になります。
library(rstan) rstan_options(auto_write = TRUE) options(mc.cores = parallel::detectCores()) d <- read.delim("input/data_20131226.txt", sep="\t") d[d$Month == 8, ]$Sale <- NA Idx.na <- which(is.na(d$Sale)) Idx.obs <- which(!is.na(d$Sale)) Sale.obs <- d$Sale[Idx.obs] data <- list( N = nrow(d), N_obs = length(Idx.obs), N_na = length(Idx.na), Idx_obs = Idx.obs, Idx_na = Idx.na, Wday = d$Wday, D1 = d$D1, D2 = d$D2, Rain_val = d$Rain_val, Event_val = d$Event_val, Sale_obs = Sale.obs ) fit <- stan(file='b5.stan', data=data, pars=c('trend', 's', 'b1', 'b2', 'b3', 'c_rain', 'c_event', 'ar', 'c_ar', 's_trend', 's_s', 's_event', 's_ar', 's_r'), iter=6000, warmup=1000, seed=123) save.image("output/result.Rdata")
結果以降は発表資料と重複しますのでやめておきます。
勉強会ではモデルの予測精度はどう評価するのか？という質問がありました。それは、例えば時系列の最後の方のデータだけは推定に使わないで、あとから評価に使います。具体的にはStanのgenerated quantitiesブロックの中で推定されたパラメータ＋分布にあてはめます。すると尤度が得られるのでそれ対数とったりして使うのがよくあるやり方です。尤度の代わりに二乗誤差とかでもよいと思います。ARIMAなどの手法と比べても精度はよいと思います。
統計モデリングに興味のある方はとても分かりやすい久保先生の緑本を読んで、BUGSやStanからはじめてみるのがよいと思います。また、2015/10から岩波データサイエンスが出版されました。vol.1とvol.6はここでとりあげた状態空間モデルについて入門と実装を扱った記事が豊富です。併せて読むとよいかと思います。
- 発売日: 2015/10/08
- メディア: 単行本（ソフトカバー）
2016.03.06追記
上記のStanコードを徹底的にvector化して高速化したのでStanコードをアップしておきます。参考になれば。
data { int<lower=1> N; int<lower=1> N_obs; int<lower=0> N_na; int<lower=1> Idx_obs[N_obs]; int<lower=1> Idx_na[N_na]; int<lower=0, upper=1> D1[N]; int<lower=0, upper=1> D2[N]; int<lower=0, upper=6> Wday[N]; vector[N] Rain_val; vector[N] Event_val; vector<lower=0>[N_obs] Sale_obs; } parameters { vector[N] trend; vector[N] s; real<lower=0, upper=1> b1; real<lower=0, upper=1> b2; real<lower=0, upper=1> b3; real c_rain; vector[N] c_event; vector[N] ar; real c_ar[2]; real<lower=0, upper=100> s_trend; real<lower=0, upper=100> s_s; real<lower=0, upper=100> s_event; real<lower=0, upper=100> s_ar; real<lower=0, upper=100> s_r; } transformed parameters { vector[N-6] sum_seg_s; vector[N] week; vector[N] rain; vector[N] event; vector[N] sale_mu; for (i in 7:N) sum_seg_s[i-6] <- sum(s[i-6:i]); for(i in 1:7) week[i] <- s[i] + D1[i]*b1*(s[i-Wday[i]]-s[i]); for(i in 8:N) week[i] <- s[i] + D1[i]*b1*(s[i-Wday[i]]-s[i]) + D2[i]*(b2*(s[i-Wday[i]-2]-s[i]) + b3*(s[i-Wday[i]-1]-s[i])); rain <- c_rain * Rain_val; event <- c_event .* Event_val; sale_mu <- trend + week + rain + event + ar; } model { trend[3:N] ~ normal(2*trend[2:N-1] - trend[1:N-2], s_trend); sum_seg_s ~ normal(0, s_s); c_event[2:N] ~ normal(c_event[1:N-1], s_event); ar[3:N] ~ normal(c_ar[1]*ar[2:N-1] + c_ar[2]*ar[1:N-2], s_ar); Sale_obs ~ normal(sale_mu[Idx_obs], s_r); } generated quantities { vector[N_na] sale_na; for (i in 1:N_na) sale_na[i] <- normal_rng(sale_mu[Idx_na[i]], s_r); }
2017.07.15追記
週効果を微小な変動ありから変動なし（合計が0になる制約）に置き換えた場合のコードを追加しました。少し柔軟さは落ちますが、高速になり収束も速くなります。また、lowerとupperを使った範囲制限よりもmodelブロックでnormalで縛った方がオススメです（「StanとRでベイズ統計モデリング」10章参照）。 そこでb1, b2, b3の事前情報分布を設定します。また、各標準偏差にも弱情報事前分布を設定します。これで1chainあたり40分ぐらいで計算できるようになりました（iter=6000, warmup=1000）。
data { int<lower=1> N; int<lower=1> N_obs; int<lower=0> N_na; int<lower=1> Idx_obs[N_obs]; int<lower=1> Idx_na[N_na]; int<lower=0, upper=1> D1[N]; int<lower=0, upper=1> D2[N]; int<lower=0, upper=6> Wday[N]; vector[N] Rain_val; vector[N] Event_val; vector<lower=0>[N_obs] Sale_obs; } parameters { vector[N] trend; vector[6] s_raw; real b1; real b2; real b3; real c_rain; vector[N] c_event; vector[N] ar; real c_ar[2]; real<lower=0> s_trend; real<lower=0> s_event; real<lower=0> s_ar; real<lower=0> s_r; } transformed parameters { vector[7] s; vector[N] week; vector[N] rain; vector[N] event; vector[N] sale_mu; s[1:6] = s_raw; s[7] = -sum(s_raw); for(i in 1:7) week[i] = s[Wday[i]+1] + D1[i]*b1*(s[1]-s[Wday[i]+1]); for(i in 8:N) week[i] = s[Wday[i]+1] + D1[i]*b1*(s[1]-s[Wday[i]+1]) + D2[i]*(b2*(s[6]-s[Wday[i]+1]) + b3*(s[7]-s[Wday[i]+1])); rain = c_rain * Rain_val; event = c_event .* Event_val; sale_mu = trend + week + rain + event + ar; } model { b1 ~ normal(0.5, 0.5); b2 ~ normal(0.5, 0.5); b3 ~ normal(0.5, 0.5); s_trend ~ normal(0, 10); s_event ~ normal(0, 10); s_ar ~ normal(0, 10); s_r ~ normal(0, 10); trend[3:N] ~ normal(2*trend[2:N-1] - trend[1:N-2], s_trend); c_event[2:N] ~ normal(c_event[1:N-1], s_event); ar[3:N] ~ normal(c_ar[1]*ar[2:N-1] + c_ar[2]*ar[1:N-2], s_ar); Sale_obs ~ normal(sale_mu[Idx_obs], s_r); } generated quantities { vector[N_na] sale_na; for (i in 1:N_na) sale_na[i] = normal_rng(sale_mu[Idx_na[i]], s_r); }
