はじめに

先日、Twitterを見ていたら、「StanとかJAGSは、MCMCを自力で書けない人が使うべきでない」というような指摘をするツイートを見つけました。それを見て、私自身「ギクッーー」となりました。この指摘は本当におっしゃる通りであり、Stanなんかは事前分布と尤度を書くだけで勝手にベイズ推定ができちゃいます。その他背後でどのような計算が行われているかは、やはり自分の手で書いてみないと何もわかりません。何が起こっているかわからないものを使うのは個人的にちょっと気持ち悪い感じがします。何でもかんでもStanにぶっこんでベイズベイズ～とか言ってましたが、確かに自分の手でMCMCを実装するべきだなと痛感しました。

そこで、この記事では前学期に受講したベイズの授業で扱った教科書である``Introduction to Bayesian econometrics" (Greenberg, 2011) の練習問題を基にして、自力でMCMCをしてみよう！という記事になります。

www.amazon.co.jp

ここでは何を扱うかといえば、通常のlinear regressionだけではあまり面白くないので、``Heteroskedasticity"を考慮したHeteroskedastic linear regressionをRで実装します。ちなみに、Greenberg (2011)の 8.4 Excersisesから8.1と8.3の練習問題に基づいています。今回は理論編です。実装編は以下から飛んでください。また、間違いや改善点等ございましたら是非指摘してください。

sunaninattahito.hatenablog.com

はじめに
Heteroskedasticityについて
Posteriorの導出
まとめ

Heteroskedasticityについて

まず、そもそもHeteroskedasticityとは何かについて説明していきます。通常のlinear regressionにおいては、誤差項 $\epsilon_i$ の分散である $Var(\epsilon_i)$ がすべての $i$ について一定であることを仮定します。しかし、実際に扱う（現実の）データは、すべての観察 $i$ において $Var(\epsilon_i)$ が一定であるとは限りません*1。

ここで、出てくるのがHeteroskedasticityの問題です。もし、分散が均一でないのなら、通常のlinear regressionによる推定はバイアスを持つことになります。そこで、このような問題に対処しようとするモデルが、今回実装しようとするHeteroskedastisic linear regressionです。

通常のlinear regressionは

$Y_i \sim \mathcal{N}(X_i\beta, σ^{2})$ .

といったgenerative modelですが、heteroskedastic linear regressionでは

$\begin{align*} Y_i &\sim \mathcal{N}(X_i\beta, \lambda_{i}^{-1}σ^{2}), \\ \lambda_i &\sim Ga\left(\frac{\nu}{2}, \frac{\nu}{2}\right), \ \ \text{where }\nu \text{ known} \end{align*}$

というような形で、表すことになります。ここでは、分散 $\sigma^{2}$ を $i$ ごとのパラメータである $\lambda_i^{-1}$ で補正するような形になってます。この $\lambda_i$ の導入によって、いわゆるheteroskedasticityの問題に対処します。ここで、 $\nu$ とは自由度（degree of freedom）を表し、また簡単のために自由度 $\nu$ は既知とします。

Posteriorの導出

次に、それぞれのパラメーターについてのposteriorを導出していきたいと思います。このセクションはあくまで私の理解度確認テストのようなものなので、「早くRで実装したいンだが」という方は実装編に飛んじゃってください。

sunaninattahito.hatenablog.com

Setting

まず、以下のようなモデルを仮定します。

\begin{align*} Y_i &= X_i\beta + \epsilon_i, \\ \text{where } &\epsilon_i \sim \mathcal{N}(0, σ^{2}). \end{align*}

ここで、 $X_i$ は $1 × K$ かつ 1 列目が切片である（=1）行列であり、 $\beta$ は $K × 1$ の係数の行列であるとします。

さらに、 $Y$ を $(Y_1, \cdots ,Y_n)$ となるような $N × 1$ 行列とし、 $X$ を観察$i$の全数と対応するように $N × K$ の行列とします。また、 $\epsilon$ を $(\epsilon_1, \cdots, \epsilon_n)$ となるような $N × 1$ 行列とします。そうした場合、上式を以下のように書き換えることができます。

$Y = X\beta + \epsilon$ .

ここで、Heteroskedastic errors を得るために、以下を仮定します。

$Y_i \sim \mathcal{N}(X_i\beta, \lambda_{i}^{-1}σ^{2})$ .

Likelihood

次に、Likelihood function (尤度関数) を見ていきます。ここで、$\Lambda$を$\lambda_i$についての対角行列であるとして

$$ \Lambda = \left( \begin{array}{cccc} \lambda_1& 0 & \ldots & 0 \\ 0 & \lambda_2 & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \ldots & \lambda_n \end{array} \right) $$

とします。また、$Y=(Y_1,\cdots, Y_n)$とすると、$\beta$, σ², $\lambda_i$を所与とした場合の尤度関数は

$$ f(Y|\beta, σ^2, \lambda_i) = \prod\limits_{i=1}^{n}\left(\frac{\lambda_i}{2\pi σ^2}\right)^{1/2} \exp\left[-\frac{1}{2σ^2}(Y-X\beta)^{T}\Lambda(Y-X\beta)\right]. $$

となります。

Priors

次にpriorを見ていきます。

βについて

ここで、$\beta$に以下の多変量正規分布のpriorを置きます。

$$ \beta \sim \mathcal{N}_K(\beta_0, B_0). $$

この時、$\beta$に関する事前分布の確率密度 $\pi (\beta)$は以下のように表すことができます。

$$ \pi(\beta) = \frac{1}{(2\pi)^{k/2}|B_0|^{1/2}} \exp\left[-\frac{1}{2}(\beta-\beta_0)^{T}B_0^{-1}(\beta-\beta_0)\right]. $$

σ²について

次に、σ²のpriorを見ていきます。ここでは、σ²に以下の逆ガンマ分布のpriorを置きます。

$$ σ^2 \sim \mathcal{IG}\left(\frac{\alpha_0}{2}, \frac{\delta_0}{2}\right). $$

この時、σ²に関する事前分布の確率密度 $\pi (σ^{2})$は以下のように表すことができます。

$$ \pi(σ^2) =\frac{(\delta_0 / 2)^{\alpha_0/2}}{\Gamma(\alpha_0/2)} \left(\frac{1}{σ^2}\right)^{(\alpha_0/2) + 1} \exp\left[-\frac{\delta_0}{2σ^2}\right]. $$

λiについて

ここでは、$\lambda_i$に以下のガンマ分布のpriorを置きます。

$$ \lambda_i \sim Ga\left(\frac{\nu}{2}, \frac{\nu}{2}\right), \ \ \nu \text{ known}. $$

ここで、自由度$\nu$は既知とします。この時、$\lambda_i$に関する事前分布の確率密度$\pi(\lambda_i)$は

$$ \pi(\lambda_i) = \frac{(\nu/2)^{\nu/2}}{\Gamma(\nu/2)}\lambda_i^{(v/2)-1} \exp\left[-\sum\limits_{i=1}^{n}\frac{\nu\lambda_i}{2}\right]. $$

となります。また、$\pi(\lambda) = \pi(\lambda_1, \cdots, \lambda_n)$と定めると

$$ \begin{align} \pi(\lambda) &= \prod\limits_{i=1}^{n}\pi (\lambda_i)\\ &= \frac{(\nu/2)^{\nu/2}}{\Gamma(\nu/2)} \prod\limits_{i=1}^{n}\lambda_i^{(v/2)-1} \exp\left[-\sum\limits_{i=1}^{n}\frac{\nu\lambda_i}{2}\right]. \end{align} $$

として表します。

Posterior

最後にPosteriorを導出していきます。これはめちゃくちゃダルいので心が折れそうになりますが頑張ります。

Joint posterior

Joint posterior $\pi (\beta, σ^{2}, \lambda|Y)$ はベイズの定理より以下のように書き表すことができます。

$$ \begin{align} \pi(\beta, σ^2, \lambda|Y) &\propto \pi(\beta)\pi(σ^2)\pi(\lambda)f(Y|\beta, σ^2, \lambda)\\ &= \frac{1}{(2\pi)^{k/2}|B_0|^{1/2}} \exp\left[-\frac{1}{2}(\beta-\beta_0)^{T}B_0^{-1}(\beta-\beta_0)\right]\\ &\times \frac{(\delta_0 / 2)^{\alpha_0/2}}{\Gamma(\alpha_0/2)} \left(\frac{1}{σ^2}\right)^{(\alpha_0/2) + 1} \exp\left[-\frac{\delta_0}{2σ^2}\right]\\ &\times \frac{(\nu/2)^{\nu/2}}{\Gamma(\nu/2)} \prod\limits_{i=1}^{n}\lambda_i^{(v/2)-1} \exp\left[-\sum\limits_{i=1}^{n}\frac{\nu\lambda_i}{2}\right]\\ &\times \prod\limits_{i=1}^{n}\left(\frac{\lambda_i}{2\pi σ^2}\right)^{1/2} \exp\left[-\frac{1}{2σ^2}(Y-X\beta)^{T}\Lambda(Y-X\beta)\right] \end{align} $$

Conditional posterior of σ²

ここから先は、上記のjoint posteriorからそれぞれのconditional posteriorを求めることになります。まずは、導出のし易いσ²から見ていきます。σ²のconditonal posteriorを求めたいときは、上記のjoint posteriorからσ²に関する項のみに着目します。この時

$$ \begin{align} \pi(σ^2|\beta, \lambda, Y) &\propto \left(\frac{1}{σ^2}\right)^{(\alpha_0/2) + 1} \exp\left[-\frac{\delta_0}{2σ^2}\right]\\ &\times \left(\frac{1}{σ^2}\right)^{n/2} \exp\left[-\frac{1}{2σ^2}(Y-X\beta)^{T}\Lambda(Y-X\beta)\right]\\ &= \left(\frac{1}{σ^2}\right)^{[(\alpha_0/ + n)/ 2] + 1}\exp\left[-\frac{1}{2σ^2}\{(\delta_0 + Y-X\beta)^{T}\Lambda(Y-X\beta)\}\right]. \end{align} $$

となります。また、

$$ \begin{align} \alpha_1 &= \alpha_0 + n, \\ \delta_1 &= \delta_0 + (Y-X\beta)^{T}\Lambda(Y-X\beta). \end{align} $$

とすると、上式は

$$ \begin{align} \pi(σ^2|\beta, \lambda, Y) &\propto \left(\frac{1}{σ^2}\right)^{(\alpha_1/2) + 1} \exp\left[-\frac{\delta_1}{2σ^2}\right] \end{align} $$

と書き換えられます。この時、σ²のconditional posteriorは逆ガンマ分布からのサンプリングとなります。よって、

$$ \begin{align} &σ^2 | \beta, \lambda, Y \sim \mathcal{IG} \left(\frac{\alpha_1}{2}, \frac{\delta_1}{2}\right)\\ \text{where }\ &\alpha_1 = \alpha_0 + n, \\ &\delta_1 = \delta_0 + (Y-X\beta)^{T}\Lambda(Y-X\beta). \end{align} $$

Conditional posterior of λi

Joint posteriorの式より、$\lambda_i$に関する式を抜き出します。$\lambda_i$のconditional posteriorは

$$ \begin{align} \pi(\lambda_i|\beta, σ^2, Y_i) &\propto \lambda_i^{(v/2)-1} \exp\left[-\frac{\nu\lambda_i}{2}\right]\\ &\times \lambda_i^{1/2}\exp\left[-\frac{\lambda_i}{2σ^2}(Y_i - X_i\beta)^2\right]\ \ \cdots\ \ \text{Likelihoodの}Y, X\text{を}i\text{についての式に書き換えるとこうなる。} \\ &= \lambda_i^{[(\nu+1)/ 2]-1} \exp\left[-\frac{\lambda_i}{2}\{\nu + σ^{-2} (Y_i - X_i\beta\}\right]. \end{align} $$

また、

$$ \begin{align} \nu_1 &= \nu + 1\\ \nu_{2, i} &= \nu + σ^{-2}(Y_i - X_i\beta) \end{align} $$

とすると、上式は

$$ \pi(\lambda_i|\beta, σ^2, Y_i) \propto \lambda_i^{[(\nu_1)/ 2]-1} \exp\left[-\frac{\nu_{2, i}\lambda_i}{2}\right] $$

と書き換えられます。この時、$\lambda_i$のconditional posteriorはガンマ分布からのサンプリングとなります。よって

$$ \begin{align} &\lambda_i | \beta, σ^2, Y_i \sim Ga\left(\frac{\nu_1}{2}, \frac{\nu_{2, i}}{2}\right), \\ \text{where }\ \ &\nu_1 = \nu + 1,\\ &\nu_{2, i} = \nu + σ^{-2}(Y_i - X_i\beta). \end{align} $$

Conditional posterior of β

最後です。丁寧に解説するのはかなり面倒なので、ところどころ端折ってやります。Joint posteriorより、$\beta$のconditional posteriorは

$$ \begin{align} \pi(\beta| \lambda, σ^2, Y) &\propto \exp\left[-\frac{1}{2}(\beta-\beta_0)^{T}B_0^{-1}(\beta-\beta_0)\right]\\ &\times \exp\left[-\frac{1}{2σ^2}(Y-X\beta)^{T}\Lambda(Y-X\beta)\right]\\ &= \exp \left[-\frac{1}{2}\{(\beta-\beta_0)^{T}B_0^{-1}(\beta-\beta_0) + σ^{-2}(Y-X\beta)^{T}\Lambda(Y-X\beta)\}\right] \end{align} $$

となります。また、上式の$\exp$内を色々やって整理すると

$$ \begin{align} \pi(\beta| \lambda, σ^2, Y) &\propto \exp\left[-\frac{1}{2}\{(\beta-\tilde{\beta})^{T}B_1^{-1}(\beta-\tilde{\beta})\}\right]\\ \text{where}\ \ B_1 &= [B_0^{-1} + σ^{-2}X^{T}\Lambda X]^{-1}, \\ \tilde{\beta} &= B_1[B_0^{-1}\beta_0 + σ^{-2}X^{T}\Lambda Y]. \end{align} $$

となる。この時、$\beta$のconditional posteriorは正規分布からのサンプリングとなります。よって

$$ \begin{align} &\beta | \lambda, σ^2, Y \sim \mathcal{N}_K(\tilde{\beta}, B_1), \\ \text{where}\ \ &B_1 = [B_0^{-1} + σ^{-2}X^{T}\Lambda X]^{-1}, \\ &\tilde{\beta} = B_1[B_0^{-1}\beta_0 + σ^{-2}X^{T}\Lambda Y]. \end{align} $$

まとめ

ここまでで求められたposteriorを以下で並べてみましょう。

$$ \begin{align} \beta|\lambda, σ^2, Y &\sim \mathcal{N}_K(\tilde{\beta}, B_1), \\ σ^2|\beta, \lambda, Y &\sim \mathcal{IG}\left(\frac{\alpha_1}{2}, \frac{\delta_1}{2}\right), \\ \lambda_i|\beta, σ^2, Y_i &\sim \textit{Ga}\left(\frac{\nu_1}{2}, \frac{\nu_{2i}}{2}\right), \\ \\ \text{where}\ \ B_1 &= [B_0^{-1} + σ^{-2}X^{T}\Lambda X]^{-1}, \\ \tilde{\beta} &= B_1[B_0^{-1}\beta_0 + σ^{-2}X^{T}\Lambda Y],\\ \alpha_1 &= \alpha_0 + n, \\ \delta_1 &= \delta_0 + (Y-X\beta)^{T}\Lambda (Y-X\beta), \\ \nu_1 &= \nu + 1, \\ \nu_{2,i} &= \nu + σ^{-2}(Y_i-X_i\beta). \end{align} $$

どうでも良い話ですが、はてなブログ上でtex記法を使って\sigmaを打とうとすると、なぞにキーワードにリンクされてしまい、数式がうまく表示されないのでクソすぎると思いました。

*1:実際には、Breusch-Pagan testやWhite testを用いてこのHeteroskedasticityを確認します。

砂になった人

砂になった人の趣味や勉強についての備忘録

【R】Bayesian heteroskedastic linear regresssionをRで実装してみる（理論編）

はじめに

Heteroskedasticityについて

Posteriorの導出

Setting

Likelihood