ill-identified diary

所属組織の見解などとは一切関係ありません

[R] fb Prophet の解剖で学ぶベイズ時系列モデリング

初めに

昨年, KFAS, bsts と, いくつか R の時系列モデリングパッケージを紹介記事を書いた. FaceBook によって開発されたという prophet パッケージも紹介したかったところだが, 日本語での説明は既に公開されている hoxo_m 氏のものが網羅的であり, 使い方の解説としてはこれ以上やることがほぼないと言っていい.

あとはあるとすれば紹介論文やヘルプの全訳くらいだが, そんな面倒 (かつ退屈) なことはしたくない.

そこで, prophet のソースコードを読むことで prophet ひいてはベイズ時系列モデリングについての理解を深めようという変化球な名目で書くことにする.

大きな変化はないと思うが, 検証するバージョンは 0.2.1 である. この記事を読むには R と統計モデリングまたは機械学習の知識も多少必要である.

prophet のモデルの説明

ソースコードだけを見て完全に理解するのは大変なので, モデルの基本的なことだけを説明しておく. Prophet のしくみについて説明したTaylor & Letham (2017, Forecasting at scale) に基づいた説明であり, 上で挙げた Prophet入門【理論編】Facebookの時系列予測ツール でも同様の説明をしているので, これらを見ているならこのセクションは読み飛ばしてかまわない.

prophet では, 時系列 \{y_{t}\} を次のように g(t),s(t),h(t),\varepsilon_{t} という4つの要素に分解して表現する.

\begin{align}
y_{t}= & g(t)+h(t)+s(t)+\varepsilon_{t}\end{align}

g(t)growth つまり成長関数で, 恒常的な上昇/下降トレンドを意味する. 線形トレンドは ARIMA などでも合わせて使われるが, Prophet では線形トレンドだけでなく (つい最近話題になった) ロジスティック曲線を使用できる. 今回読んだのは線形トレンドのほうで, もう一方はロジスティック曲線を利用したものになる. ここまでの話が理解できていれば簡単に解読できると思うので, ロジスティックトレンドバージョンの説明は省略する. 線形トレンドの場合, 式は, 成長率をk とすれば, kt という風に表せる. ここに変化点の要因を加える.

\begin{align}
a_{j}(t)= & \begin{cases}
1 & \text{if }t\geq s_{j}\\
0 & \text{otherwise}
\end{cases}\end{align}

というダミー変数を作ると, s_{j} 時点で変化が発生したということを表現できる. S 個の変化点があるとして, \{a_{j}(t)\}_{j=1}^{S}をまとめて \boldsymbol{a}_{t}(t) というベクトルができる. 変化点での変化の大きさを, \begin{gathered}S\end{gathered}
個のパラメータ\boldsymbol{\delta} で表せば,

\begin{align}
g(t):= & kt+\mathbf{a}(t)^{\intercal}\boldsymbol{\delta}t\\
= & (k+\boldsymbol{a}(t)^{\intercal}\boldsymbol{\delta})t\end{align}

と書き換えられる. しかし, ここでは変化点まわりでの平滑化のために, さらに切片パラメータ m\boldsymbol{\gamma}というパラメータを追加して,
\begin{align}
g(t):= & (k+\boldsymbol{a}(t)^{\intercal}\boldsymbol{\delta})t+(m+\boldsymbol{a}(t)^{\intercal}\boldsymbol{\gamma})\end{align}
としている*1.

ロジスティックトレンドなら,

\begin{align}
g(t):= & \frac{C(t)}{1+\exp\left( (k+\mathbf{a}^{\intercal}\boldsymbol{\delta})(t-m-\boldsymbol{a}(t)^{\intercal}\boldsymbol{\gamma})\right)}\end{align}

となる. C(t) はロジスティックトレンドの上限を表す値を意味し, C(t)=1 なら上限は1になる. C(t) となっているように, 各時点ごとに異なる上限値を設定できる.

t の係数が複雑なのは, ここに変化点の要因を含んでいるからである.

s(t) は seasonalily, 季節変化を表す項である. という用語を訳しただけで, 四季だけでなく毎月とかもっと短い変化も設定できる. 「周期性変化 (periodic changes)」と言ったほうがわかりやすい. Prophet では, 周期性変化を表現するために, s(t)フーリエ級数と定義している. N 次のフーリエ級数には2N 個のパラメータが必要で, そのパラメータが \beta である. 周期を P とすると,

\begin{align}
s(t):= & \begin{bmatrix}\cos(2(1)\pi t/P)\\
\sin(2(1)\pi/P)\\
\vdots\\
\cos(2(N)\pi/P)\\
\sin(2(N)\pi/P)
\end{bmatrix}^{\intercal}\begin{bmatrix}\beta_{1}\\
\vdots\\
\beta_{2N}
\end{bmatrix}\\
= & \boldsymbol{x}(t)^{\intercal}\boldsymbol{\beta}\end{align}

と表せる. 周期 P を事前に設定*2すれば,\boldsymbol{x}(t) を事前に計算できるので, \boldsymbol{x}(t) はパラメータではなく, データとして扱える. よって, s(t) の部分は線形回帰モデルとして扱える.

h(t) は holiday, 休日効果である. 祝日など, 季節性とは別に, 特定のタイミングで一時的に発生することがわかっている効果を表す. 祝日のリストを \{D_{1},D_{2},\cdots,D_{L}\} とすると,

\begin{align}
h(t):= & 1(t\in D_{1})\kappa_{1}+\cdots1(t\in D_{L})\kappa_{L}\\
= & \begin{bmatrix}1(t\in D_{1})\\
\vdots\\
1(t\in D_{L})
\end{bmatrix}^{\intercal}\begin{bmatrix}\kappa_{1}\\
\vdots\\
\kappa_{L}
\end{bmatrix}\\
= & \boldsymbol{d}(t)^{\intercal}\boldsymbol{\kappa}\end{align}

と表せる.

この\boldsymbol{z}(t)は, 単純にt に応じてダミー変数を与えるだけで, 事前に計算できる. よって, ここも線形モデルとして扱える.

最後の \varepsilon_{t} は誤差項である. 平均ゼロの独立正規分布が仮定される.

というわけで, 実はモデルのほとんどが線形である. 特に g(t) を線形トレンドにした場合は,

\begin{align}
y_{t}= & k+\boldsymbol{a}(t)^{\intercal}t+\boldsymbol{x}(t)^{\intercal}\boldsymbol{\beta}+\boldsymbol{z}(t)^{\intercal}\boldsymbol{\kappa}+\varepsilon_{t}\end{align}

となり, 事実上線形回帰モデルと同じである. 加えて, 0.2 では他の説明変数を追加する機能も登場した. ただし, 線形回帰しかできない.

ハリボテの R言語

では, ソースコードを覗いてみよう. 本当は github のものを見るのが良いのだが, 今回は prophet パッケージを読み込んだ状態で, prophet と入力して, チュートリアルで使われている prophet() を確認する. この関数だけで, モデルのパラメータ推定ができるらしい. 覗いてみると, この関数は入力された引数を整形して, fit.prophet() に渡しているだけだとわかる. 名前からして, fit.prophet() はパラメータ推定の処理をしているのだろう*3. では, fit.prophet() はどのような処理をしているのだろうか.

fit.prophet() 関数の内容はは少し長いが, 最初のほうは例外判定と, 入力データの整形の処理である. 結局, 重要なのは, 以下の部分である.

  if (min(history$y) == max(history$y)) {
    m$params <- stan_init()
    m$params$sigma_obs <- 0
    n.iteration <- 1
  }
    stan.fit <- rstan::sampling(model, data = dat, init = stan_init,
                                 iter = m$mcmc.samples, ...)
    m$params <- rstan::extract(stan.fit)
    n.iteration <- length(m$params$k)
  }
  else {
    stan.fit <- rstan::optimizing(model, data = dat, init = stan_init,
                                   iter = 10000, as_vector = FALSE, ...)
    m$params <- stan.fit$par
    n.iteration <- 1
  }

ここでは rstan::sampling()rstan::optimizing() という, rstan パッケージの関数を呼び出している.

実は, prophet パッケージは, モデルのパラメータ推定を R でやっていない. この部分はベイズ統計モデル用パッケージ, rstan に完全に依存している*4. グラフはさすがに R の ggplot2 を使わざるを得ないが, R 言語で書かれた部分は, データフレームを stan の入出力のために変換する処理や, 推定したパラメータをもとに予測値を計算する処理だけで, 事実上 API としての機能がほとんどである. これは Python でも同様である*5.

というわけで, rstan でどのような処理をしているかを見るために, stan の言語で書かれたソースコードを読み込む過程を, 順を追って説明する*6. prophet_linear_growth.stan, prophet_linear_logistic.stan という2つのファイルがあるが, まず前者を見てみる*7.

stan ソースコードは, 主に data{...}, parameter {...}, model{...} ブロックで構成される. それぞれ順に, 計算に使うデータ (=計算を通して変わらない値) の変数名と型を宣言するブロック, パラメータに使う変数と型を宣言するブロックで, そして事前分布や尤度を記述し, モデルを指定するブロックである. data ブロックは, 以下のようになっている.

data {
  int T;                                // Sample size
  int<lower=1> K;                       // Number of seasonal vectors
  vector[T] t;                            // Day
  vector[T] y;                            // Time-series
  int S;                                // Number of changepoints
  matrix[T, S] A;                   // Split indicators
  real t_change[S];                 // Index of changepoints
  matrix[T,K] X;                // season vectors
  vector[K] sigmas;              // scale on seasonality prior
  real<lower=0> tau;                  // scale on changepoints prior
}

R と違い, stan では変数を作成するときに型を指定する. 整数なら int, 実数 (小数) なら real, (列) ベクトルなら vector, 行列なら matrix である. 例えば最初の int T; というのは, 時系列の長さを指定している. 実際の値は, rstan::sampling() などから与えられたものが代入される. vector[ ]matrix[ , ] などの型名の後の [ ] は次元を指定してる. また, 最後の real tau; というのは, 下限がゼロの実数型変数 tau という意味である. コメントによると, 変化点パラメータのスケールらしいので, 正値に限定する必要があるようだ. parameter ではなく data ブロックにあるのは, tau は最初から任意の設定するため, ここでは推定する必要がないからである.

次に, parameter を確認する. 構文は data ブロックとほぼ同じだが, stan では int 型のパラメータを計算することができない.

parameters {
  real k;                            // Base growth rate
  real m;                            // offset
  vector[S] delta;                       // Rate adjustments
  real<lower=0> sigma_obs;               // Observation noise (incl. seasonal variation)
  vector[K] beta;                    // seasonal vector
}

ここから, k, m, delta, sigma_obs, beta という5種類のパラメータがあることがわかる. 親切なことに, こちらも各パラメータの意味もコメントしてあり, 前節での式と概ね一致していることがわかる.

なお, この後さらに transformed parameter というブロックがあるが, これは他のパラメータと連動して値が決まるようなパラメータを定義するブロックである. ここでは成長関数 g(t) のパラメータ \boldsymbol{\gamma}が定義されている.

今すでにわれわれは, 各変数のどれがパラメータでどれがデータか知っているので, あとは model {...} ブロックで各パラメータにどのような事前分布が設定され, パラメータでモデルがどう表現されているかを見れば, モデルのおおまかなことがわかる.

ここで, ベイズ統計についておさらいしておこう. ベイズ統計は, 確率で表現したモデルの対数 (=対数尤度) \ln f(\{y_{t}\},\theta)と, パラメータに対して設定した事前確率の対数 \ln p(\theta) の和である, 対数事後確率 \ln p(\theta\mid\{y_{t}\}) を評価してパラメータを決める:

\begin{align}
\ln p(\theta\mid\{y_{t}\})= & \ln f(\{y_{t}\},\theta)+\ln p(\theta)\end{align}

一方で, 最尤推定の場合は, 事前分布を使わず, 対数尤度 \ln f(\{y_{t}\};\theta) だけを最大化するパラメータを求める.

以下が model ブロックの記述である.

model {
 //priors
  k ~ normal(0, 5);
  m ~ normal(0, 5);
  delta ~ double_exponential(0, tau); // = Laplace
  sigma_obs ~ normal(0, 0.1);
  beta ~ normal(0, sigmas);
  // Likelihood
  y ~ normal((k + A * delta) .* t + (m + A * gamma) + X * beta, sigma_obs);
}

stan をある程度使ったことのある人間なら, 単純すぎて拍子抜けするかもしれない. あれだけ自由度のある時系列モデルだから, さぞかし複雑なモデルを設定しているのだろう, と思ったらこれである. ◯ii 威力棒を分解した画像を思い越させる.

model ブロックで使う x ~ ... ; という構文は, X\sim\mathcal{N}(\mu,\sigma^{2}) のような数式と同じで, 事後分布に記述した確率分布を足すということである (実際には確率の対数値で計算する). そのため, 構文の上では, 尤度も事前分布も ~ を使って記述できてしまう*8. しかし, 既にどの変数がパラメータなのか我々はわかっているし, ここでも親切なことにコメントが書かれている. //prior と書かれた行以下が事前分布の設定で, // Likelihood 以下が尤度である.

よって, 最初の k ~ normal(0, 5); は, 事前分布の対数確率,

\ln f(k;\mu=0,\sigma=5)

を足している.

次に //Likelihood 以下に注目する. y ~ normal(...) ; と書かれているので, 正規分布だとわかる. 変数 y はパラメータではなくデータだから, コメントに書いてあるとおり, これは事前分布ではなくモデルの記述である. ここから,

\begin{align}
y_{t} & \sim\mathcal{N}(\cdots,\sigma^{2})\end{align}

のような正規分布の時系列モデルだと見当がつく.

Stan は R のようにベクトルを意識した書き方ができるので, ynormal(...) 内のパラメータがベクトル型であれば, 上記の構文で y_{1},y_{2},\cdots,y_{T} に対して, 正規分布の確率が計算され, 合計して対数尤度に加算される. つまり, y ~ normal(...) ;
\begin{align}
\ln f(\{y_{t}\},\theta)= & \sum_{t=1}^{T}\ln p(y_{t};\cdots,\sigma^{2})\end{align}
を対数事後確率に足すという意味である. R 風に書くと,

log_p <- dnorm(y, mean=..., sd=sigma, log=T)
log_likelihood <- sum(log_p)

としていることになる.

では, この正規分布の平均パラメータ内は何だろうか.

(k + A * delta) .* t + (m + A * gamma) + X * beta

は, パラメータのコメントから, (k + A * delta) .* t + (m + A * gamma) の部分は成長トレンド項 g(t) に相当するとわかる. kk, A\boldsymbol{a}(t) を積み上げた行列, delta\boldsymbol{\delta} に対応する. R と違い, 行列やベクトルどうしの積は単に * と書く. 一方で, t は行列の積としてではなく, 要素毎の積 (アダマール積) を計算したいため, .* と書く*9. よって,

\left[k+\begin{bmatrix}\boldsymbol{a}(1)^{\intercal}\\
\vdots\\
\boldsymbol{a}(T)^{\intercal}
\end{bmatrix}\boldsymbol{\delta}\right]\circ\begin{bmatrix}1\\
\vdots\\
T
\end{bmatrix}+m+\begin{bmatrix}\boldsymbol{a}(1)^{\intercal}\\
\vdots\\
\boldsymbol{a}(T)^{\intercal}
\end{bmatrix}\boldsymbol{\gamma}

となる. 残る X * beta は, 変数名から周期性トレンド項 s(t) と考えられる. では, 休日効果項 h(t) はどこだろうか. これは, 少しわかりにくいのだが, s(t)h(t) も線型結合なので, 一つの行列に集約している. この処理は R 側の fit.prophet() 関数でなされている.

\begin{bmatrix}\cos\frac{2\pi t}{P}\\
\vdots\\
\sin\frac{2(N)\pi t}{P}\\
1(t\in D_{1})\\
\vdots\\
1(t\in D_{L})
\end{bmatrix}^{\intercal}\begin{bmatrix}\beta_{1}\\
\vdots\\
\beta_{2N}\\
\kappa_{1}\\
\vdots\\
\kappa_{L}
\end{bmatrix}

よって, X * betas(t)+h(t) に相当する*10.

prophet_linear_logistic.stan は, 成長トレンド項 g(t) がロジスティック曲線の場合である. 内容は 尤度の記述いがいはほとんど変わらないため, 説明を省略する.

なぜベイズ推定する必要があるのか

Prophet のモデルは, 一般化加法モデル (GAM) の部類に入る. この種のモデルは単純なので, (準) ニュートン法のような古典的な最適化アルゴリズムでも最適パラメータを計算できてしまう. にも関わらずなぜ, 尤度だけを使った最尤推定ではなく, 事前分布を設定し, stan を使ったベイズ推定をするのか. というわけで, prophet の事前分布は決して「雰囲気でやっている」ものではないことを説明する.

重要なポイントは「変化点の自動検知」という機能である. Prophet では変化点の発生回数上限さえ設定すれば, 変化点の数と位置を自動で決めてくれる. しかし, このようにするのは変化点がどこで起こったのかはわからないからで, そのため変化点の上限回数を多めに指定することになる. そうなると, 事前分布を使わずに最尤法で推定すれば, 推定されたモデルは過剰適合を引き起こす可能性がある.

そこで, 過剰適合しないように事前分布を設定している. ソースコードを見ると,

delta ~ double_exponential(0, tau); // = Laplace

とある. 二重指数分布とはなにか, と思うが, 親切にもラプラス分布のことだとコメントされている. ラプラス分布は平均パラメータ \mu とスケールパラメータ \sigma を持ち, その密度関数は,

\begin{align}
p(\delta;\mu,\sigma):= & \frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{\left\Vert \delta-\mu\right\Vert _{1}}{2\sigma}\right)\end{align}

と定義される. 参考までに正規分布と比較した図 1をあげてみる*11. このような尖った確率分布は, 正規分布と比べて, ゼロ, もしくは絶対値の大きな値がでやすいことがわかる.


f:id:ill-identified:20180528013747p:plain

図1: ラプラス分布とガウス (正規) 分布


別の視点からラプラス事前分布の意味を見てみる. \mu=0,\sigma=\tau のときに対数をとると,

\begin{align}
\ln p(\delta;\mu,\sigma)= & -\ln\sqrt{2\pi}\tau-\frac{1}{2\tau}\left\Vert \delta\right\Vert _{1}\end{align}

である. ベイズ推定における目的関数である対数事後確率は, 対数尤度と, この対数事前確率の和となる.

機械学習の知識が多少ある人なら \frac{\tau}{2}\left\Vert \delta\right\Vert _{1} に見覚えがあるだろう. 1/2\tau正則化パラメータだとすると, これはパラメータ \delta に対して L_{1}正則化をしていることと同じである. つまり, ラプラス分布を事前分布にするということは, L_{1}スパース推定と同じ効果があり, パラメータ \delta がゼロになるように制約がかかる.

というわけで, ラプラス事前分布は, ベイズ統計モデル的に見れば, 変化点の数をあらかじめ特定しないノンパラメトリックなモデルとして説明できるし, 機械学習の文脈で見れば, スパース学習として説明できる.

残りのパラメータについては, どれも平均ゼロの正規事前分布を与えている. 正規分布の密度関数がいずれも

\begin{align}
p(\theta;0,\sigma)= & \frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{\left\Vert \theta\right\Vert _{2}^{2}}{2\sigma^{2}}\right)\end{align}

であることを思い出せば, ラプラス分布と同様に対数をとると,

\frac{1}{2\sigma^{2}}\left\Vert \theta\right\Vert _{2}^{2}

が現れるとわかる. つまり, 正規事前分布は L_{2} 正則化と同じ性質であることがわかる. Prophet では, 成長関数の基準値 k, 切片 m の事前分布は分散パラメータが固定されているが, s(t) の周期性効果と h(t) の休日効果の係数パラメータの正規事前分布については, それぞれ個別に異なる分散パラメータを与えられるので, 両者の効果のスケールを差別化したモデリングができる*12.

ところで, やや蛇足なことを言うが, 「最適化問題なら古典的な統計モデル, モンテカルロ法ならベイズ」という区分は正しくない. 実際, prophet は, デフォルトでは rstan::optimizing() を呼び出して計算している. この関数は, 事後確率を最大化するパラメータを導く, MAP 推定をおこない, 実際に計算に用いている最適化アルゴリズム準ニュートン法である. よって, prophet のモデルは, モンテカルロ法を使わないが, 事前分布を使うという意味ではベイズ統計モデルということになる.

結論

モデルというと, つい複雑なものを作ってしまうが, この prophet で採用されているモデルはとてもシンプルで, それでいて表現の自由度がある. ソースコードも stan を使うことでシンプルに記述されている. ユーザー側に求められる操作も, 時系列データの入力, 祝日やイベント発生日情報の入力, 成長トレンドの指定くらいであり, 「専門知識のないユーザでも使えるスケーラブルなパッケージ」という謳い文句に偽りはない.

参考文献



Taylor, S. J., & Letham, B. (2017). Forecasting at scale. PeerJ Preprints, 5, e3190v2. doi:Forecasting at scale [PeerJ Preprints]


*1:詳しい意味は最初に挙げた参考文献をあたること.

*2:ただし, 年周期は 365.25 に設定されているため, うるう年効果がないことに注意する.

*3:厳密には fit=TRUE のときだけ呼び出される

*4:ここでは, rstan の構文は必要最小限の説明で済ます. 詳しく知りたい場合は公式ドキュメントや, 松浦健太郎『StanとRでベイズ統計モデリング』(通称, 松浦本, アヒル本) などを参照せよ. 統計モデリングについても初心者という人は, 後者が統計モデリング全般の理解についても役立つ説明をおこなっているので推奨する.

*5:Python の場合, stan ソースコードWindows OS 用と UNIX OS 用で分かれているが, 本質的な違いはない. OS ごとに依存しているライブラリが違うため, 行列オブジェクトの型が異なっているというだけである.

*6:Prophet はドキュメントが充実してるので, それを先に読んだほうが早く理解できる. しかしここではプログラムから解き明かすという体裁にしたい.

*7:RStudio には stan の構文をハイライトしてくれる機能があるので, RStudio で開くと見やすい.

*8:最近では, わかりにくいからか target += ...; という記述が推奨されている. target は予約変数で, 対数事後確率にこの値が足される.

*9:つまり, matlab と同じである.

*10:他の説明変数を追加する場合は, それもこの項に統合される.

*11:線形回帰モデルが正規分布の場合, 最小二乗法となるのに対して, ラプラス分布にすると, ラプラスが最初に考案したとされる最小絶対誤差法になる. これがラプラス分布という別名の理由である.

*12:しかし, stan のソースコードを見ると, seasonalities と holidays をひとまとめで扱っているので, スケールを同一にするしかなくなる, ということに気づく人もいるかもしれない. これについては, stan にデータを与える前に, データの方のスケールを調整することで表現している.