Spike-and-Slab事前分布を理解する 〜理論＆PyMCによる実装編〜

この記事は、確率的プログラミング言語 Advent Calendar 2023 24日目の記事になります。

https://qiita.com/advent-calendar/2023/ppl

先に公開していた下記のメモの完成版です（投稿が大幅に遅れ大変申し訳ございません）。

https://zenn.dev/tatamiya/scraps/9d14f3d1ee5a37

この記事は？

本記事では、ベイズ推定において用いられる Spike-and-Slab 事前分布について解説します。

Spike-and-Slab 事前分布は線形回帰の際に変数選択を行うために導入するもので、時系列データによる施策効果検証でよく用いられる CausalImpact にも組み込まれています。^[1]

本記事では主に CausalImpact とその関連論文 [Brodersen et al. (2015)] [Scott and Varian (2014)] に沿って Spike-and-Slab 事前分布の理論的な側面を紹介したのち、同事前分布を用いた際の挙動を視覚的・直感的に理解しやすいように PyMC を用いてモデルを作成しました。
そのうえで、Spike-and-Slab 事前分布を用いないベイズ線形回帰モデルと推定結果の比較を行いました。

サマリー

• Spike-and-Slab 事前分布は、各説明変数をモデルに含めるか含めないかをベルヌーイ分布に従う潜在変数によって決定することで、確率的に変数選択を行う。
• 同事前分布を用いることで寄与のない説明変数がモデルから除外されやすくなるため、通常のベイズ線形回帰モデルと比べてノイズにオーバーフィットしにくく滑らかな予測結果が得られる。

Spike-and-Slab 事前分布とは？

Spike-and-Slab 事前分布は、線形回帰モデル作成の際に各説明変数をモデルに含めるか含めないかをベルヌーイ分布に従って決定するものである。
この際、予測への寄与が小さい変数ほどモデルに採り入れられる（事後）確率が小さくなる。

したがってこれを用いてモデルを作成すると確率的に変数選択が行われる。
もう少し詳しく説明すると、寄与の大きい説明変数を優先的に採り入れつつ異なる説明変数を組み合わせた複数のモデルが作成され、それらのいわば平均値が予測・推定結果として得られる。

設定

以下のような $pp$ 次元の説明変数 $xi=(x1i,x2i,...,xpi)⊤\boldsymbol{x}_i = (x_{1i}, x_{2i}, ..., x_{pi})^\top$ により目的変数 $yiy_i$ を予測する線形回帰モデルを考える：

$$yi=β1x1i+β2x2i+...+βpxpi+εi=xi⊤β+εi. \begin{aligned} y_i &= \beta_1 x_{1i} + \beta_2 x_{2i} + ... + \beta_p x_{pi} + \varepsilon_i\\ &= \boldsymbol{x}_i^\top \boldsymbol{\beta} + \varepsilon_i. \end{aligned}$$

ここで、$εi\varepsilon_i$ は期待値 0 分散 $σ2\sigma^2$ の独立同時な正規分布 $N(0,σ2)N(0, \sigma^2)$ に従う誤差項、$β=(β1,β2,...,βp)⊤\boldsymbol{\beta} = (\beta_1, \beta_2, ..., \beta_p)^\top$ は回帰係数である。

サンプルサイズを $nn$ として、上記のモデルは目的変数ベクトル $Y=(y1,y2,...,yn)⊤\boldsymbol{Y} = (y_1, y_2, ..., y_n)^\top$、計画行列 $X=(x1,x2,...,xn)⊤\boldsymbol{X} = (\boldsymbol{x}_1, \boldsymbol{x}_2, ..., \boldsymbol{x}_n)^\top$、誤差項ベクトル $ε=(ε1,ε2,...,εn)⊤\boldsymbol{\varepsilon} = (\varepsilon_1, \varepsilon_2, ..., \varepsilon_n)^\top$ を用いて以下のように表すことができる：

$$Y=Xβ+ε. \boldsymbol{Y} = \boldsymbol{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon}.$$

これは、$β,σ2\boldsymbol{\beta}, \sigma^2$ による条件付き分布として書き表すこともできる：

$$Y∣β,1/σ2∼N(Xβ,σ2In). \boldsymbol{Y}\vert \boldsymbol{\beta}, 1/\sigma^2 \sim N(\boldsymbol{X}\boldsymbol{\beta}, \sigma^2 \boldsymbol{I}_n).$$

$In\boldsymbol{I}_n$ は $n×nn \times n$ の単位行列である。

以降、上記のようなモデル設定をもとにパラメータ $β,σ2\boldsymbol{\beta}, \sigma^2$ をベイズ推定することを考える。

ベイズ線形回帰

まずは Spike-and-Slab 事前分布を用いずに推定を行うことを考える。

この際、パラメータ $β,σ2\boldsymbol{\beta}, \sigma^2$ の事前分布は以下のように設定する：

$$β∣1/σ2∼N(b,σ2Σ)1/σ2∼Ga(ν2,s2) \begin{aligned} \boldsymbol{\beta} \vert 1/\sigma^2 &\sim N(\boldsymbol{b}, \sigma^2 \boldsymbol{\Sigma})\\ 1/\sigma^2 &\sim {\rm Ga}\left(\frac{\nu}{2}, \frac{s}{2}\right) \end{aligned}$$

$N(μ,Σ)N(\boldsymbol{\mu}, \boldsymbol{\Sigma})$ は期待値が $μ\boldsymbol{\mu}$ 分散共分散行列が $Σ\boldsymbol{\Sigma}$ の多変量正規分布、$Ga(a,b){\rm Ga}(a, b)$ は期待値が $a/ba/b$ で分散が $a/b2a/b^2$ のガンマ分布である。

これらをグラフィカルモデルで表すと以下のようになる：

このとき、同時分布が

$$P(Y,β,1/σ2)=P(Y∣β,1/σ2)P(β∣,1/σ2)P(1/σ2) P(\boldsymbol{Y}, \boldsymbol{\beta}, 1/\sigma^2) = P(\boldsymbol{Y} \vert \boldsymbol{\beta}, 1/\sigma^2) P(\boldsymbol{\beta} \vert , 1/\sigma^2) P(1/\sigma^2)$$

と表せることから、事後分布は以下のようになる：

$$1/σ2∣Y∼Ga(n+ν2,SY2),β∣Y,1/σ2∼N(β~,σ2Ω). \begin{aligned} 1/\sigma^2 \vert \boldsymbol{Y} &\sim {\rm Ga}\left( \frac{n+\nu}{2}, \frac{S_{Y}}{2}\right),\\ \boldsymbol{\beta} \vert \boldsymbol{Y}, 1/\sigma^2 &\sim N(\tilde{\boldsymbol{\beta}}, \sigma^2 \boldsymbol{\Omega}). \end{aligned}$$

ただし、

$$Ω=(X⊤X+Σ−1)−1,β~=Ω(Σ−1b+X⊤Y),SY=12(s+Y⊤Y+b⊤Σ−1b−β~⊤Ω−1β~). \boldsymbol{\Omega} = (\boldsymbol{X}^\top \boldsymbol{X} + \boldsymbol{\Sigma}^{-1})^{-1}, \quad \tilde{\boldsymbol{\beta}} = \boldsymbol{\Omega} (\boldsymbol{\Sigma}^{-1} \boldsymbol{b}+ \boldsymbol{X}^\top \boldsymbol{Y}),\\ S_Y = \frac{1}{2}(s + \boldsymbol{Y}^\top \boldsymbol{Y} + \boldsymbol{b}^\top \boldsymbol{\Sigma}^{-1}\boldsymbol{b} - \tilde{\boldsymbol{\beta}}^\top \boldsymbol{\Omega}^{-1}\tilde{\boldsymbol{\beta}}).$$

Spike-and-Slab 事前分布の導入

次に、$pp$ 個の説明変数に対応して以下のような 0 or 1 の値を取る $pp$ 個の互いに独立なベルヌーイ分布に従う確率変数 $ρ=(ρ1,ρ2,...,ρp)⊤\boldsymbol{\rho} = (\rho_1, \rho_2, ..., \rho_p)^\top$ を定義する：

$$P(ρ)=∏k=1pπkρk(1−πk)1−ρk(0≤πk≤1). P(\boldsymbol{\rho}) = \prod_{k=1}^p \pi_k^{\rho_k} (1 - \pi_k)^{1 - \rho_k} \qquad (0\le\pi_k\le 1).$$

これを用いて、$ρk=1\rho_k=1$ なら $kk$ 番目の説明変数 $xkx_k$ をモデルに含め、$ρk=0\rho_k=0$ なら含めないようにする。
たとえば、$p=6p=6$ で $ρ=(1,0,0,0,0,1)\boldsymbol{\rho} = (1,0,0,0,0,1)$ の場合、線形回帰モデルは以下のような形になる：

$$yi=xρ,i⊤βρ+εi=β1x1i+β6x6i+εi \begin{aligned} y_i &= \boldsymbol{x}_{\rho, i}^\top \boldsymbol{\beta}_\rho + \varepsilon_i\\ &= \beta_1 x_{1i} + \beta_6 x_{6i} + \varepsilon_i \end{aligned}$$

以降、$X,β\boldsymbol{X}, \boldsymbol{\beta}$ から $ρk=1\rho_k=1$ に対応する成分だけを残したものを $Xρ,βρ\boldsymbol{X}_\rho, \boldsymbol{\beta}_\rho$ と表す。

このような設定下では、$β,1/σ2,ρ\boldsymbol{\beta}, 1/\sigma^2, \boldsymbol{\rho}$ で条件づけた目的変数ベクトル $Y\boldsymbol{Y}$ の分布は以下のように表すことができる：

$$Y∣β,1/σ2,ρ∼N(Xρβρ,σ2In) \boldsymbol{Y} \vert \boldsymbol{\beta}, 1/\sigma^2, \boldsymbol{\rho} \sim N(\boldsymbol{X}_\rho \boldsymbol{\beta}_\rho, \sigma^2 \boldsymbol{I}_n)$$

ここで回帰係数 $βρ\boldsymbol{\beta}_\rho$ と誤差項分散の逆数 $1/σ21/\sigma^2$ について、先ほどと同様に以下のような事前分布を置く：

$$βρ∣1/σ2,ρ∼N(bρ,σ2(Σρ−1)−1)1/σ2∼Ga(ν2,s2) \begin{aligned} \boldsymbol{\beta}_\rho \vert 1/\sigma^2, \boldsymbol{\rho} &\sim N(\boldsymbol{b}_\rho, \sigma^2 (\boldsymbol{\Sigma}_\rho^{-1})^{-1})\\ 1/\sigma^2 &\sim {\rm Ga}\left(\frac{\nu}{2}, \frac{s}{2}\right) \end{aligned}$$

ただし、$Σρ−1\boldsymbol{\Sigma}_\rho^{-1}$ は、$Σ−1\boldsymbol{\Sigma}^{-1}$ のうち $ρk=1\rho_k=1$ となる行・列のみを取り出して構成した行列である。^[2]

グラフィカルモデルで表すと以下のようになる：

このとき同時分布は

$$P(Y,β,1/σ2,ρ)=P(Y∣β,1/σ2,ρ)P(β∣1/σ2,ρ)P(1/σ2)P(ρ) P(\boldsymbol{Y}, \boldsymbol{\beta}, 1/\sigma^2, \boldsymbol{\rho}) = P(\boldsymbol{Y} \vert \boldsymbol{\beta}, 1/\sigma^2, \boldsymbol{\rho}) P(\boldsymbol{\beta} \vert 1/\sigma^2, \boldsymbol{\rho}) P(1/\sigma^2) P(\boldsymbol{\rho})$$

であるので、事後分布は以下のように求まる：

$$1/σ2∣Y,ρ∼Ga(n+ν2,SY,ρ2),βρ∣Y,1/σ2,ρ∼N(β~ρ,σ2(Ωρ−1)−1),P(ρ∣Y)=1P(Y)1πn2∣Σρ−1∣12∣Ωρ−1∣12Γ(n+ν2)Γ(ν2)sν2SY,ρn+ν2P(ρ).(1) \begin{aligned} 1/\sigma^2 \vert \boldsymbol{Y}, \boldsymbol{\rho} &\sim {\rm Ga}\left( \frac{n+\nu}{2}, \frac{S_{Y,\rho}}{2}\right),\\ \boldsymbol{\beta}_\rho \vert \boldsymbol{Y},1/\sigma^2, \boldsymbol{\rho} &\sim N(\tilde{\boldsymbol{\beta}}_\rho, \sigma^2 (\boldsymbol{\Omega}_\rho^{-1})^{-1}),\\ P(\boldsymbol{\rho} \vert \boldsymbol{Y}) &= \frac{1}{P(\boldsymbol{Y})} \frac{1}{\pi^{\frac{n}{2}}} \frac{|\boldsymbol{\Sigma}_\rho^{-1}|^{\frac{1}{2}}}{|\boldsymbol{\Omega}_\rho^{-1}|^{\frac{1}{2}}} \frac{\Gamma\left( \frac{n+\nu}{2}\right)}{\Gamma\left( \frac{\nu}{2}\right)} \frac{s^{\frac{\nu}{2}}}{S_{Y,\rho}^{\frac{n+\nu}{2}}} P(\boldsymbol{\rho}). \tag{1} \end{aligned}$$

ただし、

$$(Ωρ−1)−1=(Xρ⊤Xρ+Σρ−1)−1,β~ρ=(Ωρ−1)−1(Σρ−1bρ+Xρ⊤Y),SY,ρ=s+Y⊤Y+bρ⊤Σρ−1bρ−β~ρ⊤Ωρ−1β~ρ (\boldsymbol{\Omega}_\rho^{-1})^{-1} = (\boldsymbol{X}_\rho^\top \boldsymbol{X}_\rho + \boldsymbol{\Sigma}^{-1}_\rho)^{-1}, \quad \tilde{\boldsymbol{\beta}}_\rho = (\boldsymbol{\Omega}_\rho^{-1})^{-1} (\boldsymbol{\Sigma}_\rho^{-1} \boldsymbol{b}_\rho + \boldsymbol{X}_\rho^\top \boldsymbol{Y}),\\ S_{Y,\rho} = s + \boldsymbol{Y}^\top \boldsymbol{Y} + \boldsymbol{b}_\rho^\top \boldsymbol{\Sigma}_\rho^{-1}\boldsymbol{b}_\rho - \tilde{\boldsymbol{\beta}}_\rho^\top \boldsymbol{\Omega}_\rho^{-1}\tilde{\boldsymbol{\beta}}_\rho$$

とした。

PyMC を用いた数値実験

以下では、Spike-and-Slab 事前分布を用いた推定結果の振る舞いを理解しやすいように、シンプルな人工データを用いて数値実験を行った。

Spike-and-Slab 線形回帰の実装は式(1)をもとに Gibbs サンプリングを行う方法が知られているが [George and McCulloch (1993, 1997)] [ホフ (2022)]、本記事では実装のしやすさと Spike-and-Slab なしの推定結果との比較のしやすさの点から PyMC を用いた類似実装を行った。^[3]

実験データ

目的変数

本記事の実験におけるサンプルデータでは、目的変数 ${y}i=1n\{y\}_{i=1}^n$ として以下のように2種類の異なる周期を持つ sin, cos 波を重ね合わせた波形信号にガウスノイズを加えて生成した人工データを用いた。

$$yi=2sin⁡(2πxi8)−cos⁡(2πxi2)+εi,εi∼N(0,1) y_i = 2 \sin \left(2\pi \frac{x_i}{8}\right) - \cos \left(2\pi \frac{x_i}{2}\right) + \varepsilon_i, \quad \varepsilon_i \sim N(0, 1)$$

なお、サンプルサイズは $n=17n=17$ とし、${xi}i=1n\{x_i\}_{i=1}^n$ を $[0,8][0,8]$ の区間から0.5刻みで等間隔に採った。

import numpy as np

np.random.seed(0)

# x は指定区間に均一に生成する
n_samples = 16

x_min = 0.0
x_max = 8.0

x = x_min + (x_max - x_min) * np.arange(0, n_samples + 1, 1) / n_samples

# 誤差項分散
sigma = 1.0
sigma2 = sigma **2

# 目的変数 y の生成
def target_func(x):
  return 2 * np.sin(2 * np.pi / 8 * x) - 1.0 * np.cos(4 * 2 * np.pi / 8 * x)

y =  target_func(x) + sigma * np.random.randn(x.shape[0])

説明変数とモデル設定

説明変数としては、以下のように3つの異なる周期を持つ sin, cos 関数を用いた。

$$xi=(x1i,x2i,x3i,x4i,x5i,x6i)⊤=(sin⁡(2πxi8),cos⁡(2πxi8),sin⁡(2πxi4),cos⁡(2πxi4),sin⁡(2πxi2),cos⁡(2πxi2))⊤ \begin{aligned} \boldsymbol{x}_i &= (x_{1i}, x_{2i}, x_{3i}, x_{4i}, x_{5i}, x_{6i})^\top\\ &= \left(\sin\left(2\pi \frac{x_i}{8}\right), \cos\left(2\pi \frac{x_i}{8}\right), \sin\left(2\pi \frac{x_i}{4}\right), \cos\left(2\pi \frac{x_i}{4}\right), \sin\left(2\pi \frac{x_i}{2}\right), \cos\left(2\pi \frac{x_i}{2}\right) \right)^\top \end{aligned}$$

これらを用いて $p=6p=6$ の定数項なしの線形回帰モデル

$$yi=xi⊤β+εi=β1x1i+β2x2i+β3x3i+β4x4i+β5x5i+β6x6i+εi \begin{aligned} y_i &= \boldsymbol{x}_i^\top \boldsymbol{\beta} + \varepsilon_i\\ &= \beta_1 x_{1i} + \beta_2 x_{2i} + \beta_3 x_{3i} + \beta_4 x_{4i} + \beta_5 x_{5i} + \beta_6 x_{6i} + \varepsilon_i \end{aligned}$$

を作成した。

目的変数 $yiy_i$ の生成には $sin⁡(2πxi/8),cos⁡(2πxi/2)\sin\left(2\pi x_i/8\right), \cos\left(2\pi x_i / 2\right)$ 成分しか用いていないので、$β1=2,β6=−1\beta_1=2, \beta_6=-1$ のほかは $β2=β3=β4=β5=0\beta_2=\beta_3=\beta_4=\beta_5=0$ となる想定である。

PyMC によるモデル実装

PyMC で Spike-and-Slab 事前分布を用いた線形回帰モデルを実装するにあたり、[George and McCulloch (1993, 1997)] を参考にしつつ前述したモデルに多少改変を加えた。

なぜ CausalImpact で用いられているものと同様の Spike-and-Slab 回帰を実装しなかったかというと、PyMC では特定のパラメータ（ここでは回帰係数 $β\boldsymbol{\beta}$）の次元を別のパラメータ $ρ\boldsymbol{\rho}$ の値に応じて動的に変えることが難しかったからである。^[4]

そこで以下に示すように、$ρk=0\rho_k=0$ に対応する説明変数をモデルから除外するのではなく、回帰係数 $βk\beta_k$ に 0 付近に局在した（$ρk=1\rho_k=1$ の場合と比べると分散の小さい）事前分布を与えた。

$$Y∣β,1/σ2∼N(Xβ,σ2In)β∣1/σ2,ρ∼N(Iρb,σ2(Σ~ρ−1)−1)1/σ2∼Ga(ν2,s2) \begin{aligned} \boldsymbol{Y} \vert \boldsymbol{\beta}, 1/\sigma^2 &\sim N(\boldsymbol{X} \boldsymbol{\beta}, \sigma^2 \boldsymbol{I}_n)\\ \boldsymbol{\beta} \vert 1/\sigma^2, \boldsymbol{\rho} &\sim N(\boldsymbol{I}_\rho \boldsymbol{b}, \sigma^2 (\tilde{\boldsymbol{\Sigma}}_\rho^{-1})^{-1})\\ 1/\sigma^2 &\sim {\rm Ga}\left(\frac{\nu}{2}, \frac{s}{2}\right) \end{aligned}$$

$$Σ~ρ−1=IρΣ−1Iρ+n2(Ip−Iρ)diag(Σ−1) \tilde{\boldsymbol{\Sigma}}_\rho^{-1} = \boldsymbol{I}_\rho \boldsymbol{\Sigma}^{-1} \boldsymbol{I}_\rho + n^2 (\boldsymbol{I}_p - \boldsymbol{I}_\rho) {\rm diag}(\boldsymbol{\Sigma}^{-1})$$

$Iρ\boldsymbol{I}_\rho$ は対角行列で $ρk=1\rho_k=1$ に対応する成分が1でそれ以外は0となっている。
$diag(Σ−1){\rm diag}(\boldsymbol{\Sigma}^{-1})$ は $Σ−1\boldsymbol{\Sigma}^{-1}$ と対角成分が同一の対角行列である。
また $Σ~ρ−1\tilde{\boldsymbol{\Sigma}}_\rho^{-1}$ は、$\boldsymbol{\Sigma}^{-1}$ のうち $\rho_k=0$ に対応する行・列について対角成分だけ $n^2$ 倍してそれ以外の成分は0になるように作成している。

事前分布についてのハイパーパラメータは、CausalImpact と合わせるために論文 [Brodersen et al. (2015)] [Scott and Varian (2014)] の記述に基づいて以下のように設定した：

$$\boldsymbol{b} = \boldsymbol{0}, \quad \boldsymbol{\Sigma}^{-1} = \frac{1}{n} \frac{1}{2} \left\{ \boldsymbol{X}^\top \boldsymbol{X} + {\rm diag}(\boldsymbol{X}^\top \boldsymbol{X})\right\},\\ \nu=0.01, \quad s = \nu (1 - R^2) s_y,\\ R^2 = 0.5, \quad s_y = \sum_{i=1}^n (y_i - \bar{y})^2 / (n-1).$$

$\bar{y}$ は $\{y_i\}_{i=1}^n$ のサンプル平均である。

これを実装したものが以下である：

import pymc as pm
import pytensor.tensor as pt
import numpy as np

def get_model(y, X):
    n = y.shape[0]
    p = X.shape[1]

    model = pm.Model()

    Sigma = 0.5 * np.matmul(X.T, X) / n
    Sigma += np.diag(np.diag(Sigma))

    s_y2 = np.var(y, ddof=1)

    nu = 0.01
    R2 = 0.5
    s = nu * (1 - R2) * s_y2

    with model:
        rho = pm.Bernoulli('rho', .5, shape=p)

        tau = pm.Gamma('tau', nu/2, s/2)

        restricted_Sigma = pt.diag(rho) * Sigma * pt.diag(rho) + n **2 * (pt.identity_like(Sigma) - pt.diag(rho)) * np.diag(np.diag(Sigma))
        beta = pm.MvNormal('beta', 0, tau=tau * restricted_Sigma, shape=p)

        mean = pm.math.dot(X, beta)
        y_obs = pm.Normal('y_obs', mean, tau=tau, observed=y)
    return model

# 説明変数をまとめた計画行列（design matrix）の生成
def generate_design_matrix(x):
  return np.array([
      np.sin(2 * np.pi / 8 * x), np.cos(2 * np.pi / 8 * x),
      np.sin(2 * 2 * np.pi / 8 * x), np.cos(2 * 2 * np.pi / 8 * x),
      np.sin(4 * 2 * np.pi / 8 * x), np.cos(4 * 2 * np.pi / 8 * x)
      ]).T
x_des = generate_design_matrix(x)

# MCMC によるモデル推定
model = get_model(y, x_des)
with model:
    # draw 1000 posterior samples
    idata = pm.sample()

モデル推定結果

以降では、Spike-and-Slab 事前分布を用いたモデルの結果を、同事前分布を用いずに作成したベイズ線形回帰モデルと比較しながら紹介する。

事後予測値の比較

以下のグラフは、目的変数 $y$ の事後予測値を Spike-and-Slab 事前分布なし・ありで比較したものである。
青の実線はモデルによる事後予測平均値、その周りの薄青の領域は95%区間を表している。
黒点は観測データ、灰色の点線は誤差項を含まない $y$ の真の値である。

これを見ると、Spike-and-Slab 事前分布を用いない場合（左）では事後予測平均値（青線）が $x=1\sim2$ 付近などで観測誤差に引きずられて予測値が大きく上振れており全体的に変動が大きいものの、同事前分布を用いた場合（右）では予測値が真の値（灰色点線）の間を縫うように滑らかに変動している。

パラメータ事後分布の比較

次に、パラメータの事後分布を比較する。

以下は、回帰係数 $\beta_1, \beta_2, ..., \beta_6$ と誤差項分散の逆数 $\tau = 1/\sigma^2$ の事後分布 $P(\beta_k \vert \boldsymbol{Y}), \, P(1/\sigma^2 \vert \boldsymbol{Y})$ を可視化したものである。青線で表した区間は95%の High Density Interval (HDI) ^[5]であり、線が太くなっている部分は50% HDI を表している。また、白抜きの◯印は分布の中央値（50%点）、赤い星印はサンプルデータ生成時に指定した真の値である。
実装の都合上、回帰係数のインデックスが1つずつずれていることにご注意いただきたい（例：beta[0] は $\beta_1$ に対応）。

これを見ると、$\beta_1$ (beta[0]) については Spike-and-Slab 事前分布を用いない場合（左）と用いた場合（右）で後者の方が HDI がやや広いことを除けば事後分布に大きな差は見られなかった。
一方で、真の回帰係数が0の $\beta_2, \beta_3, \beta_4, \beta_5$ (beta[1], [2], [3], [4]) については Spike-and-Slab 事前分布を用いた方（右）が事後分布の中心（中央値）が0に近く、また 50% HDI の幅も短くなっていることから分布自体も0付近に局在していることが見て取れる。

なお、$\beta_6$ (beta[5]) と $\tau=1/\sigma^2$ (tau) については、Spike-and-Slab 事前分布を使わないモデルの方が事後分布の中央値が真の値に近かった。
このうち後者については同事前分布を用いた方が $\tau=1/\sigma^2$ の推定値が全体的に小さくなっており、これは誤差項分散が本来の値より大きく推定される傾向にあることを示している。

Spike-and-Slab を用いた場合の事後分布の詳細

上記では Spike-and-Slab 事前分布を使った場合では真の値が0の回帰係数 $\beta_2, \beta_3, \beta_4, \beta_5$ について事後分布が0に局在していた。

このことをより詳細に調べるために、パラメータ $\boldsymbol{\rho}$ による変数選択がどのように行われていたかと、$\rho_k=0,1$ 別に見た回帰係数 $\beta_k$ の事後分布を見てみる。

変数選択の事後確率

下記の棒グラフは、各変数が選択された事後確率に相当する $P(\rho_k =1 \vert \boldsymbol{Y})$ を可視化したものである。

これを見ると、$x_1 = \sin\left(2\pi x /8\right)$ はほぼ100%モデルに取り込まれているのに対して、真の回帰係数を0に設定した $x_2, x_3, x_4, x_5$ については25~30%程度であることがわかる。
なお、真の回帰係数を $\beta_6=-1$ に設定した $x_6 = \cos\left(2\pi x / 2\right)$ については選択事後確率が49%であり、 $x_2,x_3, x_4, x_5$ と比べてやや高い程度に留まった。^[6]

回帰係数事後分布の詳細

次に、$\rho_k=0, 1$ のそれぞれの場合に分けて回帰係数 $\beta_k$ の事後分布 $P(\beta_k \vert \boldsymbol{Y}, \rho_k)$ を見てみたものが下記のグラフである。
オレンジ色が $\rho_k=1$ （説明変数 $x_k$ が選択された場合に相当）における $\beta_k$ の事後分布、青が $\rho_k=0$ における事後分布である。縦の破線はサンプルデータ作成時に与えた $\beta_k$ の真の値を表している。
なお、$\rho_k=0,1$ のどちらの事後確率が高かったかが視覚的にわかるように、各分布を選択事後確率で重み付けし $P(\beta_k \vert \boldsymbol{Y}, \rho_k)P(\rho_k \vert \boldsymbol{Y})$ に相当する内容を描画している。したがって、分布単体で見ると規格化されていない。$\rho_k=0,1$ の分布を足し合わせたものが、先ほど forest plot で見ていたものと対応する。

これを見ると、$\beta_2, \beta_3, \beta_4, \beta_5$ では $\rho_k=0$ における事後分布が0を中心とした分布になっており、さらにその形状は $\rho_k=1$ の場合と比べて細く鋭いことがわかる。

まとめ・考察

Spike-and-Slab 事前分布は線形回帰の際に不要な説明変数をなるべくモデルから除外する目的で導入されたものであった。

実際に本記事で適用した例では予測への寄与のない説明変数 $(x_2, x_3, x_4, x_5)$ の選択事後確率 $P(\rho_k=1 \vert \boldsymbol{Y})$ が他と比べて低くなることが確認できた。
さらにそれらの回帰係数 $\beta_2, \beta_3, \beta_4, \beta_5$ の事後分布 $P(\beta_k \vert \boldsymbol{Y})$ を見ると、Spike-and-Slab 事前分布を用いなかった場合と比べて事後分布の中央値が 0 に近くかつ50% HDI の幅も狭かった。

目的変数 $y$ の事後予測結果で比べると Spike-and-Slab 事前分布を用いた場合は用いなかった場合と比べて観測ノイズへのオーバーフィットが少なく真値の大域的な特徴を反映したものとなっていた。
このことは、上記に述べた Spike-and-Slab 事前分布の性質により予測への寄与のない説明変数をなるべく含めずに（実際には0に近い回帰係数を割り当てて）モデルを作成したことで、モデルの複雑さが抑えられたためと考えられる。

不要な説明変数がない場合の挙動

一方で、寄与のない説明変数が含まれない場合に Spike-and-Slab 事前分布を導入した場合はあまり効果的でないことも下記の実験から伺うことができた。

以下は、$(\beta_1, \beta_2, \beta_3, \beta_4, \beta_5, \beta_6) = (2, 1.0, 0.5, -0.5, -1.5, -1.0)$ といずれの回帰係数も 0 でない場合に同様の実験を行った結果である。

この場合、目的変数 $y$ の事後予測値の平均値については Spike-and-Slab 事前分布を使用した場合としなかった場合とで先ほどと比べて大きな差は見られなかった。
その一方で、95%区間の幅の広さは Spike-and-Slab 事前分布を入れた方が広くなっており、予測の不確実性が大きいことがわかる。

回帰係数の事後分布についても同様に、Spike-and-Slab 事前分布を用いたモデルの方が信用区間の幅が広くなっていることがわかる。

参考までに、各説明変数の選択事後確率は以下のようになっていた。

結論

以上から、Spike-and-Slab 事前分布の導入は、モデルに取り入れたい説明変数に予測に寄与しないものが含まれている場合において、確率的な変数選択により学習データに対するオーバーフィットを抑える効果を発揮するものと考えられる。

一方、あらかじめ全説明変数が予測に一定の寄与をすることがわかっている場合はわざわざ同事前分布を用いる必要はないだろう。

実行コード

こちら

参考文献

• Brodersen, Kay H., Fabian Gallusser, Jim Koehler, Nicolas Remy, and Steven L. Scott. "Inferring causal impact using Bayesian structural time-series models." (2015): 247-274.

  • CausalImpact の論文

• Scott, Steven L., and Hal R. Varian. "Predicting the present with Bayesian structural time series." International Journal of Mathematical Modelling and Numerical Optimisation 5, no. 1-2 (2014): 4-23.

• George, Edward I., and Robert E. McCulloch. "Variable selection via Gibbs sampling." Journal of the American Statistical association (1993): 881-889.
  

• George, Edward I., and Robert E. McCulloch. "Approaches for Bayesian variable selection." Statistica sinica (1997): 339-373.

• ピーター・D・ホフ（入江 薫・菅澤 翔之助・橋本 真太郎 訳）、「標準 ベイズ統計学」（朝倉書店、2022）

  • 第9章（線形回帰）に、本稿で紹介した内容とほぼ同一のモデルが、 R によるギブスサンプラーを用いた実装と共に詳述されている。

https://amzn.asia/d/5dqGr9G

• Batı Şengül, Spike and slab: Bayesian linear regression with variable selection (2018)

  • PyMC による実装時に参考にした。

• PyMC の公式ドキュメント

脚注

1. CausalImpact ではまず説明変数を含めずにトレンド成分・季節成分を推定し、観測値からそれらの成分を差し引いた残成分に対して Spike-and-Slab 事前分布を用いた線形回帰を行う。[Brodersen et al. (2015)] [Scott and Varian (2014)]

2. 一般に $\boldsymbol{\Sigma}$ から $\rho_k=1$ となる行・列のみを取り出して構成した行列についての逆行列を作成しても、本文中で述べた $\boldsymbol{\Sigma}_\rho^{-1}$ と一致しない。$\boldsymbol{\beta}_\rho$ についての事前分布の中で $(\boldsymbol{\Sigma}_\rho^{-1})^{-1}$ を $\boldsymbol{\Sigma}_\rho$ と書き表さなかったのはそのような理由からである。

3. なお、CausalImpact の作成者の一人 Scott 氏による Spike-and-Slab 回帰の R パッケージ BoomSpikeSlab が公開されているものの、筆者の手元で動かそうとしたところインストールに失敗した。その際のメモはこちら。

4. こちら にて Matsuura さんから指摘いただいたように、離散変数の数え上げ消去を用いることで CausalImpact と同様の Spike-and-Slab 回帰の実装を行えそうである。こちらについても別途試す予定である。

5. High Density Interval (HDI) については、次の記事が詳しい https://www.pagumi-bayesian.com/2021/09/11/hdi/

6. 前節で触れたように $\beta_6$ の事後分布中央値が Spike-and-Slab 事前分布を用いたモデルの方が真の値から離れてしまっていたのは、このように変数選択の事後確率 $P(\rho_6 = 1 \vert \boldsymbol{Y})$ が 49% 程度でありおおよそ半々の確率で係数が 0 に近い値と判断されたためと考えられる。そもそも $\beta_6$ に対応する説明変数 $x_6 = \cos(2\pi x / 2)$ は周期が 2 と短く、$x$ の間隔が 0.5 刻みの学習データでは回帰係数推定が困難だったと考えられる。