はじめに

ベイズ統計について勉強していたのですが、実際に使いこなせるよう実装力をつける意味でPyMCを使い始めようとしたところ、はじめの「Introductory Overview of PyMC」内の以下の例で詰まりました。

ケーススタディ 1: 聴覚障害のある子どもたちの教育成果

Introductoryと書いておきながら、聞いたことのない分布や数理モデルなどに思わず面を食らってしまいました。全体を理解するにあたり様々な記事を参考にさせていただいたのですが、参照した記事に加え、通常の馬蹄分布の改良版であったり、潜在パラメータを使った実装であったりと、調べてもぱっと出てこないような部分もあり本当に詰まるところが多かったです。（よく見ると参考文献にしっかり出典がありました。ライブラリのイントロダクションということもあり、平易なはずだと思って読んだいたため気づくのに時間がかかってしまいました。）

以下で紹介しますが、様々な知識を必要としますので、一つも知らないという場合は理解するのが少々骨折りかと思います。私はスパースベイズ推定や馬蹄分布について何の知識もなかったため、とても複雑なことをたくさんしているのでは？とかなり構えていました。しかし、数式的なテクニックとしては大きくまとめると「Horseshoe事前分布（馬蹄分布）を用いたスパースベイズモデル」という一つのテクニックを軸にしているだけなので案外シンプルです。

本記事は、全体像とポイントとなる箇所の紹介に留めています。細かな説明などは参考文献の方で十分に解説されていますので、そちらを参照ください。
まず初めに参考文献の方をざっとご参照ください。参考文献に載っている内容につきましては、重複を避けるためこちらの記事では詳しく解説しておりません。
誰かの参考になれば幸いです。

参考文献

トピック

ベイズモデル
再パラメータ化
階層的回帰モデリング
スパースモデリング
馬蹄分布(Horseshoe Distribution)
・正則化馬蹄分布（regularized horseshoe）
縮小係数

ベイズモデル

数理モデルの全体像になります。全体の見通しが良くなるように数式をまとめて以下に書いておきました。

回帰係数$\beta_j$は単純なベイズ線形回帰では、それ自体が何らかの分布に従いますが、階層モデリングにより$\beta_j$が従う分布のパラメータ自体が、多数の他の事前分布に従うパラメータになっています。
また、Horseshoe系のモデルではシュリンケージパラメータとして、全体を一括で押さえつけるグローバルパラメータ（$τ$）と、個々の係数ごとに自由度を与えるローカルパラメータ（$λ$）の2段構えになっているのが特徴です。

(1) 事前分布

\begin{aligned}
&\sigma \sim \mathrm{HalfNormal}(0,\,25),\\
&\tau \sim \mathrm{HalfStudentT}\!\Bigl(2,\;\frac{D_0}{D-D_0} \cdot \frac{\sigma}{\sqrt{N}}\Bigr),\\
&\lambda_j \sim \mathrm{HalfStudentT}(5),\quad j=1,\dots,D,\\\\
&c^2 \sim \mathrm{InverseGamma}(1,\,1),\\\\
&\tilde{\lambda}_j^2 = \frac{c^2 \lambda_j^2}{c^2 + \tau^2 \lambda_j^2} \\ \\
\end{aligned}

(2.1) 回帰係数（理論式）

\begin{aligned}
\beta_j \sim N\bigl(\;0,\ \tau^2 \cdot \tilde{\lambda}_j^2\;\bigr) \qquad \qquad \qquad \qquad \\ \\
\end{aligned}

(2.2) 回帰係数（再パラメータ化を考慮した、実際の実装）

\hspace{4em}
\begin{align*}
\\
z_j &\sim \mathcal{N}(0,\,1),\quad j=1,\dots,D,\\ \\
\beta_j &\;=\; z_i \;\cdot\; \,\tau\, \;\cdot\; \tilde{\lambda}_i \\
& \;=\; z_j \;\cdot\; \tau \;\cdot\; \lambda_j \;\cdot\; \sqrt{\frac{c^2}{c^2 + \tau^2 \lambda_j^2}},\quad j=1,\dots,D,\\
\end{align*}

(3) 切片

\begin{aligned}
&\beta_0 \sim \mathcal{N}(100,\,25^2), \qquad \qquad \qquad \quad \qquad 
\end{aligned}

(4) 観測モデル

\hspace{3em}
\begin{aligned}
&y_i \;\sim\; \mathcal{N}\!\Bigl(\;\beta_0 + \sum_{j=1}^D x_{ij} \beta_j,\;\sigma^2\Bigr),\quad i=1,\dots,N.
\end{aligned}

再パラメータ化（reparameterization）

(2.1)式と(2.2)式についての説明になります。潜在パラメータ$z$の正体です。
ベイズモデルを実装するときに、うまくモデルが収束しない、サンプリングがなかなか進まない、効率的にサンプリングできないといった問題を経験することは少なくないと思います。そうした際に有用なのが「再パラメータ化（reparameterization）」です。
「再パラメータ化」とは、同じ確率モデルを別のパラメータで書き直すことを指します。数学的には同値であっても、サンプリングの効率（MCMCのミキシングや勾配計算の安定性）が大きく変わるケースがあります。
再パラメータ化は同じモデルを「よりサンプリングしやすい形」に書き直すことで、サンプリング効率や収束性、数値安定性などを改善させる重要なテクニックです。「ノンセンタリング (non-centered parameterization)」とも呼ばれています。

階層的回帰モデリング

階層的回帰モデリング、階層型線形回帰などというと、マルチレベルモデル（またの名をHierarchical Linear Model (HLM)）を指すことが一般的なようです。マルチレベルモデルとは、データ自体がグループに分かれておりデータ全体として階層構造を取っているような場合に、数理モデルのパラメータも入れ子状態にし、上位のパラメータで束ねることで、データのクラスや階層構造に関連付けた推定を行います。
しかし、ベイズ統計の文脈で階層的回帰モデリングというと、これに限らないようです。
今回の例のようにデータのグループ分割がないモデルでも、パラメータが多層的に定義されていれば「階層モデル」と呼んでいます。
「階層～」という言葉に対して、厳密な定義はないようなので、あまり気にしても意味が無いように思います。

スパースモデリング

その名の通り、解がスパースであるという前提があるとき、数理モデルのパラメータに正則化を施すことです。回帰係数の多くがゼロ近傍（不要な特徴量）で、一部のみが有意（重要な特徴量）だと想定される状況を「スパースな状況」と呼びます。
このパラメータを抑えつけるような正則化により、「回帰」と「変数選択」を同時に行うことが出来ます。

馬蹄分布(Horseshoe Distribution)

正則化馬蹄分布（regularized horseshoe）

馬蹄分布の新しい拡張版である「正則化されたHorseshoe」を使用しています。これは、元の馬蹄事前分布に以下の特性を追加したものです。
・スパース性を明確に制御できるフレームワークを提供。
・正則化されたHorseshoeにより、大きな係数を適切に制御可能。

馬蹄分布に関しては上記の参考文献に説明がありますが、正則化馬蹄分布に関しては、参考となる記事などが見当たらなかったため、論文を参照しましょう。

縮小係数

直接モデルに登場するパラメータの事前分布というわけではない「縮小係数」なるものの分布を考えることがなぜ重要なのかですが、回帰係数の事後分布平均が、縮小係数の事後分布平均にダイレクトに依存するためだそうです。数式の設計上の動機を理解するためとも言えます。

おわりに

やはりチュートリアルでも、手元で動かしてみると見えてくるものがたくさんあります。
また今回記事を書くにあたり、正直に申し上げますと、重要なポイントはChatGPTがほとんど教えてくれました。本当に何でも知っているなぁと改めて感心させられました。
それにしてもやはり、イントロダクションでこれって…（ムズカシ過ぎないか？）

間違いなどお気づきになりましたら、ご連絡いただけますと幸いです。

続きはこちら