「Rで学ぶベイズ統計学入門」の例題を解きながら勉強していくシリーズです。今回はこんなあるある系の問題です。
ある販売店の営業マンA君とB君はそれぞれ別のエリアを担当している。彼らは朝から晩まで担当エリアのお客さんを回っており、その営業成績は「今週の受注件数」として数値化し管理されている。
今年こそはマネージャに昇格したい野心に燃えたA君は、さらに売上を加速させるべく、彼のエリア限定で新たな広告戦略を展開することを、上司である営業部長に提案した。
A君は普段からB君よりも平均1.5倍ぐらいの注文を獲っていたが、「4週間広告を投入してB君の1.6倍以上の結果を出す」ことを部長に宣言した。部長は快くこれを承認した。
ところが4週間の合計受注数は、A君が260、B君が165であり、1.58倍という、なんとも微妙な結果に終わってしまった。このままでは部長との約束をコミットできない。
もしかしたら広告の効果は実際にあったのに、たまたま何かのタイミングで買い控えがあったり、逆にB君の受注がたまたまいつもより多かった可能性もある。何とかして部長を納得させるうまい理由を考えられないだろうか？（※題意と数字の一部を改変させていただきました。）
こうなったらもう潔く諦めて次の手を考えればいいと思うのですが、いろんな大人の事情でそうもいかない場合もあります。ここはA君のためにうまい言い訳を考えてあげましょう。こんな言い訳はどうでしょうか？
「普段がB君の1.5倍ぐらいなので、1.6倍をひとつのベンチマークとして上回れば広告効果ありと考えました。お客さんからもそれなりの反響があり手ごたえはあったのですが、今回は1.58倍という惜しい結果になりました。ただ分析をしたところ、広告の効果は確実に出ていることがわかりました！」
では、どういう分析をしたらそういうことが言えるのでしょうか？
過去の実績を分析して事前分布を決める
A君、B君の過去1年間（広告投入前）の週間受注件数はおそらくこんな感じです。A君は平均60ぐらい、B君は平均40ぐらいなので、ざっくりA君はB君の1.5倍ぐらいで推移しています。
例題では、これを踏まえて事前分布および尤度を以下のように設定することになっています。
事前分布
尤度
- A君およびB君の「実際の受注数」$y_A, ~ y_B$は、それぞれ平均受注数$\lambda_A, ~ \lambda_B$を平均としたポアソン分布にしたがう
ところで、「週あたりの平均受注数」が「確率分布にしたがう」というのは違和感があるかもしれませんが、これは、たとえば「A君の平均的パフォーマンス」あるいは「担当エリアにおける需要の水準」みたいな概念をあらわすための数字と言い換えてみると解釈しやすいと思います。
この数字を期待値として、いろいろなランダム要因（天気やイベントやお客さんのその日の気分やそのほか説明しきれないさまざまな事象）が作用することで「実際の受注数」という結果になると考えられます。「実際の受注数」が発生する仕組みをここではポアソン分布と仮定しています。
一般人的な感覚では「過去実績の平均値のほうがまぎれもない事実」と思うかもしれません。しかし、ベイズ推定の考え方では、「データそのもの」よりも「データが発生する仕組み」のほうに興味があります。データはランダム作用を受けた結果に過ぎません。それよりも「需要水準」みたいな「直接観測できないけれども重要な値」を解き明かすことが大事です。ベイズ統計ではこういう観測できない値は「確率分布」で表現することになっています。
$\lambda_A, \lambda_B$の結合事前分布をプロットするとこんな感じです。黄色い線はそれぞれの周辺分布の95%信用区間です。
広告投入前、A君の平均受注はB君の何倍か？
$A君の平均受注数 / B君の平均受注数 = \lambda_A / \lambda_B = u$と表記することにします。広告投入前後で$u$を比較して「アップした確率が何パーセント」みたいに言えれば、たぶん部長も納得するはずです。
$\lambda_A, \lambda_B$は確率変数なので$u$も確率変数です。まずは、広告投入前の$u$がどういう分布になっているか確認する必要があります。これは事前分布と同じと考えていいと思います。
ガンマ分布にしたがう独立した確率変数を割り算するとどういう形の確率分布になるのか計算するのは大変そうです。数学音痴の私の手には負えません。そこでここは乱数を使った力技で求めたいと思います。
# ガンマ事前分布のパラメータ
alpha <- c(144, 100)
beta <- c(2.4, 2.5)
# 1000個のuをランダムサンプリング
u.prior <- rgamma(1000, alpha[1], beta[1]) / rgamma(1000, alpha[2], beta[2])
# サンプルのパーセンタイル値
quantile(u.prior, c(0.025, 0.5, 0.975))
# 2.5% 50% 97.5%
# 1.172460 1.501580 1.914928
# サンプルのヒストグラム = uの事前確率分布
hist(u.prior, freq = FALSE, breaks = 50,
main = expression(paste(lambda[1] / lambda[2], "の事前分布")),
xlab = expression(lambda[1] / lambda[2]), ylab = "確率密度")
B君と比較したA君の平均受注は、中央値で1.5倍、95%信用区間は1.2倍から1.9倍ぐらいになります。
広告投入した結果どうなった？
4週間広告を投入した結果を、ベイズの定理を使ってA君、B君の平均受注数に反映したいと思います。
実はこの計算はめちゃくちゃ簡単です。$データy$が$ポアソン分布(t\lambda)$にしたがい、$\lambda$が$ガンマ事前分布(\alpha, \beta)$にしたがう場合、λの事後分布は$ガンマ分布(\alpha + y, \beta + t)$になります。ですので、嬉しいことに計算が必要なのは$\alpha + y$と$\beta + t$だけです。
ポアソン分布に対するガンマ事前分布みたいに、事後分布が事前分布と同じ分布型になることを「共役事前分布」といいます。これを使えると計算が楽です。
で、先ほどと同様に$\lambda$の事後分布から$u$の分布を求めてみます。
# 4週間のA君, B君の受注数量
y <- c(260, 165)
# 事後分布から1000個のuをランダムサンプリング
u.post <- rgamma(1000, alpha[1] + y[1], beta[1] + 4) /
rgamma(1000, alpha[2] + y[2], beta[2] + 4)
# サンプルのパーセンタイル値
quantile(u.post, c(0.025, 0.5, 0.975))
# 2.5% 50% 97.5%
# 1.319093 1.552309 1.805708
# サンプルのヒストグラム = uの事後確率分布
hist(u.post, freq = FALSE, breaks = 50,
main = expression(paste(u == lambda[1] / lambda[2], "の事後分布")),
xlab = expression(u == lambda[1] / lambda[2]), ylab = "確率密度")
中央値は1.55倍、95%信用区間は1.3倍から1.8倍になりました。
それで、効果あったの？
さて佳境です。部長を説得するためには、広告投入前と投入後の$u$を比較してアップしたことを示す必要があります。
そこでたとえば、$v = (広告投入後のu) / (広告投入前のu)$という指標を考えて、$v > 1$なら広告効果ありと主張するのはどうでしょうか。
広告投入前の$u$は事前分布で置き換えて考えていいと思うので、また力技でそれぞれの$u$からランダムサンプルから$v$を計算してみます。
v <- sample(u.post, 1000, replace = TRUE) / sample(u.prior, 1000, replace = TRUE)
# vのパーセンタイル値
quantile(v, c(0.025, 0.5, 0.975))
# 2.5% 50% 97.5%
# 0.7620471 1.0206326 1.3906756
# v > 1の確率
sum(v > 1) / length(v)
# [1] 0.555
「効果ありの確率55%」。言い訳に使うにはちょっと弱すぎです。残念ですが、A君には潔く諦めてもらいましょう。
まとめ
問題をやってみて、「B君の何倍」という主張には違和感を感じました。B君と比較しないと物を申せないA君の人間性にも問題を感じますが、それ以上に、
- 広告効果を金額に換算できないので、費用に対してリターンが十分かどうかという議論ができない
- トレンドや季節変動をモデル化できていないので、それらの要素があった場合、B君の何倍という指標はロバストでない
- B君のエリアで何か（特需、B君の急成長、競合 etc.）が起こった場合、B君と比較することに意味がない
などぱっと思いつくだけでいろいろ問題があります。
これらをうまくモデル化するなら、広告効果をダミー変数として階層モデルやカルマンフィルタなどで推定することになるのでしょうか。このあたり今後勉強していきたいと思います。詳しい方はぜひ教えてください。