(35連休)11日目：ARIMAモデルとdlm（動的線形モデル）による販売量予測モデルの検証④ ＋データサイエンティストという胡散臭い肩書きについて（余談）

これまでやってきたこと

下記の記事通り、Rのdlmとforecastパッケージを使って、BJsalesという150日間の日次販売データを用いて、週次の予測ローリングをやってみています。
http://hsuetsugu.hatenablog.com/entry/2014/08/25/152818
http://hsuetsugu.hatenablog.com/entry/2014/08/23/125238
http://hsuetsugu.hatenablog.com/entry/2014/08/22/164140

本日やりたいこと

「予測にいかす統計モデリングの基本」をぱらっと見返していたら、あるレストランでの２年間の日次販売金額データをインプットにしてモデリングしているのですが、欠損値のある期間についても推定されたトレンド成分が欠損値にふられずに算出されていたので、ちょっとトレースしてみることにしました。

予測にいかす統計モデリングの基本―ベイズ統計入門から応用まで (ＫＳ理工学専門書)

下記のサイトを主に参考にさせていただきましたので、正確にはトレースのトレースですが・・。（コードはここのを参照したので、この記事には貼付けません。）
RStanで『予測にいかす統計モデリングの基本』の売上データの分析をトレースしてみた

余談ですが・・・

今回のデータは仕事で扱ってるデータじゃないのに、結構仕事で使うようなデータとよく似ています。結構リアルな感じです。書籍から作ってくれた人すごい・・・！！
業務ではそもそも周辺イベントのデータって誰がどう取得するんだっけ、とか、データ収集の話が重いんだろうな・・。小売のチェーンストアとかになると店舗数が多くて商品数も多いので、各ノード（店×商品）について周辺イベントとか特に考えずに在庫でバッファリングするか突発的なイベントは現場業務でなんとかやりくりしてしまう（少なくとも過去実績が適切な切り口で可視化さえされていれば十分）だろうけど。そもそも店別カテゴリ別の適切な販売計画をどうたてるか、ということが肝になるので、使うとしたら人口動態みたいなマクロデータをもとにして、本部側と店側のデータドリブンなPDCAをどうまわしていくか、みたいな話が中心になりますが、その辺がいわゆる"データ分析"と"コンサル"との違いなのでしょう。そしてそのへんの違いが、"データサイエンティスト"とかいう肩書きが極めて胡散臭く聞こえる理由なのでしょうか。
自分としては、エンジニアリングおよび統計的な理解の欠如から、これを名乗れないと認識しているし、それ以前に、上記のような理由から名乗りたくない。同じように"データサイエンティスト"と名乗りたくない人は多くいるのではないかと思います。なんかかっこ悪いし、サイエンスかどうかなんて議論はほんとにどうでもいいし。
ただコンサルしている人の多くはデータ分析を大上段に構えすぎていて、定性的なヒアリングと目的は同じであることをあまりわかっていないし、分析が得意な人の多くは業務を知らなさすぎる上に客と話せなすぎる、というのは経験上間違えではなさそうで、これはコンサルとエンジニアの違いと近いものがありそうです。じゃあ自分はなんなのかというとはっきり表したい肩書きがないのですが、少なくとも事業に近いところをやろうとすると今はコンサルスキルよりも実装力だとか分析力が圧倒的に不足しているので、そこを今の期間にのばしておきたいと思っています。
それからソフトウェアからハードウェアへ。ここしばらくRaspberry Piから離れていますが、「予測にいかす統計モデリングの基本」にも最後に出てくる題材は粒子フィルタを使ってロボット制御する、というお話なので、この辺うまく組み合わせたものを9月中には何か形にしたいと思っています。

RStanのインストールとテスト(MAC OS X)

インストール

特にひっかかったところはなく、下記コマンドをRで実行して問題なくインストールできました。
参考URL:RStan Getting Started · stan-dev/rstan Wiki · GitHub

Sys.setenv(MAKEFLAGS = "-j4") 
source('http://mc-stan.org/rstan/install.R', echo = TRUE, max.deparse.length = 2000)
install_rstan()

実行テスト

これも上記のURL通り実行してみて、計算されているようでした。

library(rstan)

schools_code <- '
  data {
    int<lower=0> J; // number of schools 
    real y[J]; // estimated treatment effects
    real<lower=0> sigma[J]; // s.e. of effect estimates 
  }
  parameters {
    real mu; 
    real<lower=0> tau;
    real eta[J];
  }
  transformed parameters {
    real theta[J];
    for (j in 1:J)
      theta[j] <- mu + tau * eta[j];
  }
  model {
    eta ~ normal(0, 1);
    y ~ normal(theta, sigma);
  }
'

schools_dat <- list(J = 8, 
                    y = c(28,  8, -3,  7, -1,  1, 18, 12),
                    sigma = c(15, 10, 16, 11,  9, 11, 10, 18))

fit <- stan(model_code = schools_code, data = schools_dat, 
            iter = 1000, chains = 4)

「予測にいかす統計モデリングの基本」のトレースのトレース

データまで作成して公開していただいたのは本当にありがたいです・・・。今回のデータは、天気、周辺イベントなど外的要因も多く含まれてますが、まずは結果としての売上時系列がどうなっているのかを下記のような感じで把握しました。

まずは時系列の動きを把握

library(ggplot2)

# データ読み込み
data <- read.table('/[Directory]/data_20131226.txt',
                   sep='\t',header=T)
＃8月は欠損
data[data$Month == 8, ]$Sale <- NA

# 日付と週番号を付与
data$hiduke <- as.Date(paste0(data$Year,'-',data$Month,'-',data$Day))
data$ww <- floor((data$hiduke - as.Date('2012-01-01'))/7)

# 日次時系列を確認
ggplot(data,aes(x=hiduke,y=Sale)) + geom_line()
# 年間の曜日変動を確認
ggplot(data,aes(x=Wday,y=Sale,color=factor(ww))) + geom_line()

• 日次の販売金額時系列

• 曜日別週別販売金額

イベント効果と天気の影響の概観確認

あとは、イベントの開催頻度だとか、その効果の度合いをなんとなくみておきます。

ggplot(data,aes(x=Head.count,y=Sale,color=factor(Weather))) + geom_point()

• 日次のイベント参加人数に対する売上金額（色は天気）

stanでモデリング

基本的な考え方として、売上データをトレンド（長期変動）と曜日変動（週効果）と天気（雨効果）と周辺イベント効果と自己回帰(AR)項と残差の和で表します。
曜日変動は単純に平日と休日というだけでなく、祝日の効果と祝前日の効果を土日、金土とどれだけ似ているのかをモデリングできるようにしていて勉強になりました。週単位だったらMDカレンダで週番号さえあわせておけば業務的にはここまでしなくてもいいかもしれません。

とりあえず上記のURLからコードをもってきて実行してみましたが、計算実行に時間がかかるようので、今日はここまでとします。
ただ、下記の部分でRからstanに欠損値の期間を明示的に渡しており、それが必要となると使いにくい・・・。

Idx.na <- which(is.na(d$Sale))
Idx.obs <- which(!is.na(d$Sale))
Sale.obs <- d$Sale[Idx.obs]

欠損値なのか、ほんとに売上がなかったのか判別が難しい場合はどのようにしたらいいのでしょうか。今回のケースだと欠損値であることは売上データからみて明確ですが、グローバルで多店舗で展開しているようなケースで、店舗ごとに店頭欠品がある期間がそれぞれいくつかある場合になにかやりようがないでしょうか。まあ店頭欠品は各店で在庫が0になった時とかって判別できるので、そういうデータを渡す、ということにすれば可能ですが、そんな感じでやるのかな。

また続きます・・・。