7/15 multiple change points detection of CVR by pymc3

この記事では、多分一番コードがシュッとする感じに書けるpymc3を用いた複数の変化点検出を紹介します。
また今回採用する題材は、いわゆる"率"の変化点検出であるため、web業界などでのCVRやCTRの時系列的変化として比較的汎用なものです。
たとえば、施策を打った際のCVRの変化などをこういった方法で検出することは非常に自然な考え方です。

なお、変化点の数は分析者が置きで決めるタイプのものとして、変化点の数も一緒にベイズ推定するさらにハードコアなやつは今後の課題とします。

使うデータは
xn--v6q193b2rby9ymh4b.realtime-chart.info
の支持率データをお借りすることにします。

まず、簡単にプロットしてみましょう。

いくつか急激に上がったり下がったりしている部分、つまりイベント効果が見受けられます。
こいつを、

df = df.resample('1W').sum()

とすれば、weeklyのデータにresamplingすることができます。（ただし、indexをTimeIndexにする必要があります。）
resamplingしたものはこちらです。

さっきよりギザギザは減ってますが、それでもイベント効果が乗ってる感じはしますね。

今回やりたい分析は、こういったイベント効果を取り除いて、"支持率のベースラインの変化"を取り出したいとしましょう。
本当は分析結果がイベント効果に吸い付かなかったためにこういうお題目にしているわけでは決してありません。

さて、そのためのモデルを組んでいきましょう。
1つの点の変化点検出はpymc3のチュートリアルにもあるので、詳しく動作させるのはそちらやgistを見ていただくとして、こういう感じにモデルを書けば動きます。

変化点検出はこちら。
Getting started with PyMC3 — PyMC3 3.5.rc1 documentation

weeks = df.reset_index().index
with pm.Model() as model:
    switchpoint = pm.DiscreteUniform('sp', 
                                     lower=weeks.min(),
                                     upper=weeks.max())
    p1 = pm.Beta('p1',alpha=1,beta=1)
    p2 = pm.Beta('p2',alpha=1,beta=1)
    rate = pm.math.switch(switchpoint>=weeks, p1, p2)
    k = pm.Binomial('k',observed = df['Yes'].values,
                        n = df['Count'].values, 
                        p = rate)

で、こういう感じの結果が得られます。

これをplotするとこんな感じです。

一回上がって一気に下がったあたりでベースラインが下がってるっぽいわよ、という分析結果になっています。
これだと流石にモデリングが足りてない感じがするので、もうちょっと変化点の数を増やしたいよね、というのが今回の本題です。

しかし、上のコードを見ればわかるように、pymc3の変化点検出のチュートリアルは変化点の数に対して線形にコードが長くなっていく実装になっていて、これはイケているとは言えません。
そこで、今回は変化点の数に対して線形にコードが長くならない実装を紹介します。

K = 5
with pm.Model() as model:
    switchpoint = [pm.DiscreteUniform('sp0',lower=weeks.min(),upper=weeks.max())]
    for i in range(1,K):
        switchpoint = switchpoint + [pm.DiscreteUniform(f'sp{i}',lower=switchpoint[i-1] + 1,
                                                        upper=weeks.max())]
    p = pm.Beta('p',1,1,shape = K+1)
    rates = [pm.math.switch(switchpoint[i] >= weeks, p[i], p[i+1]) for i in range(K)]
    k = pm.Binomial('k',observed = df['Yes'].values,
                        n = df['Count'].values, 
                        p = rates)

こういう感じに書けば、見事変化点がいくつ増えてもコードが長くなりません。
ハマったポイントは、最初はpm.math.switchのほうで変化点の位置をコントロールしようとしたのですが、そうでなくて事前分布のほうで変化点の位置をコントロールする実装が正しいというところです。
これはいかにもベイズっぽくてな、なるほど……となりました。

こやつにMCMCをかけた結果がこちらです。

で、plotしたのがこちらです。

うーん、まあそこそこベースラインを取り出せていると言っていいんじゃないでしょうか。
このモデルのチューニングの特徴としては、変化点の事前分布を結構強めに置かないといまいちワークしないことでした。
実際にK=8とかでやっちゃうと、後ろの方に変化点が詰まっちゃってモデルの推定が失敗します。
推定がうまくいかないときは、目に見えてMCMCが遅くなるので、そのへんはpymc3と対話してください。

gistはこちら。
gist.github.com