線形回帰モデル
前回DeepLearningを0から実装できるようにするための連載 Vol.0では確率モデルの最尤推定について説明しました。今回は最も基本的でかつ重要な確率モデルである線形回帰モデルを扱います。このモデルに対して最尤推定を行うことでおなじみの最小二乗法が出てきます。
モデル定義
線形予測
線形モデルではまず変数を２種類$(x, y)$にわけます。１つ目の$x$は独立変数と言い、２つ目の$y$は従属変数と言います。ここでは$x$が$D$次元で$y$が１次元として展開しますが、$y$が多次元の場合でも同様の議論が成立します。線形回帰モデルの目的は$y$を$x$の線型結合で予測することです。
y^{pred} = \sum_i w_ix_i + b
ここで$w$を系数、$b$を切片と言います。$y^{pred}$は$y$の予測値です。$x$の定義を$(x_{\cdot 1}, x_{\cdot 2}, ...,x_{\cdot D}, 1)$のように最後に１を追加し、$w$も$(w_1, w_2, ...,w_D, b)$のように末尾にbを追加することでよりコンパクトに
y^{pred} = \sum_i w_ix_i = w^T x
のようにベクトルの内積を用いて表現できます。そしてこの係数$w$を観測したデータから推定することが線形回帰モデルのパラメータを学習することです。
確率を導入する。
線形回帰モデルは確率モデルですが、ここまで確率がまだ登場していません。先ほど$y^{pred}$と書いたことには意味があり、ここで確率が登場します。線形回帰モデルでは実際に観測された$y$は$y^{pred}$を平均として正規分布に従っていると仮定します。
\begin{align}
P(y|y^{pred}) & = \mathcal{N(y|y^{pred}, \sigma^2)} \\
& = \frac{1}{\sqrt{2\pi}\sigma}\mathrm{exp}\left( \frac{(y-y^{pred})^2}{2\sigma^2}\right) \\
& = \frac{1}{\sqrt{2\pi}\sigma}\mathrm{exp}\left( \frac{(y-w^Tx)^2}{2\sigma^2}\right)
\end{align}
これでモデルの準備が整いました。
最尤推定
$N$つの独立なデータ$(X, Y) = [(x_n,y_n), n=1,...,N]$を観測した際の尤度は次のように単純な積で表現出来ました。
P(Y |X, w, \sigma) = \prod_n \frac{1}{\sqrt{2\pi}\sigma}\mathrm{exp}\left( \frac{(y_n-w^Tx_n)^2}{2\sigma^2}\right)
これを最大にするような$w$を探します。まずはNegetive Log-Likelihood(NNL)を計算します。
\begin{align}
E(w, \sigma) & = -\ln \prod_n \frac{1}{\sqrt{2\pi}\sigma}\mathrm{exp}\left( \frac{(y_n-w^Tx_n)^2}{2\sigma^2}\right) \\
& = \frac{1}{2}\sum_n\left[ \ln{2\pi} + \ln \sigma^2 + \frac{1}{\sigma^2}(y_n - w^Tx_n)^2 \right] \\
& = \frac{1}{2}\sum_n\left[ \ln{2\pi} + \ln \sigma^2 \right] + \frac{1}{\sigma^2}E(w)
\end{align}
ここで
E(w) = \frac{1}{2} \sum_n (y_n - w^Tx_n)^2 = \frac{1}{2}\sum_n (y_n - y^{pred}_n)^2
が$w$に依存する部分です。最小二乗法の式が登場しました。これを最小にするには$E(w)$を$w$で微分して勾配を求めます。
\nabla E(w) = \sum_n (w^T x_n - y_n)x_n
となります。紙とペンで各自確認してみてください。あとは勾配降下法で$w$を逐次更新していけばOKです。
w^{t+1} = w^t - \alpha \nabla E(w^t)
正規方程式
線形回帰モデルは非常にシンプルなモデルで実は勾配降下法で反復的に解く以外にも行列計算で一度に解くこともできます。まずは勾配の式がゼロになる連立方程式を立てます。
0 = \sum_n (w^T x_n - y_n)x_n
これを変形すると次のようになります。
w^T \left( \sum_n x_nx_n^T \right) = \sum_n x_ny_n
行列表記を用いると
X^T Xw = X^T Y
のようになります。これを正規方程式と言います。正規方程式をLAPACKなどのライブラリで解いてもいいのですが、更に変形を進めるとより直鉄的な式が得られます。
w_{ML} = \left( X^T X \right) ^{-1}X^T Y = X^{\dagger}Y
ここで$X^{\dagger}$はXの一般化逆行列(pseudo-inverse)と言います。$X$の次元は$(N,D)$で必ずしも正方行列ではないので逆行列も必ずしも存在しませんが、一般化逆行列は常に存在します。一般化逆行列を計算するときに上の式と通り計算すると速度や精度に問題のある逆行列の計算が必要なので通常は特異値分解(SVD)を用いた方法を使用します。(ここではこれ以上説明しません。)
分散の推定
線形回帰モデルのパラメータは$w$以外に分散パラメータの$\sigma$もありました。こちらは$E(w, \sigma)$に$w_{ML}$を代入し、$\sigma$で微分をすると以下の推定式が得られます。
\sigma ^2 _{ML} = \frac{1}{N} \sum_n (y_n - w_{ML}^T x_n)
多次元出力への拡張
$y = (y_1, ..., y_K)$のように従属変数が多次元の場合へのモデル拡張は非常にシンプルです。
\begin{align}
y_{\cdot k}^{pred} & = w_k^T x \\
E(w) & = \frac{1}{2} \sum_n \sum_k (y_{nk} - w_k^Tx_n)^2 = \frac{1}{2}\sum_n \sum_k (y_{nk} - y^{pred}_{nk})^2 \\
\nabla _{w_k} E(w) & = \sum_n (w_k^T x_n - y_{nk})x_n
\end{align}
係数の更新式は以下のようになります。
w_k^{t+1} = w_k^t - \alpha \nabla _{w_k} E(w^t)
次回予告
線形回帰は従属変数が連続変数なので回帰分析のための手法です。次回は線形回帰と同様に線形モデルを使用する分類のための手法であるロジスティック回帰(Logistic Regression)をやります。(回帰と名前がついていますがクラス分類の手法です。)