Prophetのモデル式を1から理解する

この記事はBASE Advent Calendar 2019の20日目の記事です。

devblog.thebase.in

DataStrategyの岡が担当します。

Prophet is 何？

ProphetはFacebook社製の時系列予測ライブラリです。RとPythonから利用でき下記gitで公開されています。 https://github.com/facebook/prophet

分析者仲間の間で「時系列予測ならまずこれを使っとけ」と言われるくらい高精度らしいのですが、私自身がイマイチ理論を把握してない & ググってもさらっとした解説の日本語ドキュメントしかない印象です。

Prophetの元となる論文は下記にて公開されています。 https://peerj.com/preprints/3190.pdf

冒頭だけ読むと、時系列分析の知見のない人でもドメイン知識を組み込みながら予測ができるようなツールを目指して開発されたようです。論文タイトルが「Forecasting at Scale」となっていて「Scale」というのはここでは「誰でも使える、スケールしやすい予測ツール」みたいな意味で使われています。

論文に書かれているProphetのモデル式をしっかり理解できるように丁寧になぞっていこう、というのがこの記事の趣旨です。

モデル式の概要

時系列データは、トレンド + 季節要因 + ノイズ などの複数の要素から成り立っていて、これらの要素に分解することで理論値の予測ができる、という考え方があります。

Prophetでは、時系列は下記のような構成要素をもつと捉えています。

• g(t):トレンド関数$g(t):トレンド関数$
• s(t):季節変化$s(t):季節変化$
• h(t):休日効果$h(t):休日効果$
• εt:誤差項$\epsilon t:誤差項$

さらに、時系列はこれらの要素の和と捉え下記のようなモデル式を組み立てています。

y(t)=g(t)+s(t)+h(t)+εt.$y(t)=g(t)+s(t)+h(t)+\epsilon t.$

このモデル式を理解するため、誤差項以外の3要素を一つずつ見ていきます。

1. g(t)$g(t)$: トレンド関数

まず g(t)$g(t)$ ですが

1. ロジスティック非線形トレンド
2. 線形トレンド

の2種類があります。1のほうから読み解いてみます。

1-1. ロジスティック非線形トレンド

ごく単純化すると下記のように表されます。

g(t)=C1+exp(−k(t−m))(1)$\begin{array}{}\text{(1)}& {\displaystyle g(t)=\frac{C}{1+\exp (-k(t-m))}}\end{array}$

C:収容上限k:成長率m:オフセット項$C:収容上限k:成長率m:オフセット項$

これはロジスティック曲線をベースに作られているそうなので比べてみます。

y(t)=11+exp(−t)(2)$\begin{array}{}\text{(2)}& {\displaystyle y(t)=\frac{1}{1+\exp (-t)}}\end{array}$

Prophetのトレンド関数(1)式と比べると、 C=1$C=1$ となり exp$\exp$ の中身がシンプルになっています。 この関数の動きを把握するため、t$t$ を動かすとき y$y$ がどのようになるか図示してみます。

• t→−∞$t\to -\mathrm{\infty }$ のとき y=0$y=0$
• t→+∞$t\to +\mathrm{\infty }$ のとき y=1$y=1$

となり出力 y$y$ の下限と上限が決まっています。この基本となる式(2)に一つずつ

• C:収容上限$C:収容上限$
• k:成長率$k:成長率$
• m:オフセット項$m:オフセット項$

の3要素を継ぎ足してトレンド関数の形に近づけていきます。

生物の個体数やプログラムのバグ発見数など、初めは少ない → 途中は多くなる → その後また少なくなる という流れを持つ現象はロジスティック曲線で説明されることがあります。 生物の個体数などには何かしら上限数があると考えられ、トレンド関数では C:収容上限$C:収容上限$ をつけることでこの上限を設定しています。

トレンド関数にならって、ロジスティック関数（式(2)）の分子を C$C$ に変えてみます。

y(t)=C1+exp(−t)${\displaystyle y(t)=\frac{C}{1+\exp (-t)}}$

この C$C$ の値を1, 2, 3,...と変化させると先ほどの曲線は下図のように変化します。

y$y$ のmax値が C$C$ の値に等しくなり、これで上限を設定できるようになりました。 続いて k:成長率$k:成長率$ を加えてみましょう。

y(t)=C1+exp(−kt)${\displaystyle y(t)=\frac{C}{1+\exp (-kt)}}$

C$C$ を固定した状態で、 k$k$ の値を-1, 0.5, 1,...と変化させると曲線は下図のように変化します。

k$k$ が大きいほど曲線が急になり、小さいほど緩やかになるのが見て取れます。k=−1$k=-1$ の時はもはや減少関数となり、 k$k$ は「傾き」のようなものだと解釈できます。

生物の個体数でいうなら、成長スピードの早い群はより早く上限に達し、成長の遅い群は上限に達するのも遅い、といった知見をパラメタ k$k$ で表現できます。

続いて、m:オフセット項$m:オフセット項$ を追加していきます。

y(t)=C1+exp(−k(t−m))${\displaystyle y(t)=\frac{C}{1+\exp (-k(t-m))}}$

これで冒頭のトレンド関数（式(1)）と同様の式になりました（y(t)→g(t)$y(t)\to g(t)$と読み替えてみてください）。これも m$m$ を変化させながら曲線の動きを見てみましょう。

m$m$ は t$t$ の値をダイレクトに減算する式になっているので単に曲線が左右に動くだけですね。言い換えると同じ t$t$ でも y$y$ の値を上下させる働きを持っています。線形回帰でいうところの切片のようなものとイメージして下さい。

以上で下記の式の説明は終わりです。

g(t)=C1+exp(−k(t−m))(2)$$\begin{gathered} \begin{array}{}\text{(2)}& {\displaystyle g(t)=\frac{C}{1+\exp (-k(t-m))} \\ }\end{array} \end{gathered}$$

ただ、これはあくまでトレンド関数を単純化した式です。 論文ではさらに、C:収容上限$C:収容上限$ も k:成長率 $k:成長率$ も一定ではなく時間によって変化するものと捉え、

C→C(t)k→k+a(t)Tδ$$\begin{gathered} {\displaystyle C\to C(t) \\ k\to k+\mathbf{a}(t{)}^{\mathrm{T}}\mathit{\delta }} \end{gathered}$$

と展開しています。

C(t)$C(t)$は各 t$t$ 時点で異なる上限を設定できるようになった、というだけの話なので説明はこれくらいに留め、k+a(t)Tδ$k+\mathbf{a}(t{)}^{\mathrm{T}}\mathit{\delta }$ のほうをじっくり見ていきます。

k:成長率$k:成長率$ は転換期のようなものを迎えた場合かなり異なったレートに変化するはず、というニュアンスを数式で表現したいとします。仮にそのような転換点がS個あったとして、そのときのt$t$時点を

sj（j=1,...,S)${\displaystyle s_j（j=1,...,S)}$

と表します。また、その sj$s_j$ における成長率を調節する変更率として、δj${\delta }_j$ を定義します。sj$s_j$ に対応してS個分あるので

δ∈RS${\displaystyle \mathbf{\delta }\in {\mathbb{R}}^S}$

というベクトル（慣例にならってベクトルは太字で表記）も定義しておきます。 ある時間 t$t$ における成長率は基本レート k$k$ と t$t$ 時点までに出現した変更率 δj${\delta }_j$ の総和として下記の式で表されます。

k+∑j:t>sjδj(3)$\begin{array}{}\text{(3)}& {\displaystyle k+\sum _{j:t>s_j}{\delta }_j}\end{array}$

論文ではさらに ∑j:t>sjδj$\sum _{j:t>s_j}{\delta }_j$ の部分を扱いやすいベクトルで表記するため、

a(t)∈{0,1}S${\displaystyle \mathbf{a}(t)\in \{0,1{\}}^S}$

というS個分の0または1の要素から成るベクトルを用意し、その各要素を

aj(t)={1,ift≧sj,0,otherwise.${\displaystyle a_j(t)=\{\begin{array}{l}1,\mathrm{i}\mathrm{f}t\geqq s_j,\\ 0,\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}.\end{array}}$

と定義しています。例えば全部で10時点のt$t$があって、そのうち3回の変更点sj$s_j$が3,5,7番目に起きたとすると a(t)$\mathbf{a}(t)$ の内訳は下記のようになります。

a(t)=⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪[000],ift=1[000],ift=2[100],ift=3[100],ift=4[110],ift=5[110],ift=6[111],ift=7[111],ift=8[111],ift=9[111],ift=10${\displaystyle \mathbf{a}(t)=\{\begin{array}{l}[\begin{array}{ccc}0& 0& 0\end{array}],\mathrm{i}\mathrm{f}t=1\\ [\begin{array}{ccc}0& 0& 0\end{array}],\mathrm{i}\mathrm{f}t=2\\ [\begin{array}{ccc}1& 0& 0\end{array}],\mathrm{i}\mathrm{f}t=3\\ [\begin{array}{ccc}1& 0& 0\end{array}],\mathrm{i}\mathrm{f}t=4\\ [\begin{array}{ccc}1& 1& 0\end{array}],\mathrm{i}\mathrm{f}t=5\\ [\begin{array}{ccc}1& 1& 0\end{array}],\mathrm{i}\mathrm{f}t=6\\ [\begin{array}{ccc}1& 1& 1\end{array}],\mathrm{i}\mathrm{f}t=7\\ [\begin{array}{ccc}1& 1& 1\end{array}],\mathrm{i}\mathrm{f}t=8\\ [\begin{array}{ccc}1& 1& 1\end{array}],\mathrm{i}\mathrm{f}t=9\\ [\begin{array}{ccc}1& 1& 1\end{array}],\mathrm{i}\mathrm{f}t=10\end{array}}$

結局 a(t)$\mathbf{a}(t)$ は各変更点 sj$s_j$ を迎えるごとにひとつずつ変更点のフラグが立つだけのベクトルと言えます。これと先ほど定義した変更率 δj${\delta }_j$ のベクトル、δ∈RS$\mathit{\delta }\in {\mathbb{R}}^S$ を組み合わせると式(3)は、

k+∑j:t>sjδj=k+a(t)Tδ${\displaystyle k+\sum _{j:t>s_j}{\delta }_j=k+\mathbf{a}(t{)}^{\mathrm{T}}\mathit{\delta }}$

と変形でき、「一定ではなく時間によって変化する成長率」を表現することができました。 ここまでの式をトレンド関数に反映してみましょう。

g(t)=C(t)1+exp(−(k+a(t)Tδ)(t−m))(4)$$\begin{gathered} \begin{array}{}\text{(4)}& {\displaystyle g(t)=\frac{C(t)}{1+\exp (-(k+\mathbf{a}(t{)}^{\mathrm{T}}\mathbf{\delta })(t-m))} \\ }\end{array} \end{gathered}$$

これでt軸の転換点ごとに曲線の傾きが δj${\delta }_j$ ぶん変化していれば良いわけですが、実はこのままだと曲線の変化がなだらかになりません。 この成長率の変化によってどのような曲線が描かれるか可視化してみましょう。

import numpy as np
from matplotlib import pyplot as plt

T = np.linspace(-6, 6, 100)
S = np.array([-3, 1, 3])  # -3, 1, 3の時点で成長率に変化が起きる 
delta = np.array([0.1, 0.3, -0.6])  # 各S時点での成長率の変更率 

def logistic_trend(T, S, delta, k=1, C=1, m=0):
    a = np.vstack([np.where(S < t, 1, 0) for t in T])
    y = C / (1 + np.exp(-(k + (a * delta).sum(axis=1)) * (T - m)))
    return y

out = logistic_trend(T, S, delta, k=0.1, m=0)
plt.plot(T, out)
plt.xlabel("t", fontsize=16)
plt.ylabel("y", fontsize=16)

# plot change point
ymax = out.max()
ymin = out.min()
plt.vlines(
    S,
    ymin=ymin, ymax=ymax,
    linestyle='dashed',
    color='gray',
    label='change point'
    )

plt.legend()
plt.show()

変化点ごとに曲線が切れてしまってます。

さきほど k:成長率$k:成長率$ は「傾き」のようなもので、m:オフセット項$m:オフセット項$は「切片」のようなものだと説明しました。y=ax+b$y=ax+b$ のようなシンプルな一次関数と同じように考えて欲しいのですが、傾きaは　

a={1x≧−2−1x<−2${\displaystyle a=\{\begin{array}{l}1x\geqq -2\\ -1x<-2\end{array}}$

のようにxの領域によって変化し切片bは一定で0として可視化すると、さきほどと同じように途切れた直線になってしまいます。

T1 = np.linspace(-6, -2, 30)
T2 = np.linspace(-2, 6, 70)

def sample_linear(T, a=1):
    y = a * T
    return y

out1 = sample_linear(T1, a=1)
out2 = sample_linear(T2, a=-1)
plt.plot(
    np.hstack([T1, T2]),
    np.hstack([out1, out2]),
)
plt.xlabel("t", fontsize=16)
plt.ylabel("y", fontsize=16)
plt.show()

この場合は切片bを調節することで連続した直線が得られますが、話をトレンド関数に戻して考えるとオフセット項mを調整すれば連続した曲線が得られそうです。

転換点S個だけ調節する値が必要なので δ$\mathit{\delta }$ と同じく

γ∈RS${\displaystyle \mathit{\gamma }\in {\mathbb{R}}^S}$

というベクトルを用意し、その内訳を下記のように定義します。

γj=(sj−m−∑l<jγl)(1−k+∑l<jδlk+∑l≦jδl)${\displaystyle {\gamma }_j=(s_j-m-\sum _{l<j}{\gamma }_l)(1-\frac{k+\sum _{l<j}{\delta }_l}{k+\sum _{l\leqq j}{\delta }_l})}$

これを用いて式(4)のオフセット項mを調整するとトレンド関数は最終的に下記のようになります。

g(t)=C(t)1+exp(−(k+a(t)Tδ)(t−(m+a(t)Tγ))(5)$\begin{array}{}\text{(5)}& {\displaystyle g(t)=\frac{C(t)}{1+\exp (-(k+\mathbf{a}(t{)}^{\mathrm{T}}\mathit{\delta })(t-(m+\mathbf{a}(t{)}^{\mathrm{T}}\mathit{\gamma }))}}\end{array}$

この修正を先ほどのコードに加えると下記のようにグラフ化できます。

T = np.linspace(-6, 6, 100)
S = np.array([-3, 1, 3])
delta = np.array([0.1, 0.3, -0.6])

def logistic_trend(T, S, delta, k=1, C=1, m=0):
    a = np.vstack([np.where(S < t, 1, 0) for t in T])
    gamma = np.zeros(S.shape)
    for j in range(0, gamma.shape[0]):
        gamma[j] = (S[j] - m - gamma[:j].sum()) * (1 - ((k + delta[:j].sum()) / (k + delta[:j + 1].sum())))
    y = C / (1 + np.exp(-(k + (a * delta).sum(axis=1)) * (T - (m + (a * gamma).sum(axis=1)))))
    return y

out = logistic_trend(T, S, delta, k=0.1, m=0)
plt.plot(T, out)
plt.xlabel("t", fontsize=16)
plt.ylabel("y", fontsize=16)

# plot change point
ymax = out.max()
ymin = out.min()
plt.vlines(
    S,
    ymin=ymin, ymax=ymax,
    linestyle='dashed',
    color='gray',
    label='change point'
    )

plt.legend()
plt.show()

これで連続した曲線が得られました。

1-2. 線形トレンド

線形トレンドは下記の式で表されます。

g(t)=(k+a(t)Tδ)t+(m+a(t)Tγ)${\displaystyle g(t)=(k+\mathbf{a}(t{)}^{\mathrm{T}}\mathit{\delta })t+(m+\mathbf{a}(t{)}^{\mathrm{T}}\mathit{\gamma })}$

式(5)と見比べると、expの中身が出てきて成長率の部分はロジスティックの時と同様の式です。 オフセット項の部分だけ異なっていて、ここでは

γj=−sjδj${\displaystyle {\gamma }_j=-s_j{\delta }_j}$

という定義の γ$\mathit{\gamma }$ ベクトルで連続した直線が得られるよう調整されています。 参考までにコードを書くと下記のようになります。

T = np.linspace(-6, 6, 100)
S = np.array([-3, 1, 3])
delta = np.array([0.1, 0.3, -0.6])

def linear_trend(T, S, delta, k=1, m=0):
    a = np.vstack([np.where(S < t, 1, 0) for t in T])
    gamma = -S * delta
    y = (k + (a * delta).sum(axis=1)) * T + (m + (a * gamma).sum(axis=1))
    return y

out = linear_trend(T, S, delta, k=0.1, m=0)
plt.plot(T, out)
plt.xlabel("t", fontsize=16)
plt.ylabel("y", fontsize=16)

# plot change point
ymax = out.max()
ymin = out.min()
plt.vlines(
    S,
    ymin=ymin, ymax=ymax,
    linestyle='dashed',
    color='gray',
    label='change point'
    )

plt.legend()
plt.show()

1-3. 変化点の自動検出

変化点 sj$s_j$ はユーザー自身で設定することもできますが、スパース推定のようなことをして自動検出も可能です。 論文によると変化点 sj$s_j$ の上限数を多めにとり、各点の変更率 δ$\mathit{\delta }$ に対し

δj∼Laplace(0,τ)${\displaystyle {\delta }_j\sim \mathrm{L}\mathrm{a}\mathrm{p}\mathrm{l}\mathrm{a}\mathrm{c}\mathrm{e}(0,\tau )}$

という事前分布を仮定すれば良いとあります。Laplace(0,τ)$\mathrm{L}\mathrm{a}\mathrm{p}\mathrm{l}\mathrm{a}\mathrm{c}\mathrm{e}(0,\tau )$ はラプラス分布のことですが、その形状を確認してみましょう。

正規分布よりも0付近の値が出現しやすくなっている、という特徴がありそうです。 さらに τ$\tau$ の部分を変化させると下記のような分布が得られます。

τ$\tau$ が0に近づくにつれ、ほとんど0の値しか出現しない（スパースな）分布になっていることが見て取れます。
ここまでの話をまとめると、、、変更率 δj${\delta }_j$ にラプラス分布を仮定すると、多めの変更点をとっておいても変更率はほぼ0になり、かつ稀に大きな変更率が発生するという事象を再現することができます。また、変更率の大きさそのものは τ$\tau$ を小さく調整することで抑えることが可能になります。

1-4. トレンド関数の予測

ここまで扱ってきた変更率 δj${\delta }_j$ ですが、実際の予測の際にはどのような値を取ると良いでしょうか。論文を読むと、時系列データを予測する際に変更率 δj${\delta }_j$ を下記のようにシミュレーションさせるとあります。

まずラプラス分布のスケール τ$\tau$ を決定するためベイズ推定で事後分布を得るか、そうでなければ最尤推定的に解いて λ=1S∑Sj=1|δj|.$\lambda =\frac{1}{S}\sum _{j=1}^S|{\delta }_j|.$ を分布のスケールとします。この場合の λ$\lambda$ は過去に出現した変更率 δj${\delta }_j$ の絶対値の平均です。

過去の時系列の長さが T$T$ 個、そのうち成長率の変更のあった時点が S$S$ 個と定義したので、変更点の発生確率は ST$\frac{S}{T}$ 、発生しない確率は T−ST$\frac{T-S}{T}$ と言えます。 これらを踏まえ、論文では将来の δj${\delta }_j$ について下記のように定義しています。

∀j>T,{δj=0w.pT−ST,δj∼Laplace(0,λ)w.pST.${\displaystyle \mathrm{\forall }j>T,\{\begin{array}{l}{\delta }_j=0\mathrm{w}.\mathrm{p}\frac{T-S}{T},\\ {\delta }_j\sim \mathrm{L}\mathrm{a}\mathrm{p}\mathrm{l}\mathrm{a}\mathrm{c}\mathrm{e}(0,\lambda )\mathrm{w}.\mathrm{p}\frac{S}{T}.\end{array}}$

左の式ですが、∀$\mathrm{\forall }$ は全称記号なので T$T$ より大きい任意の j$j$ つまり未来に起きるすべての変更点について、という意味になります。右の式は、w.p$\mathrm{w}.\mathrm{p}$ は with probabilityの略なので

• T−ST$\frac{T-S}{T}$ の確率で δj=0${\delta }_j=0$
• ST$\frac{S}{T}$ の確率で δj${\delta }_j$ は Laplace(0,λ)$\mathrm{L}\mathrm{a}\mathrm{p}\mathrm{l}\mathrm{a}\mathrm{c}\mathrm{e}(0,\lambda )$ の分布に従う乱数

という意味になります。

まとめると、未来の変更率 δj${\delta }_j$ を求めるには、まず ST$\frac{S}{T}$ の確率で δj${\delta }_j$ が0になるかラプラス分布に従う乱数となるかが決まり、ラプラス分布に従う場合は その分布の乱数が変更率 δj${\delta }_j$ となる...以上のプロセスを予測したい時点ぶん繰り返すことになります。

この一連のシミュレーションをコードで書くと下記のようになります。トレンド関数にはロジスティックのほうを用いています。

class LogisticTrendEstimator:
    def fit(self, T, S, delta, k=1, C=1, m=0):
        self._T = T
        self._S = S
        self._delta = delta
        self._k = k
        self._C = C
        self._s_freq = len(S) / len(T)
        self._mu_delta = np.abs(delta).mean()
       
        self._y, self._gamma = self._logistic_trend(T, S, delta, k, C, m)
        
    def _logistic_trend(self, T, S, delta, k=1, C=1, m=0):
        a = np.vstack([np.where(S < t, 1, 0) for t in T])
        gamma = np.zeros(S.shape)
        for j in range(0, gamma.shape[0]):
            gamma[j] = (S[j] - m - gamma[:j].sum()) * (1 - ((k + delta[:j].sum()) / (k + delta[:j + 1].sum())))
        y = C / (1 + np.exp(-(k + (a * delta).sum(axis=1)) * (T - (m + (a * gamma).sum(axis=1)))))
        return y, gamma
    
    def forecast(self, length=10, seed=None):
        np.random.seed(seed=seed)
        
        # generate future change point, and its change rate
        occurrence = np.random.binomial(n=1, p=self._s_freq, size=length)
        generated_s = np.where(occurrence  == 1)[0] + self._T.max()
        generated_delta = np.random.laplace(0, self._mu_delta, generated_s.shape[0])
        
        # predict
        future = np.arange(length) + self._T.max()
        future_y, _ = self._logistic_trend(
            T=future,
            S=generated_s,
            delta=generated_delta,
            k=self._k,
            C=self._C,
            m=self._gamma[-1]
            )
        
        # plot y
        plt.plot(self._T, self._y, c='steelblue', label='past')
        plt.plot(future, future_y, c='darkorange', label='predict')
        plt.xlabel("t", fontsize=16)
        plt.ylabel("y", fontsize=16)
        
        # plot change point
        ymax=np.max([self._y.max(), future_y.max()])
        ymin=np.min([self._y.min(), future_y.min()])
        plt.vlines(
            np.hstack([self._S, generated_s]),
            ymin=ymin, ymax=ymax,
            linestyle='dashed',
            color='gray',
            label='change point'
            )

        plt.legend()
        plt.show()
        return future_y

T = np.arange(100)
S = np.array([20, 60, 80])
delta = np.array([-0.03, 0.01, 0.02])

estimator = LogisticTrendEstimator()
estimator.fit(T=T, S=S, delta=delta, k=0.01, m=0)

pred = estimator.forecast(length=100, seed=123)

2. s(t)$s(t)$: 季節変化

次に季節変化を表現する下記の式を理解していきます。 英語の直訳で「季節変化」と書きましたが、意味的には季節を含め、週、月、年といったあらゆる周期性を s(t)$s(t)$ で扱えます。

季節による変動がある → 周期性がある → 信号処理っぽく表現できる、という発想で s(t)$s(t)$ は下記のように一般的なフーリエ級数で表現されています。

s(t)=∑n=1N(ancos(2πntP)+bnsin(2πntP))${\displaystyle s(t)=\sum _{n=1}^N(a_n\cos \left(\frac{2\pi nt}{P}\right)+b_n\sin \left(\frac{2\pi nt}{P}\right))}$

この式を理解するために、まずフーリエ級数展開の気持ちを簡単に復習します。

フーリエ級数展開について

下記のような曲線をどうにかして関数f(t)$f(t)$で表したいとします。

（これはもちろん私が作ったので事前に知ってるだけですが）調べたら下記の式で表せることがわかりました。

f(t)=3sin(t)+2sin(2t)−sin(3t)${\displaystyle f(t)=3\sin (t)+2\sin (2t)-\sin (3t)}$

どうやらこの曲線は3つの三角関数の和で表現されているようです。3つの三角関数をバラバラにプロットした図をみると下記のようになります。

このようにマクローリン展開などと違って、sinやcosなどの三角関数の和で関数近似しようというのがフーリエ級数展開の特徴です。 季節性の売上など、周期性をもった波形のデータであれば

f(t)=a1sin(t)+a2sin(2t)+a3sin(3t)+....+b1cos(t)+b2cos(2t)+b3cos(3t)+...${\displaystyle f(t)=a_1\sin (t)+a_2\sin (2t)+a_3\sin (3t)+....+b_1\cos (t)+b_2\cos (2t)+b_3\cos (3t)+...}$

というように各周波数(ここではt, 2t, 3t...のこと)の成分を追加していけばどのような波形でも表現可能になります。

ここまでの話をまとめるため、少し強引に一般化した式に直すと

f(t)=∑n=−∞∞(ancos(nt)+bnsin(nt))${\displaystyle f(t)=\sum _{n=-\mathrm{\infty }}^{\mathrm{\infty }}(a_n\cos (nt)+b_n\sin (nt))}$

のようにsin波とcos波の和で表現できます。

三角関数の周波数について

Prophetの季節関数を理解するためにあと一点だけ、周波数nt$nt$の部分について深掘っていきます。 周波数はsin波cos波の振幅数を表しています。2π$2\pi$という単位で一周するので、例えば30日のうち1週間ごとに一周するsin波を表現したい場合は

t=1,2,3,4...,30f(t)=sin(2πt7) $$\begin{gathered} {\displaystyle t=1,2,3,\mathrm{4...},30 \\ f(t)=\sin \left(\frac{2\pi t}{7}\right)} \end{gathered}$$

という式になり、周波数の分母で周期の単位（ここでは１週間なので7）を指定します。 これをグラフ化すると下記のようになり、30日の間にsin波が4周している（30日とは4週間ちょっとの期間なので）ことが見てとれます。

Prophetの季節関数では、このような週、月、年単位の周期を持つことを柔軟に表現できるよう変数 P$P$ を用いて

f(t)=∑n=−∞∞(ancos(2πntP) +bnsin(2πntP) )${\displaystyle f(t)=\sum _{n=-\mathrm{\infty }}^{\mathrm{\infty }}(a_n\cos \left(\frac{2\pi nt}{P}\right)+b_n\sin \left(\frac{2\pi nt}{P}\right))}$

という変形をしています。さらに n$n$ が無限大だと計算量が過多になる & 正の実数のみで十分に季節変化を表現できるため、

∑∞n=−∞$\sum _{n=-\mathrm{\infty }}^{\mathrm{\infty }}$ を ∑Nn=1$\sum _{n=1}^N$ に直し、f(t)→s(t)$f(t)\to s(t)$ に書き換えると

s(t)=∑n=1N(ancos(2πntP)+bnsin(2πntP))${\displaystyle s(t)=\sum _{n=1}^N(a_n\cos \left(\frac{2\pi nt}{P}\right)+b_n\sin \left(\frac{2\pi nt}{P}\right))}$

と変形でき、最初に示した季節変化の関数s(t)$s(t)$が得られます。

パラメタ推定しやすい形に変形する

このとき最適化すべきは [a1,b1,....aN,bN$[a_1,b_1,....a_N,b_N$] の部分、合計2N個のパラメタなので扱いやすいベクトルで抜き出して表現します。

β=[a1,b1,...,aN,bN]T${\displaystyle \mathit{\beta }=[a_1,b_1,...,a_N,b_N{]}^{\mathrm{T}}}$

残りの三角関数の部分もベクトル化して抜き出します。仮にN=10$N=10$の粒度までで関数近似し、年単位の季節変化に見るために P=365.25$P=365.25$とした場合には、

X(t)=[cos(2π(1)t365.25),sin(2π(1)t365.25),...,cos(2π(10)t365.25),sin(2π(10)t365.25)]$$\begin{gathered} {\displaystyle X(t)=[\cos \left(\frac{2\pi (1)t}{365.25}\right),\sin \left(\frac{2\pi (1)t}{365.25}\right),...,\cos \left(\frac{2\pi (10)t}{365.25}\right),\sin \left(\frac{2\pi (10)t}{365.25}\right)] \\ } \end{gathered}$$

というベクトルを作ります。
※ 閏年を含めると1年の平均日数は365.25

これらを用いて結局 s(t)$s(t)$ は下記のように変形できます。

s(t)=X(t)β${\displaystyle s(t)=X(t)\mathit{\beta }}$

論文ではさらに β∼Normal(0,σ2)$\mathit{\beta }\sim \mathrm{N}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{l}(0,{\sigma }^2)$ とし、フーリエ級数の各係数に正規分布を仮定しているようです。

ちなみにパラメタN$N$は多いほど周期変動にうまく当てはまるようになりますが、同時にオーバーフィッティングしてしまう問題も抱えています。論文では、年単位の周期ならN=10$N=10$、週単位ならN=3$N=3$くらいで程よくフィッティングすると書いてあります（個人的にN = 10 はやや多いのでは、という気もしますが）。

3. h(t)$h(t)$: 休日効果

最後に突発的なイベント効果をモデルに組み込む h(t)$h(t)$ について見ていきます。 これも英語の直訳で「休日効果」と書きましたが、意味的には休日含むイベント全てを h(t)$h(t)$ で扱えます。

休日やイベントの多くは事前に予見できるわりに周期といったものはなく、季節変化 s(t)$s(t)$ では取り入れにくい要素です。そのため自動検知云々は諦め、Prophetでは分析者自身がイベントカレンダーのリストを作ってモデルに組み込めるように設計しています。

この設計をProphetではどのような数式で表しているのかを見ていきます。
あるイベントをi$i$とし、それに該当する日付を全て含んだベクトル Di$D_i$ を作ります。 たとえば i=クリスマス$i=クリスマス$ なら

Di=[...,1975/12/25,1976/12/25,...,2020/12/25,...]${\displaystyle D_i=[...,1975/12/25,1976/12/25,...,2020/12/25,...]}$

というように過去と未来全ての12月25日を含むことになります。さらに各時点 t$t$ が Di$D_i$に該当するか否かを表すインディケーターとして

Z(t)=[1(t∈D1),...,1(t∈DL)]${\displaystyle Z(t)=[1(t\in D_1),...,1(t\in D_L)]}$

というベクトルを定義します。ある時点 t$t$ が各イベントi$i$ に該当するかのフラグが 0 or 1で入っています。
各イベント i$i$ に対する係数パラメタを κi${\kappa }_i$ とし、そのベクトルを κ$\mathit{\kappa }$ で表すと、最終的に休日効果 h(t)$h(t)$は

h(t)=Z(t)κ${\displaystyle h(t)=Z(t)\mathit{\kappa }}$

となり、季節変化 s(t)$s(t)$ と同様にパラメタ部分だけをベクトルに分離して表されています。 季節変化のパラメタ β$\mathit{\beta }$ と同様に κ$\mathit{\kappa }$ も下記のように正規分布が仮定されています。

κ∼Normal(0,ν2)${\displaystyle \mathit{\kappa }\sim \mathrm{N}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{l}(0,{\nu }^2)}$

まとめ

冒頭のProphetのモデル式を再掲すると下記のような、主に3つのコンポーネントからできていました。

y(t)=g(t)+s(t)+h(t)+εt.${\displaystyle y(t)=g(t)+s(t)+h(t)+\epsilon t.}$

各要素は最終的に下記のように展開できました。

• g(t)=C(t)1+exp(−(k+a(t)Tδ)(t−(m+a(t)Tγ))(Nonlinear Growth)$g(t)={\displaystyle \frac{C(t)}{1+\exp (-(k+\mathbf{a}(t{)}^{\mathrm{T}}\mathit{\delta })(t-(m+\mathbf{a}(t{)}^{\mathrm{T}}\mathit{\gamma }))}\mathtt{\text{(Nonlinear Growth)}}}$
• g(t)=(k+a(t)Tδ)t+(m+a(t)Tγ)(Linear Growth)$g(t)=(k+\mathbf{a}(t{)}^{\mathrm{T}}\mathit{\delta })t+(m+\mathbf{a}(t{)}^{\mathrm{T}}\mathit{\gamma })\mathtt{\text{(Linear Growth)}}$
• s(t)=X(t)β$s(t)=X(t)\mathit{\beta }$
• h(t)=Z(t)κ$h(t)=Z(t)\mathit{\kappa }$

...本当はこれら未知のパラメタの最適化について書かないとキリが悪いのですが、もう体力が記事のボリュームが限界なので簡単に述べます。

諸々のパラメタを定義していく中で、記事中では下記の3つについてわざわざ確率分布を仮定していました。

δj∼Laplace(0,τ)β∼Normal(0,σ2)κ∼Normal(0,ν2)$$\begin{gathered} {\displaystyle {\delta }_j\sim \mathrm{L}\mathrm{a}\mathrm{p}\mathrm{l}\mathrm{a}\mathrm{c}\mathrm{e}(0,\tau ) \\ \mathit{\beta }\sim \mathrm{N}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{l}(0,{\sigma }^2) \\ \mathit{\kappa }\sim \mathrm{N}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{l}(0,{\nu }^2)} \end{gathered}$$

パラメタごとに分布を仮定しておくと状態空間モデルとして扱うことができます。実際にProphetではこのモデル式がStanで記述され、L-BFGS法などで最適化されているようです。

最後に

私自身が数学があまり得意でないのでかなり噛み砕いて書いてみました。どなたかの理解のとっかかりになれば幸いです。

明日のアドベントカレンダーは Owners Growthチームの大木さんと、PAY株式会社のセキュリティエンジニア、クリスさんです！
お楽しみに！