はじめに
今回はARモデルのベイズ的な取り扱いを勉強していきます。扱う範囲として、尤度関数を求め無情報事前を使用し、パラメータの事後確率分布を計算そして周辺化。最後に計算した周辺分布から実際にパラメータをサンプリングしてみるとこまで行います。
導入として線形モデルを扱います。ARモデルの場合もほとんど同じように考えることができます。線形モデルでの結果を利用してARモデルも考えていくので、ぜひ線形モデルからご覧になってください。
ベイズ統計は知りたい情報に関する確率分布を求めることが目的です。
この知りたい情報に関する事後確率分布が予測分布と呼ばれるものです。
その他のデータに条件付けられた確率分布は単に事後確率分布と呼ばれます。
今回の記事はパラメータの事後確率分布を計算することを目標にしています。但し、これは本来のベイズ統計の目的とは異なっていることに注意してください。（簡単に言うと事後確率分布を求めることがゴールではないということ）
線形モデルから始めよう！
線形モデルを以下のように定義します。
y = Xb + \epsilon \quad \epsilon ～ \mathcal{N}(0, vI_n)
X = \begin{bmatrix}
f_1(x_1) & \cdots & f_r(x_1) \\
\vdots & & \vdots \\
f_1(x_n) & \cdots & f_r(x_n) \\
\end{bmatrix}, \quad
b = \begin{bmatrix}
b_1 \\
\vdots \\
b_r \\
\end{bmatrix}, \quad
y = \begin{bmatrix}
y_1 \\
\vdots \\
y_n \\
\end{bmatrix}
ここで、$y = \begin{pmatrix} y_1, y_2, ... , y_n \end{pmatrix}'$ が得られる確率（尤度関数）は次で表せます。但し、 $x_i$ は $X$ の $i$ 行目です。
\begin{align}
P(y|X, b, v) &= \Pi_{i=1}^{n} \mathcal{N} ( y_i | x_i b , v ) \\
&= ( 2 \pi v )^{- \frac{n}{2}} \exp( -\frac{(y - Xb)'(y - Xb)}{2v} )
\end{align}
ここで、$P(y|X, b, v)$ を最大化するようなパラメータ $\hat{b}$ は次のようになります。
\hat{b} = (X'X)^{-1}X'y
$\hat{b}$ の求め方は、$P(y|X, b, v)$ に対数を取り $b$ で微分しても出てきます。他の考え方として、以前記事にした直交射影からでも求めることができます(こちらのほうが簡単)。
さて後々のために $Q(y,b) = (y - Xb)'(y - Xb)$ を使って指数部を書き直します。
P(y|X, b, v) = ( 2 \pi v )^{- \frac{n}{2}} \exp( -\frac{Q(y,b)}{2v} ) \tag{1} \\
ここで、$Q(y,b)$ は次のように書き換えることができます。
\begin{align}
Q(y,b) &= (y - Xb)'(y - Xb) \\
&= y'y - 2y'Xb + b'X'Xb \\
&= y'y - 2y'Xb + b'X'Xb - 2y'X\hat{b} + 2y'X\hat{b} - \hat{b}'X'X\hat{b} + \hat{b}'X'X\hat{b} \\
&= b'X'Xb - \hat{b}'X'X\hat{b} - 2y'X(X'X)^{-1}(X'X)b + 2y'X(X'X)^{-1}(X'X)\hat{b} + (y - X\hat{b})'(y - X\hat{b}) \\
&= b'X'Xb - \hat{b}'X'X\hat{b} - 2\hat{b}'X'Xb + 2\hat{b}'X'X\hat{b} + (y - X\hat{b})'(y - X\hat{b}) \\
&= b'X'Xb + \hat{b}'X'X\hat{b} - 2\hat{b}'X'Xb + (y - X\hat{b})'(y - X\hat{b}) \\
&= b'X'Xb + \hat{b}'X'X\hat{b} - b'X'X\hat{b} - \hat{b}'X'Xb + (y - X\hat{b})'(y - X\hat{b}) \\
&= b'X'X(b - \hat{b}) - \hat{b}'X'X(b - \hat{b}) + (y - X\hat{b})'(y - X\hat{b}) \\
&= (b - \hat{b})'X'X(b - \hat{b}) + (y - X\hat{b})'(y - X\hat{b}) \\
\end{align}
$(y - X\hat{b})'(y - X\hat{b})$ は残差平方和（Residual Sum of Squares）です。これを $n - p$ で割ると v の不偏推定量が求められます。ここで $n$ はデータ数、 $p$ はパラメータの数です。
\begin{align}
s^2 &= \frac{R}{n-p} \\
&= \frac{(y - X\hat{b})'(y - X\hat{b})}{n - p} \tag{2}
\end{align}
式変形は最終形が、「それ以外」＋「残差平方和」となるように行っています。少し天下り的に変形を行っています。腑に落ちないのであれば、$Q(y,b)- 残差平方和$ から計算してみるといいかもしれません。同様の結果が得られます。
残差平方和の分解
残差平方和も少し式変形してみましょう。表記の都合上、以降、残差平方和を $R$ と置くことがあります。
\begin{align}
R &= (y - X\hat{b})'(y - X\hat{b}) \\
&= y'y - 2y'X\hat{b} + \hat{b}'X'X\hat{b} \\
&= y'y - 2y'X(X'X)^{-1}X'y + y'X(X'X)^{-1}X'X(X'X)^{-1}X'y \\
&= y'y - y'X(X'X)^{-1}X'y
\end{align}
多変数なのでわかりづらいかもしれないですが、$y'y$ が全変動（Total Sum of Squares）、$y'X(X'X)^{-1}X'y$ は回帰変動（Explained Sum of Squares）です。
$$
y'y - R = y'X(X'X)^{-1}X'y
$$
回帰変動が大きければ、残差平方和が小さい。つまりよりデータにフィットしていると解釈できます。ここらへんは1変数の場合と同じですね。
別解
$X(X'X)^{-1}X'y$ が $y$ の直交射影であることを利用します。
\begin{align}
R &= \begin{Vmatrix}y - X\hat{b} \end{Vmatrix}^2 \\
&= \begin{Vmatrix}y - X(X'X)^{-1}X'y \end{Vmatrix}^2 \\
&= \begin{Vmatrix}(I - X(X'X)^{-1}X')y \end{Vmatrix}^2 \\
\\
& ここで I は単位行列. \\
&X(X'X)^{-1}X'y は y の ImX への直交射影であるので, X(X'X)^{-1}X'は二乗しても変わらない. \\
\\
&= y'(I - X(X'X)^{-1}X')y \\
&= y'y - y'X(X'X)^{-1}X'y
\end{align}
詳しくはぜひこちらを↓
ベイズらしくしていく
ここからは事前確率分布を導入し、よりベイズらしい扱いをしていきます。事前分布は次のような分布を使用します。（一様分布ではないことに注意）
P(b, v) \propto \frac{1}{v} \tag{3}
パラメータ v が今回のように連続の場合、定義域が有界でないなら v 上の積分が発散してしまい、正しく正規化できなくなります。こうした分布は、変則事前分布（不完全事前分布; improper prior）と呼ばれています。
ただ、実際計算される事後分布が正しく正規化されるなら、こうした事前分布を使用することはよくあります。
詳しくはPRML 2.4.3へ
線形回帰モデルのパラメータに変則事前分布を適用する
(1), (3) を使ってパラメータの事後分布を求めます。
P(b, v | y, X) \propto P(y|X,b,v)P(b, v) = \frac{P(y|X, b, v)}{v} = ( 2 \pi v )^{- \frac{n}{2}} \exp( -\frac{Q(y,b)}{2v} ) × \frac{1}{v}
この事後確率分布を $v$ について積分消去すると $P(b|y,X)$ となり パラメータ $b$ についての事後確率分布を求めることができます。
後半では、この事後確率分布から実際にデータをサンプリングして、どのような感じになるか見てみます。
βについての周辺分布を考える
では事後確率分布からパラメータ $v$ を積分消去してみましょう。
\begin{align}
P(b | y, X)
& \propto \int_{0}^{\infty} \frac{P(y|X,b,v)}{v} dv \\
&= \int_{0}^{\infty} ( 2 \pi )^{- \frac{n}{2}} v^{- \frac{n}{2}-1} \exp( -\frac{Q(y,b)}{2v} ) dv \\
\end{align}
ここで、 $t = \frac{Q(y,b)}{2v}$ とおくと、$\frac{dv}{dt} = - \frac{Q(y,b)}{2t^2}$ となります。
積分区間は、
| v | 0 → ∞ |
|---|---|
| t | ∞ → 0 |
となりますが、積分区間を反転させるのでマイナスがでてきます。このマイナスは置換するときに出てくるマイナスと打ち消し合います。
\begin{align}
& \Rightarrow \int_{0}^{\infty} ( 2 \pi )^{- \frac{n}{2}} v^{- \frac{n}{2}-1} \exp( -\frac{Q(y,b)}{2v} ) dv \\
&= \int_{0}^{\infty} ( 2 \pi )^{- \frac{n}{2}} \biggr( \frac{Q(y,b)}{2t} \biggr)^{- \frac{n}{2}-1} \exp( -t ) \frac{Q(y,b)}{2t^2} dt \\
&= ( 2 \pi )^{- \frac{n}{2}} \biggr( \frac{Q(y,b)}{2} \biggr)^{- \frac{n}{2}} \int_{0}^{\infty} t^{\frac{n}{2}-1} \exp( -t ) dt \\
&= \frac{ \pi }{2}^{- \frac{n}{2}} \Gamma{(\frac{n}{2})} Q(y,b)^{- \frac{n}{2}} \\
\end{align}
$Q(y,b) = (b - \hat{b})'X'X(b - \hat{b}) + R$ であることと、(2) の結果を利用すると、
\begin{align}
& \Rightarrow Q(y,b)^{- \frac{n}{2}}\\
&= \biggl[ (b - \hat{b})'X'X(b - \hat{b}) + R \biggr]^{- \frac{n}{2}} \\
& \propto \biggl[ 1 + \frac{(b - \hat{b})'X'X(b - \hat{b})}{R} \biggr]^{- \frac{n}{2}} \\
&= \biggl[ 1 + \frac{(b - \hat{b})'X'X(b - \hat{b})}{s^2(n-p)} \biggr]^{- \frac{n}{2}}
\end{align}
事後確率分布からパラメータ $v$ を積分消去した $b$ の分布は t分布 になることがわかります。
正規化定数は、多変数のt分布と見比べて作っていきましょう。
PRMLより。$\Lambda$ は精度行列、 $v$ は自由度で $D$ は次元数です。
\mathcal{St}(x|\mu, \Lambda, v) =
\frac{ \Gamma{(\frac{v}{2} + \frac{D}{2})} }{\Gamma{(\frac{v}{2})}}
\frac{ \begin{vmatrix}\Lambda \end{vmatrix}^{\frac{1}{2}} }{(v\pi)^{\frac{D}{2}}}
\biggl[ 1 + \frac{\Delta^2}{v} \biggr]^{ - \frac{v}{2} - \frac{D}{2}} \\
\\
\Delta^2 = (x-\mu)'\Lambda(x-\mu)
まずはガンマ関数の部分から見ていきましょう。
$\frac{v}{2} + \frac{D}{2} = \frac{n-p}{2} + \frac{p}{2} = \frac{n}{2}$ なので、
$$
\frac{ \Gamma{(\frac{v}{2} + \frac{D}{2})} }{\Gamma{(\frac{v}{2})}}
\Rightarrow
\frac{ \Gamma{(\frac{n}{2})} }{\Gamma{(\frac{n-p}{2})}}
$$
ここで $p$ はパラメータ数、すなわち次元数です。
次は行列式を考えましょう。$\Delta^2$ と $s^{-2}Q(y,b)$ を見比べて、
$$
(x-\mu)'\Lambda(x-\mu)
\Rightarrow
s^{-2}(b - \hat{b})'X'X(b - \hat{b})
$$
すなわち、$|\Lambda|^{\frac{1}{2}} \Rightarrow (s^2)^{-\frac{p}{2}}|X'X|^{\frac{1}{2}}$ です。(式変形の都合上約分しません)
残った $(v\pi)^{\frac{D}{2}}$ は $\bigl( (n-p)\pi \bigr)^{\frac{p}{2}}$ となります。
全て合わせてみましょう。
\begin{align}
P(b | y, X)
&=
\frac{ \Gamma{(\frac{n}{2})} }{\Gamma{(\frac{n-p}{2})}}
\frac{ (s^2)^{-\frac{p}{2}} |X'X|^{\frac{1}{2}} }{ \bigl( (n-p)\pi \bigr)^{\frac{p}{2}} }
\biggl[ 1 + \frac{(b - \hat{b})'X'X(b - \hat{b})}{s^2(n-p)} \biggr]^{- \frac{n}{2}} \\
&=
\frac{ \Gamma{(\frac{n}{2})} }{\Gamma{(\frac{n-p}{2})}}
\frac{ |X'X|^{\frac{1}{2}} }{ \bigl( s^2 \pi (n-p) \bigr)^{\frac{p}{2}} }
\biggl[ 1 + \frac{(b - \hat{b})'X'X(b - \hat{b})}{s^2(n-p)} \biggr]^{- \frac{n}{2}} \tag{4}
\end{align}
これで事後確率分布を正しく正規化できました。基本的にt分布など確率分布の積分は大変なので、正規化定数は見比べて作成することが多いらしいです。
（計算練習をしなくていいとは言ってないですよ！）
n -> 無限大の極限を考える
t分布の n -> 無限大の極限を考えてみましょう。自然対数の極限を利用します。
例）
\begin{align}
(1 + \frac{x^2}{n})^{-\frac{n}{2}}
&= \bigl( (1 + \frac{x^2}{n})^{\frac{n}{x^2}} \bigr) ^{-\frac{n}{2} \frac{x^2}{n}} \\
& n \Rightarrow ∞とすると \\
& (1 + \frac{x^2}{n})^{\frac{n}{x^2}} \Rightarrow e\\
&= \exp{(- \frac{x^2}{2}) }
\end{align}
この結果を利用して、(4) の正規化定数以外の部分に着目すると、
\begin{align}
\biggl[ 1 + \frac{(b - \hat{b})'X'X(b - \hat{b})}{s^2(n-p)} \biggr]^{- \frac{n}{2}}
&=
\Biggl( \biggl[ 1 + \frac{(b - \hat{b})'X'X(b - \hat{b})}{s^2(n-p)} \biggr]^{\frac{s^2(n-p)}{(b - \hat{b})'X'X(b - \hat{b})}} \Biggr)^{- \frac{n}{2} \frac{(b - \hat{b})'X'X(b - \hat{b})}{s^2(n-p)}} \\
& n \rightarrow \infty とすると \\
&= \exp{\Biggl( - \frac{(b - \hat{b})'X'X(b - \hat{b})}{2s^2}\Biggr)}
\end{align}
従って、$n \rightarrow \infty$ のときのt分布は、平均が $\hat{b}$ 分散が $s^2(X'X)^{-1}$ のガウス分布になることがわかります。
また、ガウス分布の分散がt分布のスケールになっていることも、この結果からわかります。
スケールは単位によって分布を横に引き伸ばしたり縮めたりする役目。
方眼紙の1マスの大きさについて言及していると考えるとわかりやすいかも。
ARモデルのパラメータに変則事前分布を適用する
ようやく本題のARモデルについてみていきます。ARモデルは線形モデルの少し特別な場合というだけなので、基本的な結果は変わりません。
では、さっそく進めていきましょう。簡単のために、ラグ1、定数0の1変数ARモデルを扱うこととします。
$$
y_t = \phi y_{t-1} + \epsilon_t \quad \epsilon ～ \mathcal{N}(0, v) \tag{5}
$$
線形モデルの場合と同じく尤度関数を作りましょう。観測値 $y_1, ... , y_t$ が手に入る確率は、
\begin{align}
P(y_t, y_{t-1}, ... , y_2 | y_1 , v)
&= \Pi_{t=1}^{T-1} \mathcal{N}( \phi y_{t}, v) \\
&= \biggl( \frac{1}{\sqrt{2 \pi v}} \biggr)^{ T - 1 } \exp{ \biggl(
- \frac{1}{2v} \Sigma_{t=1}^{T-1}( y_{t+1} - \phi y_{t} )^2 \biggr) }
\end{align}
この尤度関数を最大化させるパラメータ $\phi$ は、線形モデルの解を利用すると、
X = \begin{bmatrix}
y_1 \\
\vdots \\
y_{t-1} \\
\end{bmatrix}, \quad
y = \begin{bmatrix}
y_2 \\
\vdots \\
y_{t} \\
\end{bmatrix}
\begin{align}
\hat{\phi}
&= (X'X)^{-1}X'y \\
&= \frac{ \Sigma_{t=2}^{T} y_t y_{t-1} }{\Sigma_{t=1}^{T-1} y_t^2}
\end{align}
と表せます。もちろん尤度関数を微分しても構いません。同じ結果が出ると思います。
φについての周辺分布を考える
周辺分布の導出も線形モデルと同様に式変形をしてもいいのですが、せっかく行列で結果を出したので、それを利用して求めてみましょう。
まずは、残差平方和の分解を行っておきましょう。
\begin{align}
R &= y'y - y'X(X'X)^{-1}X'y \\
&= y'y - y'X \hat{b} \\
&= \Sigma_{t=2}^{T} y_t^2 -
\Sigma_{t=2}^{T} y_t y_{t-1}
\frac{ \Sigma_{t=2}^{T} y_t y_{t-1} }{\Sigma_{t=1}^{T-1} y_t^2} \\
&= \Sigma_{t=2}^{T} y_t^2 -
\frac{( \Sigma_{t=2}^{T} y_t y_{t-1} )^2 }{\Sigma_{t=1}^{T-1} y_t^2} \tag{6}
\end{align}
上の計算は普通に $\Sigma_{t=1}^{T-1}( y_{t+1} - \phi y_{t} )^2$ を展開して式変形したほうが、わかりやすいかもしれないです。
ここで一旦求めたい事後確率分布について少し整理します。
P(\phi, v | y_{1:T}) \propto P(y_{2:T}| y_1, \phi, v)P(\phi, v) = \frac{P(y_{2:T}|y_1, \phi, v)}{v}
$\phi$ についての周辺分布を求めたいので、$v$ を積分消去します。
\begin{align}
P(\phi | y_{1:T})
& \propto \int_{0}^{\infty} \frac{P(y_{2:T}|y_1, \phi, v)}{v} dv \\
\end{align}
この積分結果は、線形モデルで見たように__t分布__になります。
具体的に見ていきましょう。まずは平均ですが、これは最尤推定の解がそのままなので、$\hat{\phi} = \frac{ \Sigma_{t=2}^{T} y_t y_{t-1} }{\Sigma_{t=1}^{T-1} y_t^2}$ となります。
分散はスケールを自由度で調整したものです。スケールはガウス分布の分散に等しいことを利用します。(6)の結果も使うとt分布のスケールは
\begin{align}
s^2(X'X)^{-1}
&=
\biggl( \Sigma_{t=2}^{T} y_t^2 -
\frac{( \Sigma_{t=2}^{T} y_t y_{t-1} )^2 }{\Sigma_{t=1}^{T-1} y_t^2} \biggr)
\biggl( (n-2)\Sigma_{t=1}^{T-1} y_t^2 \biggr)^{-1} \\
&=
\biggl(
\frac{ \Sigma_{t=2}^{T} y_t^2 \Sigma_{t=2}^{T} y_{t-1}^2 -
\bigr( \Sigma_{t=2}^{T} y_t y_{t-1} \bigr)^2 }{\Sigma_{t=1}^{T-1} y_t^2} \biggr)
\biggl( (n-2)\Sigma_{t=1}^{T-1} y_t^2 \biggr)^{-1} \\
&=
\frac{ \Sigma_{t=2}^{T} y_t^2 \Sigma_{t=2}^{T} y_{t-1}^2 -
\bigr( \Sigma_{t=2}^{T} y_t y_{t-1} \bigr)^2 }
{(n-2) \bigl( \Sigma_{t=1}^{T-1} y_t^2 \bigr)^2}
\end{align}
となります。
1行目から2行目で通分を行っていますが、その時に分子の和の範囲を $1 ～ T-1$ から $2 ～ T$ の範囲に変えていることに注意しましょう。
以上より、事後確率分布のパラメータ $\phi$ についての周辺分布は
P(\phi | y_{1:T}) ～ \mathcal{t_{(t-2)}}( \hat{\phi}, \frac{n-2}{(n-2)-2}C(y_{1:n})) \\
\hat{\phi} = \frac{ \Sigma_{t=2}^{T} y_t y_{t-1} }{\Sigma_{t=1}^{T-1} y_t^2}
C(y_{1:n}) =
\frac{ \Sigma_{t=2}^{T} y_t^2 \Sigma_{t=2}^{T} y_{t-1}^2 -
\bigr( \Sigma_{t=2}^{T} y_t y_{t-1} \bigr)^2 }
{(n-2) \bigl( \Sigma_{t=1}^{T-1} y_t^2 \bigr)^2}
となります。
φについての周辺分布からサンプリングを行う
最後に周辺分布からサンプリングを行ってこの記事を終わろうと思います。Juliaでは、t分布の平均と分散をTDsitでは制御できないためLocationScaleを使います。
LocationScaleの使い方です。第一引数にロケーション、第二引数にスケール、第三引数に1変数の確率分布を取ります。要は平行移動させたり伸び縮みさせられるってことです。
LocationScale(μ,σ,ρ)
A location-scale transformed distribution with location parameter μ, scale parameter σ, and given univariate distribution ρ.
If 𝑍 is a random variable with distribution ρ, then the distribution of the random variable
𝑋=μ+σ𝑍
is the location-scale transformed distribution with location parameter μ and scale parameter σ.
If ρ is a discrete distribution, the probability mass function of the transformed distribution is given by
𝑃(𝑋=𝑥)=𝑃(𝑍=𝑥−μσ).
If ρ is a continuous distribution, the probability density function of the transformed distribution is given by
𝑓(𝑥)=1σρ(𝑥−μσ).
LocationScale(μ,σ,ρ) # location-scale transformed distribution
params(d) # Get the parameters, i.e. (μ, σ, and the base distribution ρ)
location(d) # Get the location parameter
scale(d) # Get the scale parameter
dist = LocationScale( -1, 1/2, Normal( 2, 2 ) )
mean(dist), var(dist) => # (0.0, 1.0)
using Plots
using Random
rng = MersenneTwister(20220127)
ϕ = 0.9
v = 100^0.5
N = 500
# AR(1)
f( x, ϕ, v ) = ϕ*x + rand( rng, Normal( 0, v ) )
X = zeros(N)
X[1] += f( rand(rng), ϕ, v )
for n in 2:N
X[n] += f( X[n-1], ϕ, v )
end
# テストデータのプロット
plot( X, label="AR(1)" )
# ロケーション
mu = X[1:end-1]'X[2:end] / sum(X[1:end-1].^2)
# スケール
C = ( sum(X[1:end-1].^2)sum(X[2:end].^2) - sum(X[1:end-1]'X[2:end])^2 ) / sum(X[1:end-1].^2)^2
# 事後確率分布からサンプリング
posterior = LocationScale( mu, C / ( N - 2 ), TDist( N - 2 ) )
posterior_sample = rand( rng, posterior, 5_000 );
# 周辺分布のプロット
histogram( posterior_sample, label="samples from P(ϕ|y)" )
scatter!( [mu], [0], label="ϕ_ML" )
$\phi$ の最尤推定解の周辺にサンプルが集中していることがわかりますね。
データの生成モデルは $\phi = 0.9$ なので、いい感じに推定できていると思います。
参考図書, ブログ等
図書
- Time Series: Modeling, Computation, and Inference
- 経済・ファイナンスデータの計量時系列分析
- Pattern Recognition and Machine Learning
サイト