« 読了：Guo, Blundell, Wallach, Heller (2015) ベイジアン・エコー・チェンバー | メイン | 読了：Asparouhov, Hamaker, & Muthen (2017) ものどもひれ伏せよ、これが動的潜在クラス分析だ »
2017年8月18日 (金)
Scott, S.L., Varian, H. (2014) Predicting the Present with Bayesian Structural Time Series. International Journal of Mathematical Modelling and Numerical Optimisation, 5, 4-23.
Googleの中の人謹製、ベイジアン構造時系列モデリングのためのRパッケージbstsの元論文。draftのPDFで読んだ。 別に読みたかないけど、実戦投入前の儀式としてぱらぱらめくった次第。
このモデルでの主な使い道として想定されているのは、Google Trendsのようなたくさんの、そんなに長くない時系列があり、それらを予測子としてある目的変数の時系列をnowcastingしたいんだけど、でも予測子のうちほんとに効く奴はいくつかしかない、という状況である。動的因子分析なんかで多変量時系列を縮約するのではなく、ベイジアンモデル平均のアプローチで事後分布をシミュレートする。
モデルの説明。
時点$t$における観察値$y_t$について、次の状態空間表現を考える。
観察方程式 $y_t = Z^T_t \alpha_t + e_t \ \ e_t \sim M(0, H_t)$
遷移方程式 $\alpha_t = T_t \alpha_t + R_t \eta_t \ \ \eta_t \sim N(0, Q_t)$
さて、
- カルマン・フィルタとは...[略]
- カルマン・スムーザーとは...[略]
- ベイジアン・データ拡大法とは... $\mathbf{y} = y_{1:n}, \mathbf{\alpha} = \alpha_{1:n}$として$p(\mathbf{\alpha}|\mathbf{y})$からうまいことサンプリングしてシミュレーションする方法。Durbin & Koopman (2002, Biometrika)をみよ。
外的な予測子の効果$\beta^T \mathbf{x}$について。$\beta$は時変なしとし、$\alpha_t$に値1の変数を追加し、 $\beta^T \mathbf{x}$を$Z_t$に入れ込む。事前分布としてspike-and-slab事前分布を使う[メモは後述]。するとパラメータの事後分布がこんな風になる...[略]。推定にはMCMCをこんな風に使う...[略]。
事例はふたつ。(1)週次の失業保険新規申請件数をGoogle Trendsで予測する。(2)月次の全米小売売上金額をGoogle Trendsで予測する。
この方法で擬似相関が回避できるわけではないんだけど(実際、後者の例では変なキーワードの検索量時系列が売上金額の強力な予測子となっている)、分析者の主観的判断を事前分布として生かせるという点が特徴。
... bstsパッケージというのはMARSSパッケージのMCMC版みたいなものかと思ってたんだけど、蓋をあけてみたらベイジアンモデル平均の話であった。なんでも読んでみるもんね。
マニュアルをちらっとみたところでは、dlmパッケージのように要素別ヘルパー関数を組み合わせて使うようで、なんだか面白そうだ。と、自分で自分に景気をつけて...
以下、spike-and-slab事前分布についてのメモ。ついつい逐語訳になってしまった。
ベイジアン・パラダイムで[たくさん予測子があるけど効くのはわずかという]スパース性を表現する自然な方法として、回帰係数にspike and slab事前分布を与えるという方法がある。
$\gamma_k$を、$\beta_k \neq 0$のときに$1$, $\beta_k=0$のときに$0$となる変数とする。$\beta$の非ゼロ要素からなる下位集合を$\beta_\gamma$とする。
spike-and-slab事前分布とは
$p(\beta, \gamma, \sigma^2_e = p(\beta_\gamma|\gamma, \sigma^2_e) p(\sigma^2_e |\gamma) p(\gamma)$$p(\gamma)$の周辺分布は「スパイク」である。つまり、それはゼロの位置に正の確率質量を持つ。原理的には、$p(\gamma)$を調整して階層原理のようなベスト・プラクティスを実現できる(階層原理とは、高次の交互作用項が入るときは低次の項も入るという原理)。実際には、単に独立ベルヌーイ事前分布を使うのが便利である [各予測子が独立に確率$\pi_k$でモデルに入るってことね] 。
$\gamma \sim \sum_{k=1} \pi^{\gamma_k}_k (1-\pi_k)^{1-\gamma_k}$
もっと単純化して、すべての$\pi_k$を$\pi$とすることも多い。これは$\pi_k$をいちいちセットするのは大変だという主旨なのだけど、事前分布の交換可能性という基盤があれば正当化できる。$\pi$を決める自然な方法のひとつとして、分析者に「モデルの 期待されるサイズ」を尋ねるという手がある。分析者が「非ゼロの予測子は$p$個」というならば、$\mathbf{x_t}$の次元数を$K$として$\pi = p / K$とすればよい。場合によっては、特定の予測子について$\pi_k$を0か1に決め打ちするのも便利である。(我々は採用しないけど)別の戦略として、予測子を「モデルに入りそうか」でいくつかのグループに主観的に分け、それぞれのグループの予測子に主観的に決めた確率を与えるという手もある。正方行列$\Omega^{-1}$について、その行と列が$\gamma_k=1$に対応している行列を$\Omega^{-1}_\gamma$とする[←えええ...? 先生すいません、頭悪くてよくわかんないっす...]。条件つき事前分布$p(1/\sigma^2_e | \gamma)$と$p(\beta_\gamma | \sigma_e, \gamma)$は次の条件つき共役対として表現できる。
$\beta_\gamma | \sigma^2_e, \gamma \sim N(b_\gamma, \sigma^2_e (\Omega^{-1}_\gamma)^{-1})$
$\frac{1}{\sigma^2_e} | \gamma \sim Ga(\frac{\nu}{2}, \frac{ss}{2})$
ここで$Ga(r,s)$は平均$r/s$, 分散$r/s^2$のガンマ分布である。これが「スラブ」[厚板]である。なぜそう呼ぶかというと、すごく弱い情報性しか持たないように(一様に近づくように)パラメータを選べるからである。$\gamma$と同様、ここでも合理的なデフォルト値が存在する。まず、事前平均ベクトル$b$はゼロにするのが一般的である。ただし、特定の予測子が特に有用だという信念がある場合には、情報的な事前分布を指定することも容易である。$ss$と$\nu$は事前の平方和と事前のサンプルサイズだと解釈できる。これらはユーザに、回帰に期待する$R^2$と、その推測がどのくらいのサイズの観察と同じ重みをもつと思うか($\nu$)を尋ねてセットすればよい。ここで$ss/\nu = (1-R^2) s^2_y$である。$s^2_y$とは反応の周辺SDである。$s^2_y$でスケーリングするのはベイジアンのパラダイムからちょっと逸脱しているのだが(事前分布をデータで決めているわけだから)、便利な方法だし、実務的な悪影響はない。上の式でもっとも高次元なパラメータは、モデル全体の事前情報行列$\Omega^{-1}$である。$\mathbf{x}_t$を$t$行目に積みあげて得られる計画行列を$\mathbf{X}$としよう。通常の回帰モデルの尤度は情報行列$\mathbf{X}^T \mathbf{X} / \sigma^2_e$を持つ。だから、
$\Omega^{-1} = \kappa \mathbf{X}^T \mathbf{X} /n$
とすれば、事前平均$b$に観察数$\kappa$ぶんの重みを与えたことになる。これがZellnerのg事前分布である [←ぐああああ... わっかんないよう...]。実務的には、$\mathbf{X}$の列の間の完全な共線性に対する防御策が必要となる。フルランクを保証するためには、
$\Omega^{-1} = \kappa (w \mathbf{X}^T \mathbf{X} + (1-w) diag( \mathbf{X}^T \mathbf{X} ) )/n$
とするとよろしい。我々はデフォルト値として$w=1/2, \kappa=1$を採用している。
[ふひいいいいい。こんな話を理解するために生まれたわけではないわ！ とにかく情報行列についてうまいこと考えて下さったってことですね、はい、信じます先生]まとめよう。spike-and-slab事前分布を使うことで、事前パラメータ$\pi_k, b, \Omega^{-1}, \nu$を通じて事前の意見を簡単に表現することができる。 ある種の合理的想定に頼るというコストを払って単純さを選ぶ分析者は、事前情報を、期待されるモデルサイズ、期待される$R^2$、その推測に与える重みを表すサンプルサイズ$\nu$にまで切り詰めることができる。事前分布について考えること自体を避けたいという分析者のために、我々のソフトはデフォルト値$R^2=0.5, \nu=0.01, \pi_k=0.5$をご用意している。
論文：データ解析(2015-) - 読了：Scott & Varian (2014) ベイジアン構造時系列モデリング