2020-07-16

線形判別分析を理解する

線形判別分析

線形判別分析(LDA)は各クラスの説明変数の分布が正規分布に従っていること、また各クラスの共分散行列が等しいと仮定した特別な場合です。説明変数が複数ある場合は多変量正規分布になりますが、多変量正規分布の確率密度関数は

$\displaystyle{ \begin{align} f(\boldsymbol{x}) = \frac{1}{(2 \pi)^{p/2} |\boldsymbol{\Sigma}|^{1/2}} \exp \left\{- \frac{1}{2} (\boldsymbol{x} - \boldsymbol{\mu})^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} (\boldsymbol{x} - \boldsymbol{\mu})\right\} \end{align} }$

です。ここで $\boldsymbol{x}$ は $p$ 次元（説明変数の数が $p$ ）で、説明変数の平均が $\boldsymbol{\mu}$ （説明変数ごとに平均値があるためベクトル）で $\boldsymbol{\Sigma}$ は共分散行列です。

では線形判別分析によるクラス分類を考えます。各クラスの説明変数の分布が多変量正規分布に従っているとしたとき、クラス $k$ の多変量正規分布のパラメータを $(\boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)$ と書くことにします。では、分類のためクラス $k$ の事後確率を求めます。これはベイズの定理より以下のように求めます。

$\displaystyle{ \begin{align} \mathrm{P}(G = k | \boldsymbol{X} = \boldsymbol{x}) = \pi_k \frac{f_k(\boldsymbol{x})}{\sum^K_{l=1}f_l(\boldsymbol{x}) \pi_l} \end{align} }$

各項の解釈はこちらになります。

$\mathrm{P}(G = k | \boldsymbol{X} = \boldsymbol{x})$ : 事後確率。データ $\boldsymbol{x}$ が観測されたもとでのクラス $k$ が生じる確率。
$\pi_k$ : 事前確率。クラス $k$ が生じる確率。
$f_k(\boldsymbol{x})$ : 尤度。 $\mathrm{P}(\boldsymbol{X} = \boldsymbol{x} | G = k)$ 。クラス $k$ が与えられたもとでデータ $\boldsymbol{x}$ が得られる確率密度。
$\sum^K_{l=1}f_l(\boldsymbol{x}) \pi_l$ ：周辺化された尤度。 $P(\boldsymbol{X} = \boldsymbol{x})$ 。データ $\boldsymbol{x}$ が得られる確率。 $K$ はクラス数（2クラスであれば2）。

この(2)式を使って各クラスの確率を比較してみます。今考えている問題は各クラスの確率です。そのため、クラス（ここではクラス $k$ ）に依存しない項は意味を持たないため、無視することにします。分母の $\sum^K_{l=1}f_l(\boldsymbol{x}) \pi_l$ はデータが得られる確率であり、クラスに依存しないため無視します。比例記号を使って

$\displaystyle{ \mathrm{P}(G = k | \boldsymbol{X} = \boldsymbol{x}) \propto \pi_k \cdot f_k(\boldsymbol{x}) }$

と表します。つぎに対数を取り $\delta_k(\boldsymbol{x})$ と書くことにします。

$\displaystyle{ \delta_k(\boldsymbol{x}) = \log \{ \pi_k \cdot f_k(\boldsymbol{x}) \} }$

多変量正規分布の(1)式を代入し、変形していきます。

$\displaystyle{ \begin{align} \delta_k(\boldsymbol{x}) &= \log \{ \pi_k \cdot f_k(\boldsymbol{x}) \} \nonumber \\ &= \log \pi_k + \log f_k(\boldsymbol{x}) \nonumber \\ &= \log \pi_k + \log \{(2\pi)^{p/2} |\boldsymbol{\Sigma}_k^{1/2}|\}^{-1} - \frac{1}{2} (\boldsymbol{x} - \boldsymbol{\mu}_k)^{\mathrm{T}} \boldsymbol{\Sigma}_k^{-1} (\boldsymbol{x} - \boldsymbol{\mu}_k) \nonumber \\ &= \log \pi_k - \frac{p}{2} \log (2\pi) - \frac{1}{2} \log |\boldsymbol{\Sigma}_k| - \frac{1}{2} (\boldsymbol{x} - \boldsymbol{\mu}_k)^{\mathrm{T}} \boldsymbol{\Sigma}_k^{-1} (\boldsymbol{x} - \boldsymbol{\mu}_k) \end{align} \nonumber }$

$\frac{p}{2}\log(2\pi)$ も各クラス共通に現れるので無視します。

$\displaystyle{ \delta_k(\boldsymbol{x}) = \log \pi_k - \frac{1}{2} \log |\boldsymbol{\Sigma}_k| - \frac{1}{2} (\boldsymbol{x} - \boldsymbol{\mu}_k)^{\mathrm{T}} \boldsymbol{\Sigma}_k^{-1} (\boldsymbol{x} - \boldsymbol{\mu}_k) }$

線形判別分析は各クラスが共通の共分散行列を持つと仮定した場合の手法です。なので、 $\boldsymbol{\Sigma}_k = \boldsymbol{\Sigma}$ とします。

$\displaystyle{ \delta_k(\boldsymbol{x}) = \log \pi_k - \frac{1}{2} \log |\boldsymbol{\Sigma}| - \frac{1}{2} (\boldsymbol{x} - \boldsymbol{\mu}_k)^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} (\boldsymbol{x} - \boldsymbol{\mu}_k) }$

こうすると $\frac{1}{2} \log |\boldsymbol{\Sigma}|$ も各クラス共通になるので無視します。

$\displaystyle{ \delta_k(\boldsymbol{x}) = \log \pi_k - \frac{1}{2} (\boldsymbol{x} - \boldsymbol{\mu}_k)^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} (\boldsymbol{x} - \boldsymbol{\mu}_k) }$

第2項を展開します。

$\displaystyle{ \delta_k(\boldsymbol{x}) = \log \pi_k - \frac{1}{2}(\boldsymbol{x}^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} \boldsymbol{x} - \boldsymbol{x}^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu}_k - \boldsymbol{\mu}_k^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} \boldsymbol{x} + \boldsymbol{\mu}_k^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu}_k) }$

ここで、 $\boldsymbol{\mu}_k^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} \boldsymbol{x}$ はスカラなので、転置をとっても同じ結果になるので、

$\displaystyle{ \begin{align} \delta_k(\boldsymbol{x}) &= \log \pi_k - \frac{1}{2}(\boldsymbol{x}^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} \boldsymbol{x} - \boldsymbol{x}^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu}_k - (\boldsymbol{\mu}_k^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} \boldsymbol{x})^\mathrm{T} + \boldsymbol{\mu}_k^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu}_k) \nonumber \\ &= \log \pi_k - \frac{1}{2}(\boldsymbol{x}^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} \boldsymbol{x} - \boldsymbol{x}^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu}_k - \boldsymbol{x}^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu}_k + \boldsymbol{\mu}_k^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu}_k) \nonumber \\ &= \log \pi_k - \frac{1}{2}\boldsymbol{x}^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} \boldsymbol{x} + \boldsymbol{x}^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu}_k - \frac{1}{2} \boldsymbol{\mu}_k^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu}_k \nonumber \\ \end{align} }$

$- \frac{1}{2}\boldsymbol{x}^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} \boldsymbol{x}$ もクラス共通のため無視します。

$\displaystyle{ \begin{align} \delta_k(\boldsymbol{x}) = \log \pi_k + \boldsymbol{x}^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu}_k - \frac{1}{2} \boldsymbol{\mu}_k^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu}_k \end{align} }$

この式を線形判別関数といいます。この式の計算結果が最大となるクラスに分類します。

線形判別関数のパラメータ

線形判別関数のパラメータですが、(4)式を眺めると入力データである $\boldsymbol{x}^\mathrm{T}$ 以外の( $\pi_k$ , $\boldsymbol{\mu}_k$ , $\boldsymbol{\Sigma}$ )です。これらを訓練データを基に推定（学習）する必要があります。と言ってもとくに難しいことではなく、

$\displaystyle{ \begin{align} \hat{\pi}_k &= \frac{N_k}{N} \nonumber \\ & (N_k\mathrm{はクラス}k\mathrm{の数}) \nonumber\\ \hat{\boldsymbol{\mu}}_k &= \frac{1}{N_k} \sum_{g_i=k} \boldsymbol{x}_i \nonumber \\ & (\mathrm{クラス}k\mathrm{の説明変数ごとの平均値}) \nonumber \\ \hat{\boldsymbol{\Sigma}} &= \frac{1}{N - K} \sum^K_{k=1} \sum_{g_i=k} (\boldsymbol{x}_i - \hat{\boldsymbol{\mu}}_k)(\boldsymbol{x}_i - \hat{\boldsymbol{\mu}}_k)^\mathrm{T} \nonumber \\ & (\mathrm{クラスごとの共分散行列の平均値}) \nonumber \end{align} }$

で推定できます。

対数比による分類

特に2クラス分類のとき、事後確率の対数比から容易に分類が行なえます。

クラス1とクラス2を分類することを考えた場合、それらの事後確率の対数比は、

$\displaystyle{ \log \frac{\mathrm{P}(G=2 | \boldsymbol{X} = \boldsymbol{x})}{\mathrm{P}(G=1 | \boldsymbol{X} = \boldsymbol{x})} \nonumber }$

です。確率比の対数を取る理由ですが、それは確率比のパターンを分けて考えると理解できます。確率比は以下の3パターンです。

分子が大きいケース： $\mathrm{P}(G=2 | \boldsymbol{X} = \boldsymbol{x}) \gt \mathrm{P}(G=1 | \boldsymbol{X} = \boldsymbol{x})$
分母が大きいケース： $\mathrm{P}(G=2 | \boldsymbol{X} = \boldsymbol{x}) \lt \mathrm{P}(G=1 | \boldsymbol{X} = \boldsymbol{x})$
分子と分母が等しいケース： $\mathrm{P}(G=2 | \boldsymbol{X} = \boldsymbol{x}) = \mathrm{P}(G=1 | \boldsymbol{X} = \boldsymbol{x})$

それぞれのパターンについて対数を取ると以下になります。

$\log \frac{\mathrm{P}(G=2 | \boldsymbol{X} = \boldsymbol{x})}{\mathrm{P}(G=1 | \boldsymbol{X} = \boldsymbol{x})} \gt 0$
$\log \frac{\mathrm{P}(G=2 | \boldsymbol{X} = \boldsymbol{x})}{\mathrm{P}(G=1 | \boldsymbol{X} = \boldsymbol{x})} \lt 0$
$\log \frac{\mathrm{P}(G=2 | \boldsymbol{X} = \boldsymbol{x})}{\mathrm{P}(G=1 | \boldsymbol{X} = \boldsymbol{x})} = 0$

確率が高い方のクラスに分類を行いますので、つまり符号によって分類が行えることになります。なお2クラスの確率が等しい場合は0になりますので、それが決定境界を表します。

上で求めた線形判別関数(4)式を使えば

$\displaystyle{ \begin{align} \log \frac{\mathrm{P}(G=2 | \boldsymbol{X} = \boldsymbol{x})}{\mathrm{P}(G=1 | \boldsymbol{X} = \boldsymbol{x})} &= \log \mathrm{P}(G=2 | \boldsymbol{X} = \boldsymbol{x}) - \log \mathrm{P}(G=1 | \boldsymbol{X} = \boldsymbol{x}) \nonumber \\ &= \log \pi_2 + \boldsymbol{x}^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu}_2 - \frac{1}{2} \boldsymbol{\mu}_2^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu}_2 - \left\{ \log \pi_1 + \boldsymbol{x}^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu}_1 - \frac{1}{2} \boldsymbol{\mu}_1^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu}_1 \right \} \nonumber \\ &= \log \frac{\pi_2}{\pi_1} + \boldsymbol{x}^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu}_2 - \boldsymbol{x}^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu}_1 - \frac{1}{2} \boldsymbol{\mu}_2^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu}_2 + \frac{1}{2} \boldsymbol{\mu}_1^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu}_1 \nonumber \\ &= \log \frac{\pi_2}{\pi_1} + \boldsymbol{x}^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} (\boldsymbol{\mu}_2 - \boldsymbol{\mu}_1) - \frac{1}{2} ( \boldsymbol{\mu}_2 + \boldsymbol{\mu}_1)^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} (\boldsymbol{\mu}_2 - \boldsymbol{\mu}_1) \end{align} \nonumber }$

となります。あとはこの結果の符号を見て判別が行なえます。なお補足になりますが、線形判別関数を求める際に、クラス共通のため無視してきた項というのは1行目が表しているように、引き算により消える運命にあります。

2次判別分析

2次判別分析(QDA)は、線形判別分析で仮定した「各クラスの共分散行列が等しい」を仮定しない場合です。よって、仮定をする前の(3)式である以下が2次判別関数です。

この場合、 $- \frac{1}{2}\boldsymbol{x}^{\mathrm{T}} \boldsymbol{\Sigma}_k^{-1} \boldsymbol{x}$ が消えず、これは $\boldsymbol{x}$ の2次関数であるため、決定境界が2次式となります。

線形判別分析による次元削減

線形判別分析は主成分分析のように次元削減という特徴も持っています。これは別の機会に記事にします。

参考

統計的学習の基礎

2020-07-07

バイアスとバリアンスを理解する

よくバイアスがあると言うことがありますが、これが一体何を意味しているのかを説明するというのが今回の意図です。データをもとに母数を言い当てることを推定といいます。ちなみに母数とは、正規分布における平均値や分散などそのモデルが持っている特性値です。仮にデータが正規分布をもとに生成されるデータであるとすれば、平均値と分散がわかればそのデータ生成プロセスを解明できることになります。ただ、一般には大本の（母集団）の平均値や分散がわかるということはあまりなく、得られたデータからそれらを言い当てる（推定する）ということになります。当然、完璧に言い当てられるわけではありませんから「誤差」を評価する必要が出てくるわけです。ここで問題になるのがどういった推定量（推定した値）が良いのでしょうか。なるべく誤差が小さければよいというのは当たり前ですが、ちょっとあやふやそうです。というのも、推定量も確率変数なのでなんらかの分布になるのです。一方で母数は定数です。推定量は確率変数なので、期待値や分散が計算できます。よって、これらの統計量が推定量の良さを計るヒントになりそうです。

母数が定数で、推定量が確率変数であるというところを直感的にわかるよう示した図です。この場合、どちらの推定量の分布が良いと言えるでしょうか。はい右側です。推定量の分布の中心が母数に近いからです。分布の中心は期待値を使って数値的に表現できます。母数を $\theta$ 、推定量を $\hat{\theta}$ としたとき、母数と推定量の期待値の差が小さければ良いので、

$\displaystyle{ \begin{align} \mathrm{Bias}(\hat{\theta}) = \mathrm{E}[\hat{\theta}(\boldsymbol{X})]- \theta \end{align} }$

と書き、これをバイアスと言います。なので、バイアスとは真の値との平均的なズレのことです。このバイアスが0であるような場合、つまり

$\displaystyle{ \begin{align} \mathrm{E}[\hat{\theta}(\boldsymbol{X})] = \theta \end{align} }$

のとき、推定量の分布の中心と母数が一致していますので、いわゆる偏りがありません。偏りが無いので不偏性を持つと言います。また、不偏性を持つ推定量を不偏推定量と言います。不偏であるとは推定量が平均的に母数の周りに分布していることを意味します。不偏性を持てばすべてハッピーかと言えばそうでもありません。推定量の散らばりが大きすぎればダメです。ということで、散らばりを評価するため推定量の分散も必要になります。これは分散の定義どおりに、

$\displaystyle{ \begin{align} \mathrm{Var}(\hat{\theta}) = \mathrm{E}[\{\hat{\theta}(\boldsymbol{X}) - \mathrm{E}[\hat{\theta}(\boldsymbol{X})]\}^2] \end{align} }$

になります。これを分散、または英語のバリアンスと言います。これら2つを推定量の良さを計る指標をすれば良いことになります。が、できれば1つの式にまとめたいところではあります。ちょっとだけ見方を変えて、推定量と母数の差を評価することを考えてみます（念のためですが、バイアスは推定量の期待値と母数の差を評価していますので違います。）。これはいわゆる誤差の評価になりますので、評価方法はいろいろとあると思いますが、ここでは平均2乗誤差(MSE)を用います。MSEは真の値と予測値の差の2乗の平均値です。推定量は確率変数なので期待値という表現を使うと、

$\displaystyle{ \begin{align} \mathrm{MSE}(\theta, \hat{\theta}) = \mathrm{E}[\{\hat{\theta}(\boldsymbol{X}) - \theta\}^2] \end{align} }$

になります。この式を変形してみます。

$\displaystyle{ \begin{align} \mathrm{MSE}(\theta, \hat{\theta}) &= \mathrm{E}[\{\hat{\theta}(\boldsymbol{X}) - \theta\}^2] \nonumber \\ &= \mathrm{E}[ \hat{\theta}^2(\boldsymbol{X}) - 2\theta\hat{\theta}^2(\boldsymbol{X}) + \theta^2] \nonumber \\ &= \mathrm{E}[\hat{\theta}^2(\boldsymbol{X})] - 2\theta\mathrm{E}[\hat{\theta}(\boldsymbol{X})] + \theta^2 \nonumber \\ &= \mathrm{E}[\hat{\theta}^2(\boldsymbol{X})] - 2\theta\mathrm{E}[\hat{\theta}(\boldsymbol{X})] + \{\mathrm{E}[\hat{\theta}(\boldsymbol{X})] - \theta\}^2 \nonumber \\ & \space \space \space + 2\theta\mathrm{E}[\hat{\theta}(\boldsymbol{X})] - \{\mathrm{E}[\hat{\theta}(\boldsymbol{X})]\}^2 \nonumber \\ &= \{\mathrm{E}[\hat{\theta}(\boldsymbol{X})] - \theta\}^2 + \mathrm{E}[\hat{\theta}^2(\boldsymbol{X})] - \{\mathrm{E}[\hat{\theta}(\boldsymbol{X})]\}^2 \\ \end{align} }$

2行目から3行目については、母数は定数なので期待値の外に出せます。ここで(3)式は分散の公式を使うと以下のように変形できますので、

$\displaystyle{ \begin{align} \mathrm{Var}(\hat{\theta}) = \mathrm{E}[\hat{\theta}^2(\boldsymbol{X})] - \{\mathrm{E}[\hat{\theta}(\boldsymbol{X})]\}^2 \end{align} }$

よって、(5)式は(1)式と(6)式を用いると、

$\displaystyle{ \begin{align} \mathrm{MSE}(\theta, \hat{\theta}) = \{ \mathrm{Bias}(\hat{\theta}(\boldsymbol{X}))\}^2 + \mathrm{Var}(\hat{\theta}(\boldsymbol{X})) \end{align} }$

となります。MSEはバイアスとバリアンスの足し算の式で書けることがわかります。なお、不偏推定量であればバイアスが0になりますので、MSEは分散と一致します。

参考

現代数理統計学の基礎

2020-07-03

ロジスティック回帰をはじめから

ロジスティック回帰について

ロジスティック回帰は機械学習の分類問題において、おそらく基本的な手法に位置するものであると思います。この手法は線形回帰モデルをその目的変数について、カテゴリカルな量に対応できるよう拡張させたものと考えることができます。カテゴリカルな量としてはたとえば、合格と不合格、生と死などそれらが連続的ではないものを指します。たとえば、合格と不合格かを判定するロジックを考えてみます。直感的には、それを説明する変数（たとえば勉強時間や難易度など）を考慮して合格する確率はxx%だというようにすると良さそうです。予測した確率に対してしきい値を設定して、たとえば50%とすればそれを境に合格と不合格を判定することができます。ということを考えると、目的変数としてはあるカテゴリになる確率を設定すれば良さそうです。つまり、ロジスティック回帰は確率の回帰を行うモデルです。しかしながら、線形回帰モデルで出力される値の範囲は-∞から∞です。一方で、確率は0から1です。そのため、線形結合した値について0から1の範囲に収まるように変換が必要になります。この変換にロジスティック分布を用いたものをロジスティック回帰モデル（ロジットモデル）と言います。

ロジスティック回帰モデルを導出する

予測の対象が合格か不合格かなど2値を取るような場合、確率モデルとしてベルヌーイ分布を考えることができます。たとえば合格する確率を $p$ とすれば、以下のように書けます。

$\displaystyle{ P(y = 合格) = p }$

一般化するために、合格を $1$ で表し、複数回分の結果があるとすれば、

$\displaystyle{ P(y_i = 1) = p_i }$

となり、それぞれが独立であるとすれば、

$\displaystyle{ y_i \sim Ber(p_i) }$

と表現できます。合格する確率が $p_i$ のベルヌーイ分布に従います。では、勉強時間などの要因と考えられる説明変数とその回帰係数の線形結合 $\boldsymbol{x}^\mathrm{T}_i\boldsymbol{\beta}$ の結果から確率 $p_i$ を見たいとします。つまり、線形結合 $\boldsymbol{x}^\mathrm{T}_i\boldsymbol{\beta}$ が増加するにつれて、この $p_i$ も増加するように線形結合の値を変換することを考えます。なお、線形結合の範囲は $[-\infty, \infty$ ]ですが確率の範囲は $[0, 1$ ]です。よって、変換するための関数の要件は以下になります。

線形結合が増加するにつれて、変換後の値も増加する（単調増加）
$[-\infty, \infty$ ]を $[0, 1$ ]にする

これには確率分布の累積分布関数（分布関数）が使えます。累積分布関数は単調増加関数であり、

$\displaystyle{ \lim_{x \to -\infty}F(x) = 0, \lim_{x \to +\infty}F(x) = 1 }$

という性質を持っています。この数式の意味するところは、 $x$ を $-\infty$ に持っていくと結果が $0$ に近づき、 $x$ を $\infty$ に持っていくと結果が $0$ に近づくということです。この性質は要件にピッタリです。この性質は要件にピッタリです。ただし、どんな確率分布の分布関数でも良いわけではなく、たとえば正規分布のように確率変数の範囲が $[-\infty, \infty$ ]である必要はあります。下図は標準正規分布の確率密度関数（左）と分布関数（右）です。

f:id:matakanobu:20200703172330p:plain — 標準正規分布の確率密度関数と分布関数

分布関数を $F$ として表現すると、

$\displaystyle{ p_i = F(\boldsymbol{x}^\mathrm{T}_i\boldsymbol{\beta}) \tag{1} }$

となります。線形結合を分布関数を用いて確率 $p_i$ に変換しているという意味です。この分布関数に、標準正規分布の分布関数である

$\displaystyle{ F(x) = \int^x_{-\infty}\frac{1}{\sqrt{2 \pi}} \exp{\left(-\frac{x^2}{2}\right)} \mathrm{d}x }$

を用いたモデルをプロビットモデルと言います。しかし、このモデルは分布関数に積分があるのでちょっと抵抗があります。分布関数に積分が無いものが扱いやすいですが、そのような確率分布はあるのでしょうか？それがロジスティック分布です。（標準）ロジスティック分布の分布関数はこちらです。

$\displaystyle{ F(x) = \frac{1}{1 + e^{-x}} }$

f:id:matakanobu:20200703174248p:plain — ロジスティック分布の確率密度関数と分布関数

こちらは標準正規分布と違って積分が無いので扱いやすいです。分母と分子に $e^x$ を掛けると、

$\displaystyle{ F(x) = \frac{e^{x}}{1 + e^{x}} \tag{2} }$

と変形ができ、(1)式と統合すると、

$\displaystyle{ p_i = \frac{\exp(\boldsymbol{x}^\mathrm{T}_i\boldsymbol{\beta})}{1 + \exp (\boldsymbol{x}^\mathrm{T}_i\boldsymbol{\beta})} \tag{3} }$

となります。これがロジスティック回帰モデル（ロジットモデル）です。

ロジスティック回帰モデルの数式を考察する

次にロジスティック回帰モデルの(3)式について、線形結合 $\boldsymbol{x}^\mathrm{T}_i\boldsymbol{\beta}=$ の形にしてみます。これはロジスティック分布の分布関数の逆関数で下式になります。

$\displaystyle{ \log \frac{p_i}{1-p_i} = \boldsymbol{x}^\mathrm{T}_i\boldsymbol{\beta} \tag{4} }$

左辺を対数オッズ（またはロジット）と言います。オッズは

$\displaystyle{ オッズ = \frac{p}{1-p} }$

で表され、ある事象が発生する確率を発生しない確率で割った値です。なので、

odds > 1 : 発生する確率のほうが大きい
odds < 1 : 発生しない確率のほうが大きい
odds = 1 : 発生する確率としない確率が等しい

です。対数を取る理由は発生する確率が大きくなる（ $p_i$ が大きくなる）と非常に大きな値になるから、そのスケールを抑えるためです。対数オッズを計算することはロジット変換と呼ぶこともあります。

$\displaystyle{ 対数オッズ：\mathrm{logit}(p) = \log \frac{p}{1-p} \tag{5} }$

これは確率（ $[0, 1$ ]）を $[-\infty, \infty$ ]の範囲のデータに変換することを意味します。ロジット変換は確率や比率に対する特徴量エンジニアリングとしても使うことができ、これによって線形回帰モデルを推定することも可能になります。

話をロジスティック回帰モデルに戻すと、(4)式からロジスティック回帰モデルは対数オッズの形で表されるモデルであるということがわかります。

ロジスティック回帰モデルの解釈

たとえば、勉強時間から合格する確率を説明（予測）するとして、ロジスティック回帰モデルを用いることにします。得られたデータから、勉強時間の回帰係数は以下のように推定されたものとします。

回帰係数:0.05
標準誤差:0.01

勉強時間の単位は1時間であるとします。このとき、勉強時間が増えた場合に合格確率はどう変化するかを考えてみます。線形回帰モデルのときはダイレクトに回帰係数の値を解釈すればよかったですが、ロジスティック回帰モデルでは変換が行われているため、そう単純でもなさそうです。そのために、まず(4)式を変形しますが、左辺の $\log$ を外すため両辺に指数を取ります。

$\displaystyle{ \frac{p_i}{1-p_i} = \exp(\boldsymbol{x}^\mathrm{T}_i\boldsymbol{\beta}) \tag{6} }$

これは上の話で出てきたオッズの形です。ベクトル $\boldsymbol{x}_i$ の1要因（変数）に着目するので、

$\displaystyle{ \frac{p_i}{1-p_i} = \exp(x_{i1}\beta_1 + \beta_0) }$

ここで $\beta_0$ は切片項です。 $x_{i1}$ が1単位増えた場合とそうでない場合を比較したいので、それらの比を取ります。

$\displaystyle{ \frac{\exp\{(x_{i1}+1)\beta_1 + \beta_0)\}}{\exp(x_{i1}\beta_1 + \beta_0)} = \exp(\beta_1) \tag{7} }$

これはオッズの比を取っているので、オッズ比と言います。さて、この式が意味するところは、 $x_{i1}$ が1単位増えたときのオッズは $\exp(\beta_1)$ 倍になるということです。仮に1単位増えたときに変化がないとした場合（これは回帰係数が0の場合）は分子と分母が同じになり、1になります。また、回帰係数が負の場合はオッズ比は1より小さくなります。これらの関係をまとめると以下のようになります。

オッズ比 > 1 : 可能性の増加（回帰係数が正）
オッズ比 < 1 : 可能性の減少（回帰係数が負）
オッズ比 = 1 : 変化なし（回帰係数が0）

実際に値を入れてみると理解できると思いますので、先の勉強時間と合格の例で計算してみます。

$\displaystyle{ \exp(0.05) - 1 = 0.051(5.1\%) }$

よって、1時間勉強時間を増やすと合格する確率は5.1%上昇するという解釈になります。標準誤差を用いれば95%信頼区間も構成できて、

$\displaystyle{ \exp(0.05 \pm 2 \times 0.01) = (1.030, 1.073) }$

となるので、およそ3〜7%の上昇であると見積もれます。

回帰係数の推定

ロジスティック回帰モデルの回帰係数は最尤推定を用いて推定するのが一般的です。 $y_i$ をラベルとして、 $y_i=1$ のときの確率を $p_i$ とすれば $y_i=0$ のときの確率は $1-p_i$ と書けます。尤度関数は確率の積になりますので、

$\displaystyle{ L(\boldsymbol{\beta}) = \prod^n_{i=1} p^{y_i}_i(1-p_i)^{1-y_i} }$

これはベルヌーイ分布の尤度関数です。なお、 $p_i$ はロジスティック回帰モデル(3)式であることに注意してください。これに対数を取った対数尤度関数は

$\displaystyle{ l(\boldsymbol{\beta}) = \sum^n_{i=1} \left\{ y_i \log p_i + (1 - y_i) \log (1-p_i) \right\} \tag{8} }$

となります。この対数尤度が最大となる $\boldsymbol{\beta}$ が解になりますので、 $\boldsymbol{\beta}$ で偏微分して0として解きます。しかしながら、これは解析的に解くことはできず数値的に求める必要があります。その手法としてニュートン・ラフソンアルゴリズムがありますが、今回の説明はここまでといたします。

参考

現代数理統計学の基礎
統計的学習の基礎

2020-06-23

リッジ回帰の考察（主成分分析との関係）

リッジ回帰の考察

こちらの回で求めた最小二乗法とリッジ回帰の解は以下のようになりました。

$\displaystyle{ \begin{eqnarray} 最小二乗回帰 : \hat{\boldsymbol{\beta}} &=& (\boldsymbol{X}^\mathrm{T} \boldsymbol{X})^{-1} \boldsymbol{X}^\mathrm{T} \boldsymbol{y} \\ \cr リッジ回帰 : \hat{\boldsymbol{\beta}} &=& (\boldsymbol{X}^\mathrm{T} \boldsymbol{X} + \lambda\boldsymbol{I})^{-1}\boldsymbol{X}^\mathrm{T} \boldsymbol{y} \tag{1} \end{eqnarray} }$

2つの違いは逆行列を求める際に、対角成分に $\lambda$ を足し合わせる点です。 $\lambda$ は罰則項の強さを調整するハイパーパラメータです。これが一体何を意味しているか式変形を行って解明していこうと思います。今回は特異値分解を使用しますので、参考に特異値分解の記事をリンクしておきます。

固有値分解と特異値分解の説明

まず、 $\boldsymbol{X}$ を以下のように特異値分解します。

$\displaystyle{ \boldsymbol{X} = \boldsymbol{U}\boldsymbol{D}\boldsymbol{V}^\mathrm{T} \tag{2} }$

ここで $\boldsymbol{U}$ と $\boldsymbol{V}$ は、 $\boldsymbol{U}^\mathrm{T}\boldsymbol{U}=\boldsymbol{I}, \space \boldsymbol{V}^\mathrm{T}\boldsymbol{V}=\boldsymbol{I}$ となるような行列です。 $\boldsymbol{D}$ は特異値を持つ対角行列です。(1)式の最小二乗回帰の解について、

$\displaystyle{ \hat{\boldsymbol{y}} = \boldsymbol{X}\hat{\boldsymbol{\beta}} = \boldsymbol{X}(\boldsymbol{X}^\mathrm{T}\boldsymbol{X})^{-1}\boldsymbol{X}^\mathrm{T}\boldsymbol{y} }$

これを特異値分解の式(2)を使って、

$\displaystyle{ \begin{eqnarray} \boldsymbol{X}\hat{\boldsymbol{\beta}} &=& \boldsymbol{X}(\boldsymbol{X}^\mathrm{T}\boldsymbol{X})^{-1}\boldsymbol{X}^\mathrm{T}\boldsymbol{y} \\ &=& \boldsymbol{U}\boldsymbol{D}\boldsymbol{V}^\mathrm{T}(\boldsymbol{V}\boldsymbol{D}\boldsymbol{U}^\mathrm{T}\boldsymbol{U}\boldsymbol{D}\boldsymbol{V}^\mathrm{T})\boldsymbol{V}\boldsymbol{D}\boldsymbol{U}^\mathrm{T}\boldsymbol{y} \\ &=& \boldsymbol{U}\boldsymbol{D}\boldsymbol{V}^\mathrm{T}\boldsymbol{V}\boldsymbol{D}^{-2}\boldsymbol{V}^\mathrm{T}\boldsymbol{V}\boldsymbol{D}\boldsymbol{U}^\mathrm{T}\boldsymbol{y} \\ &=& \boldsymbol{U}\boldsymbol{U}^\mathrm{T}\boldsymbol{y} \tag{3} \end{eqnarray} }$

一方リッジ回帰は、

$\displaystyle{ \begin{eqnarray} \boldsymbol{X}\hat{\boldsymbol{\beta}} &=& \boldsymbol{X}(\boldsymbol{X}^\mathrm{T}\boldsymbol{X} + \lambda\boldsymbol{I})^{-1}\boldsymbol{X}^\mathrm{T}\boldsymbol{y} \\ &=& \boldsymbol{U}\boldsymbol{D}\boldsymbol{V}^\mathrm{T}(\boldsymbol{V}\boldsymbol{D}\boldsymbol{U}^\mathrm{T}\boldsymbol{U}\boldsymbol{D}\boldsymbol{V}^\mathrm{T} + \lambda \boldsymbol{I})^{-1}\boldsymbol{V}\boldsymbol{D}\boldsymbol{U}^\mathrm{T}\boldsymbol{y} \\ &=& \boldsymbol{U}\boldsymbol{D}(\boldsymbol{D}^2 + \lambda \boldsymbol{I})^{-1}\boldsymbol{D}\boldsymbol{U}^\mathrm{T}\boldsymbol{y} \\ &=& \sum^p_{j=1} \boldsymbol{u}_j \frac{d^2_j}{d^2_j + \lambda} \boldsymbol{u}^\mathrm{T}_j \boldsymbol{y} \tag{4} \\ \end{eqnarray} }$

$\displaystyle{ (\boldsymbol{U} = [\boldsymbol{u}_1,\ldots, \boldsymbol{u}_p]) }$

なお、 $\lambda \ge 0$ なので、 $\frac{d^2_j}{d^2_j + \lambda} \le 1$ です。

最小二乗回帰とリッジ回帰との違いは $\frac{d^2_j}{d^2_j + \lambda} \le 1$ の部分です。この分だけ $\boldsymbol{U}^\mathrm{T}\boldsymbol{y}$ が縮小されます。ハイパーパラメータである $\lambda$ の値を大きくするか、 $d^2_j$ の値が小さくなるとより強く縮小されることを示しています。 $d^2_j$ は特異値の各成分の2乗値ですが、これは一体何を意味しているのでしょうか。

リッジ回帰と主成分分析の関係

実はリッジ回帰は主成分分析と関係があります。主成分分析は標本共分散行列を固有値分解して求めますが、標本共分散行列は下式で求めます。

$\displaystyle{ \boldsymbol{S} = \frac{\boldsymbol{X}^\mathrm{T}\boldsymbol{X}}{N} }$

ただし、 $\boldsymbol{X}$ は中心化（平均値で引かれている）されているとします。これを特異値分解の式(2)を使って変形すると、

$\displaystyle{ \begin{eqnarray} \frac{\boldsymbol{X}^\mathrm{T}\boldsymbol{X}}{N} &=& \frac{1}{N} \boldsymbol{V}\boldsymbol{D}\boldsymbol{U}^\mathrm{T}\boldsymbol{U}\boldsymbol{D}\boldsymbol{V}^\mathrm{T} \\ &=& \boldsymbol{V}\frac{\boldsymbol{D}^2}{N}\boldsymbol{V}^\mathrm{T} \end{eqnarray} }$

となります。 $\boldsymbol{V}$ は直交行列であり、 $\boldsymbol{D}$ は対角行列です。つまりこの式は（ $N$ を無視すれば） $\boldsymbol{X}^\mathrm{T}\boldsymbol{X}$ の固有値分解の形になっています。なので $\boldsymbol{V}$ の列ベクトル $\boldsymbol{v}_j$ は固有ベクトルであり、 $\boldsymbol{D}^2$ の各値は固有値です。各固有ベクトルは $\boldsymbol{X}$ の主成分方向を表しています。たとえば1列目の固有ベクトル $\boldsymbol{v}_1$ との内積

$\displaystyle{ \boldsymbol{z}_1 = \boldsymbol{X}\boldsymbol{v}_1 \tag{5} }$

を第1主成分と言います。第1主成分は他の主成分と比較して、最も大きい標本分散を持っています。つまり $\boldsymbol{v}_j$ は分散の方向を表しています。そして分散の大きさは固有値で表されています。固有値は $\boldsymbol{D}^2/N$ なので、第1主成分の分散は、

$\displaystyle{ \mathrm{Var}(\boldsymbol{z}_1) = \frac{d^2_1}{N} \tag{6} }$

となります。ここで(4)式の話に戻します。問題は $d^2_j$ が何であるのかということでした。これは(6)式が示しているように分散です。そしてリッジ回帰では $d^2_j$ が小さい変数についてより強く縮小するという話でした。 $d^2_j$ が小さいとはつまり、主成分方向にデータを射影したときの分散が小さくなることです。したがいまして、リッジ回帰は小さい分散方向の成分を持つ変数の係数をより強く縮小させることがわかります。

参考

現代数理統計学の基礎
統計的学習の基礎

2020-06-17

固有値分解と特異値分解の説明

固有値分解と特異値分解は統計学や機械学習において頻繁に登場する手法です。とくに、主成分分析やリッジ回帰をより深く理解する上で大いに役立ちます。

固有値分解(Eigenvalue Decompositions)

一言でいうと、与えられた行列を直交するベクトルとその大きさに分解することです。

では数式で説明してみます。固有値分解とは $\textbf{A}$ という次元が $n \times n$ の正方行列があったとき、下式のように分解することです。

$\displaystyle{ \textbf{A} = \textbf{P}\textbf{D}\textbf{P}^\mathrm{T} \tag{1} }$

ここで $\textbf{P}$ は直交行列、 $\textbf{D}$ は対角行列です。直交行列とは、

$\displaystyle{ \textbf{P}\textbf{P}^\mathrm{T} = \textbf{P}^\mathrm{T}\textbf{P} = \textbf{I} }$

となるような行列です。ここで $\textbf{I}$ は単位行列です。自身の内積 $\textbf{P}^\mathrm{T}\textbf{P}$ が単位行列になるとは $\textbf{P}$ の各ベクトル（基底）が直交しているということです。たとえば $\textbf{P}$ を

$\displaystyle{ \textbf{P} = \begin{bmatrix} p_{11} & p_{12} \cr \\ p_{21} & p_{22} \end{bmatrix} }$

というような行列としたときに、自身の内積は

$\displaystyle{ \textbf{P}^\mathrm{T}\textbf{P} = \begin{bmatrix} p^2_{11} + p^2_{21} & p_{11}p_{12} + p_{21}p_{22} \cr \\ p_{11}p_{12} + p_{21}p_{22} & p^2_{12} + p^2_{22} \end{bmatrix} }$

となり、これが単位行列になるので、

$\displaystyle{ \begin{bmatrix} p^2_{11} + p^2_{21} & p_{11}p_{12} + p_{21}p_{22} \cr \\ p_{11}p_{12} + p_{21}p_{22} & p^2_{12} + p^2_{22} \end{bmatrix} = \begin{bmatrix} 1 & 0 \cr \\ 0 & 1 \end{bmatrix} }$

これを各成分ごとに書くと、

$\displaystyle{ \begin{cases} p_{11}p_{12} + p_{21}p_{22} = 0 \cr \\ p^2_{11} + p^2_{21} = 1 \cr \\ p^2_{12} + p^2_{22} = 1 \tag{2} \end{cases} }$

となります。ところで、行列 $\textbf{P}$ をベクトルを使って表すと以下になります。

$\displaystyle{ \textbf{P} = \begin{bmatrix} \textbf{p}_{1} & \textbf{p}_{2} \tag{3} \end{bmatrix} }$

列成分をベクトルと置きました。さて(2)式と(3)式を見比べてみます。(2)式の各式は何を意味しているでしょうか。(2)式の1行目はベクトル[tex:\textbf{p}1, \textbf{p}2]の内積が0であると言っています。

$\displaystyle{ \textbf{p}_1 \cdot \textbf{p}_2 = p_{11}p_{12} + p_{21}p_{22} = 0 }$

内積が0であるとは2つのベクトルが直交していることを意味しています。これが直交行列という名前が付いている理由になります。また、(2)式の2行目と3行目はそれぞれベクトル[tex:\textbf{p}1, \textbf{p}2]の大きさが1であると言っています。つまり単位ベクトルです。これらのベクトルを固有ベクトルと言います。

直交行列を統計学や機械学習で見ると、列ベクトルは変数（特徴量）に該当しますので、各変数が無相関であるような行列です。

話を固有値分解に戻します。対角行列である $\textbf{D}$ は、対角成分以外が0である行列で、以下のような形状です。

$\displaystyle{ \begin{bmatrix} \lambda_1 & 0 & \ldots & 0 \\ 0 & \lambda_2 & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \ldots & \lambda_n \end{bmatrix} }$

この $\lambda_1, \ldots,\lambda_n$ を固有値と言います。固有ベクトルが方向を示していたのに対して、こちらはその大きさを表現しています。なぜ大きさであるのかは固有値分解の式である(1)を展開するとわかります。簡単のために $\textbf{A}$ を $2 \times 2$ の行列とします。

$\displaystyle{ \begin{aligned} \textbf{A} &= \textbf{P}\textbf{D}\textbf{P}^\mathrm{T} \cr \\ &= \begin{bmatrix} p_{11} & p_{12} \cr \\ p_{21} & p_{22} \end{bmatrix} \begin{bmatrix} \lambda_{1} & 0 \cr \\ 0 & \lambda_{2} \end{bmatrix} \begin{bmatrix} p_{11} & p_{21} \cr \\ p_{12} & p_{22} \end{bmatrix} \cr \\ &= \begin{bmatrix} \lambda_1 p_{11} & \lambda_2 p_{12} \cr \\ \lambda_1 p_{21} & \lambda_2 p_{22} \end{bmatrix} \begin{bmatrix} p_{11} & p_{21} \cr \\ p_{12} & p_{22} \end{bmatrix} \cr \\ &= \begin{bmatrix} \lambda_1 p^2_{11} + \lambda_2 p^2_{12} & \lambda_1 p_{11}p_{21} + \lambda_2 p_{12}p_{22} \cr \\ \lambda_1 p_{11}p_{21} + \lambda_2 p_{12}p_{22} & \lambda_1 p^2_{21} + \lambda_2 p^2_{22} \end{bmatrix} \cr \\ &= \lambda_1 \begin{bmatrix} p^2_{11} & p_{11}p_{21} \cr \\ p_{11}p_{21} & p^2_{21} \end{bmatrix} + \lambda_2 \begin{bmatrix} p^2_{12} & p_{12}p_{22} \cr \\ p_{12}p_{22} & p^2_{22} \end{bmatrix} \cr \\ &= \lambda_1 \textbf{p}_1\textbf{p}_1^\mathrm{T} + \lambda_2 \textbf{p}_2\textbf{p}_2^\mathrm{T} \end{aligned} }$

行列 $\textbf{A}$ は直交行列の列ベクトル（固有ベクトル）に対角行列の成分（固有値）をかけ合わせ、それらを固有ベクトルの分だけ足し合わせた形に分解できることを意味しています。固有ベクトルは大きさが1で方向のみを表しており、固有値はその方向の大きさを表していることがわかります。また、繰り返しになりますが、各固有ベクトルは直交しているため足し算の形に書き表すことができるわけです。

特異値分解(Singular Value Decompositions ; SVD)

固有値分解は正方行列であるという制限がありますが、一般的には特徴量の数よりデータの数のほうが多いと考えられます（行列が横より縦のほうが長い）。そこで特異値分解です。特異値分解は、ランクが $r$ 、次元が $n \times m$ の行列 $\textbf{A}$ を下式のように分解することです。

$\displaystyle{ \textbf{A} = \textbf{P}\textbf{\Lambda}\textbf{O}^\mathrm{T} \tag{4} }$

ここで各行列は以下のような性質をを持ったものです。

$\textbf{\Lambda}$ : 次元は $r \times r$ , 対角行列である
$\textbf{P}$ : 次元は $n \times r$ , $\textbf{P}^\mathrm{T}\textbf{P} = \textbf{I}$ となるような行列
$\textbf{O}$ : 次元は $m \times r$ , $\textbf{O}^\mathrm{T}\textbf{O} = \textbf{I}$ となるような行列

なお、 $\Lambda$ は $\lambda$ の大文字です。

$\textbf{\Lambda}$ は以下のように固有値の時と同じ形になっており、各値を特異値と言います。

$\displaystyle{ \begin{bmatrix} \lambda_1 & 0 & \ldots & 0 \\ 0 & \lambda_2 & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \ldots & \lambda_r \end{bmatrix} }$

固有値分解の時は $\textbf{P}$ が直交行列であるという性質を持っていましたが、特異値分解では正方行列から一般の行列に拡張を行ったためその性質は失われています。しかしながら、特異値分解対象である $\textbf{A}$ にフルランクであるという制約を設けることで一部その性質を得ることができます。フルランクであるとはランクと列数が同じ( $r = m$ )ということです。その制約により、(4)式自体には変化はありませんが、各行列は以下の性質のように変化します。

$\textbf{\Lambda}$ : 形状は $m \times m$ , 対角行列である
$\textbf{P}$ : 形状は $n \times m$ , $\textbf{P}^\mathrm{T}\textbf{P} = \textbf{I}$ となるような行列
$\textbf{O}$ : 形状は $m \times m$ , $\textbf{O}^\mathrm{T}\textbf{O} = \textbf{O}\textbf{O}^\mathrm{T} = \textbf{I}$ となるような直交行列

参考

現代数理統計学の基礎

2020-06-10

リッジ回帰の数式を導出する

最小二乗法による線形回帰モデル

リッジ回帰を説明する上で最小二乗法による線形回帰モデルの推定方法の理解が必要になりますので、まずこちらから説明します。

説明変数が2以上の線形回帰モデル（重回帰モデル）は以下で表されます。

$\displaystyle{ \boldsymbol{y} = \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{u} }$

ここで $\boldsymbol{y}$ は目的変数、 $\boldsymbol{X}$ は説明変数、 $\boldsymbol{\beta}$ は切片項を含む回帰係数（偏回帰係数）、 $\boldsymbol{u}$ は誤差項です。これの回帰係数 $\boldsymbol{\beta}$ を推定する方法として最小二乗法があります。これは観測値 $\boldsymbol{y}$ とモデルによる予測値 $\boldsymbol{\hat{y}}$ の差の二乗が $\boldsymbol{\beta}$ について最小になるような $\beta$ を見つける方法です。観測値と予測値の差の二乗を数式で書くと、

$\displaystyle{ h(\boldsymbol{\beta}) = \sum_{i=1}^{n} \left \{ y_i - (\beta_0 + \beta_1x_{i1} + \ldots + \beta_p x_{ip}) \right \}^2 }$

この式をベクトルと行列を用いて書き直すと、

$\displaystyle{ \begin{aligned} h(\boldsymbol{\beta}) &= \left \{\boldsymbol{y} - (\beta_0 + \beta_1 \boldsymbol{x}_1 + \ldots + \beta_p \boldsymbol{x}_p) \right \}^2 \\ &= \left \{\boldsymbol{y} - (\beta_0 + \beta_1 \boldsymbol{x}_1 + \ldots + \beta_p \boldsymbol{x}_p) \right \}^\mathrm{T} \left \{\boldsymbol{y} - (\beta_0 + \beta_1 \boldsymbol{x}_1 + \ldots + \beta_p \boldsymbol{x}_p) \right \} \\ &= (\boldsymbol{y} - \boldsymbol{X}\boldsymbol{\beta})^\mathrm{T}(\boldsymbol{y} - \boldsymbol{X}\boldsymbol{\beta}) \\ &= \|\boldsymbol{y} - \boldsymbol{X} \boldsymbol{\beta} \|^2 \end{aligned} }$

となります。なお、１行目から２行目の式変形はベクトルの2乗は自身の内積になるので $\boldsymbol{x}^2 = \boldsymbol{x} \cdot \boldsymbol{x} = \boldsymbol{x}^\mathrm{T} \boldsymbol{x}$ という関係を用いています。この観測値と予測値の差の二乗を残差平方和(RSS)と呼びます。ではこのRSSが $\boldsymbol{\beta}$ について最小となるようにするので、 $\boldsymbol{\beta}$ で偏微分して $\boldsymbol{0}$ と置きます。ベクトルで微分しているので、0ベクトルとなることに注意してください。

$\displaystyle{ \frac{\partial h(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = -2\boldsymbol{X}^\mathrm{T} (\boldsymbol{y} - \boldsymbol{X}\boldsymbol{\beta}) = \boldsymbol{0} }$

式変形をして $\boldsymbol{\beta}$ について解きます。

$\displaystyle{ \begin{aligned} \boldsymbol{X}^\mathrm{T}(\boldsymbol{y} - \boldsymbol{X} \boldsymbol{\beta}) &= \boldsymbol{0} \\ \boldsymbol{X}^\mathrm{T} \boldsymbol{y} - \boldsymbol{X}^\mathrm{T} \boldsymbol{X} \boldsymbol{\beta} &= \boldsymbol{0} \\ \boldsymbol{X}^\mathrm{T} \boldsymbol{X} \boldsymbol{\beta} &= \boldsymbol{X}^\mathrm{T} \boldsymbol{y} \\ \end{aligned} }$

$\beta$ イコールの形にするためには $\boldsymbol{X}^\mathrm{T} \boldsymbol{X}$ の逆行列を両辺の左からかける必要があります。逆行列はすべての行列において計算できるわけではなく、行列がフルランクである必要があります。 $\boldsymbol{X}$ がフルランクであるとは、言い換えると説明変数同士で完全相関（相関係数が1）している組が無いということです。なお、フルランクでない行列の場合、特異と呼びます。よって、ここではフルランクであることを仮定します。

$\displaystyle{ \begin{aligned} (\boldsymbol{X}^\mathrm{T} \boldsymbol{X})^{-1} \boldsymbol{X}^\mathrm{T} \boldsymbol{X} \boldsymbol{\beta} &= (\boldsymbol{X}^\mathrm{T} \boldsymbol{X})^{-1} \boldsymbol{X}^\mathrm{T} \boldsymbol{y} \\ \therefore \hat{\boldsymbol{\beta}} &= (\boldsymbol{X}^\mathrm{T} \boldsymbol{X})^{-1} \boldsymbol{X}^\mathrm{T} \boldsymbol{y} \\ \end{aligned} }$

これが最小二乗法を用いた場合の回帰係数（偏回帰係数）の推定量になります。

最小二乗法とバイアスについて

上の最小二乗法で求めた回帰係数は不偏推定量といって、バイアス（偏り）がありません。これは回帰係数の期待値（平均）が母数（真の値）に一致することを意味しています。しかしながら、回帰係数の選び方は無数にあるため、バイアスはあるもののもっと誤差を減らせる（精度を上げられる）ような回帰係数の選び方もあるかもしれません。そこで、バイアスを増やすことによってこれを実現することを考えます。変数選択もその一つであり、今回取り上げるリッジ回帰やラッソ回帰などの縮小推定があります。縮小推定では推定する回帰係数について制限を設けて縮小または0にする手法です。いずれにしても、縮小推定を用いて推定した回帰係数はバイアスを持つ可能性があるため、解釈には注意が必要になると考えられます。

# 最小二乗法による線形回帰をリッジ回帰に拡張する

最小二乗法による線形回帰は下式で表される残差平方和(RSS)を最小とするような $\boldsymbol{\beta}$ を求めることでした。

$\displaystyle{ h(\boldsymbol{\beta}) = \sum_{i=1}^{n} \left \{ y_i - (\beta_0 + \beta_1x_{i1} + \ldots + \beta_p x_{ip}) \right \}^2 }$

リッジ回帰を含む縮小推定では回帰係数 $\boldsymbol{\beta}$ に制限を制限を設けると上で書きました。縮小推定という名前が示しているように回帰係数 $\boldsymbol{\beta}$ が小さくなるようにします。残差平方和を最小とするというアプローチは変わりませんが、回帰係数 $\boldsymbol{\beta}$ を小さくしたいので回帰係数 $\boldsymbol{\beta}$ の大きさに応じてペナルティを課すようにします。具体的には、上の残差平方和の式に回帰係数 $\boldsymbol{\beta}$ の大きさを足し合わせます。単純に回帰係数 $\boldsymbol{\beta}$ を足し合わせないのは回帰係数 $\boldsymbol{\beta}$ はマイナスになることもあり得るからです。大きさという表現は数学的には様々な表現方法がありますが、回帰係数 $\boldsymbol{\beta}$ はベクトルであるためノルムになり、絶対値を取る方法(L1ノルム)、2乗する方法（L2ノルム）などがあります。前者を用いるとラッソ回帰に、後者だとリッジ回帰になります。

回帰係数 $\boldsymbol{\beta}$ の大きさを表す式を罰則項と呼ぶことにします。 $\beta_0$ は切片項に該当しますが、これは $\boldsymbol{y}$ の原点であり今回の関心から外れるものであり、罰則項には含めません。また、説明変数を中心化することで切片の推定量は単純に $\boldsymbol{\bar{y}}$ になります。これは切片の推定量 $\hat{\alpha}$ が

$\displaystyle{ \hat{\alpha} = \bar{y} - \hat{\beta} \bar{x} }$

から求めることができるので、説明変数が中心化されていることは説明変数の平均値が0になり、第2項が0になることからわかります。

以降では説明変数が中心化されていることを前提とします。つまり説明変数 $\boldsymbol{X}$ の次元が $n \times p$ （ $n$ 個のデータと $p$ 個の説明変数）です。罰則項を以下のように表現します。

$\displaystyle{ \lambda \sum_{j=1}^{p} \beta^2_j }$

ここで $\lambda$ は罰則項の強さを調整するハイパーパラメータです。ではこの項を残差平方和に足し合わせます。

$\displaystyle{ h(\boldsymbol{\beta}) = \sum_{i=1}^{n} \left \{ y_i - (\beta_1x_{i1} + \ldots + \beta_p x_{ip}) \right \}^2 + \lambda \sum_{j=1}^{p} \beta^2_j }$

最小二乗法のときと同様にベクトルと行列を用いて書き換えます。

$\displaystyle{ \begin{aligned} h(\boldsymbol{\beta}) &= (\boldsymbol{y} - \boldsymbol{X}\boldsymbol{\beta})^\mathrm{T}(\boldsymbol{y} - \boldsymbol{X}\boldsymbol{\beta}) +\lambda\boldsymbol{\beta}^\mathrm{T}\boldsymbol{\beta} \\ &= \|\boldsymbol{y} - \boldsymbol{X} \boldsymbol{\beta} \|^2 + \lambda \|\boldsymbol{\beta}\|^2 \end{aligned} }$

これを最小化するので偏微分して $\boldsymbol{0}$ とおきます。

$\displaystyle{ \frac{\partial h(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = -2 \boldsymbol{X}^\mathrm{T}(\boldsymbol{y} - \boldsymbol{X} \boldsymbol{\beta}) + 2\lambda \boldsymbol{\beta} = \boldsymbol{0} }$

$\boldsymbol{\beta}$ について解きます。

$\displaystyle{ \begin{aligned} -2 \boldsymbol{X}^\mathrm{T}(\boldsymbol{y} - \boldsymbol{X} \boldsymbol{\beta}) + 2\lambda \boldsymbol{\beta} &= \boldsymbol{0} \\ -\boldsymbol{X}^\mathrm{T} \boldsymbol{y} + \boldsymbol{X}^\mathrm{T} \boldsymbol{X} \boldsymbol{\beta} + \lambda \boldsymbol{\beta} &= \boldsymbol{0} \\ \boldsymbol{X}^\mathrm{T} \boldsymbol{X} \boldsymbol{\beta} + \lambda \boldsymbol{\beta} &= \boldsymbol{X}^\mathrm{T} \boldsymbol{y} \\ (\boldsymbol{X}^\mathrm{T} \boldsymbol{X} + \lambda\boldsymbol{I}) \boldsymbol{\beta} &= \boldsymbol{X}^\mathrm{T} \boldsymbol{y} \\ \end{aligned} }$

$\boldsymbol{I}$ は $p \times p$ の単位行列です。最小二乗法の時と同様に逆行列を用いて $\boldsymbol{\beta}$ イコールの形にしますが、今回のリッジ回帰では $\lambda\boldsymbol{I}$ つまり $\lambda$ が対角成分に足されています。これによって、必ずフルランクになりますので逆行列を持ちます。（個人的な感想ですが、これを始めて見たときにとても感動しました。）よって、

$\displaystyle{ \begin{aligned} (\boldsymbol{X}^\mathrm{T} \boldsymbol{X} + \lambda\boldsymbol{I})^{-1}(\boldsymbol{X}^\mathrm{T} \boldsymbol{X} + \lambda\boldsymbol{I}) \boldsymbol{\beta} &= (\boldsymbol{X}^\mathrm{T} \boldsymbol{X} + \lambda\boldsymbol{I})^{-1}\boldsymbol{X}^\mathrm{T} \boldsymbol{y} \\ \therefore \hat{\boldsymbol{\beta}} &= (\boldsymbol{X}^\mathrm{T} \boldsymbol{X} + \lambda\boldsymbol{I})^{-1}\boldsymbol{X}^\mathrm{T} \boldsymbol{y} \end{aligned} }$

これがリッジ回帰における回帰係数の推定量です。最小二乗法のものと比べてみます。

$\displaystyle{ \begin{aligned} 最小二乗法 : \hat{\boldsymbol{\beta}} &= (\boldsymbol{X}^\mathrm{T} \boldsymbol{X})^{-1} \boldsymbol{X}^\mathrm{T} \boldsymbol{y} \\ \cr リッジ回帰 : \hat{\boldsymbol{\beta}} &= (\boldsymbol{X}^\mathrm{T} \boldsymbol{X} + \lambda\boldsymbol{I})^{-1}\boldsymbol{X}^\mathrm{T} \boldsymbol{y} \end{aligned} }$

その違いは $\lambda\boldsymbol{I}$ が $\boldsymbol{X}^{\mathrm{T}}\boldsymbol{X}$ に足されていることだけです。ハイパーパラメータ $\lambda$ が $0$ であれば最小二乗法の式と等価に、 $\lambda \to \infty$ で回帰係数が $0$ に縮小されます。

この数式が一体何を意味するのかは次回考察します。

参考

現代数理統計学の基礎
統計的学習の基礎

2020-05-31

線形単回帰モデルを数式で理解する

線形単回帰モデル

線形単回帰(Simple linear regression)モデルは1つの入力（変数）を用いて予測をするものであり、回帰モデルの中でも最もシンプルなものです。シンプルであるがゆえに入力と出力との関係を解釈しやすいモデルです。また、特に訓練データの数が少ない場合においては、非線形モデルよりも良い性能を発揮することもあります。そのため、回帰モデルにおいて基礎となるモデルと考えられますので、その仕組みについて把握しておくこと、また使いこなせることは非常に重要であると考えております。

今回はこの線形単回帰モデルを数式を用いて解説を行っていこうと思います。

線形単回帰モデルの概観

入力を $x$ 、出力を $y$ としたときに、この2つの変数の間に直線で表せるような関係性があるとします。たとえば以下の図のような関係です。

f:id:matakanobu:20200528174617p:plain — データの散布図

これは直線の関係性があるといって良いと考えられますので、線形単回帰モデルを用いることで $x$ で $y$ を予測することが可能になります。ここで、直線を引くよりも、たとえば各点を結ぶジグサグの線を引くほうがより良い性能のモデルになるのではないのかという考えもあるかもしれません。しかしながら、それは間違いになります。機械学習においてはいわゆる過学習というものになり、それは交差検証を行うことで間違いであることがわかります。そして、間違いである理由は線形単回帰モデルはいわゆる統計モデルの一種になりますが、観測されるデータにはランダムの誤差 $\epsilon$ が含まれているからです。つまり統計モデルで考えた場合、以下のような関数からデータは生成されたと考えます。

$y = f(x) + \epsilon$

今回は、この $f(x)$ を線形単回帰モデルで推定することにします。これはつまり $y$ には誤差が足されているというわけですから、線形単回帰モデルによって推定された直線上にすべてのデータが乗ることはありえないわけで、直線付近にデータが散らばるのが適切なのです。この誤差 $\epsilon$ は予測な不可能な変数（確率変数）です。もう少し正確に言うと、モデルで予測をしたときに、予測しきれない部分を誤差としているということになります。さて、線形単回帰モデルは1変数の直線ですので、

$y = ax + b$

という中学校で習った1次関数を用います。 $a$ が傾きで $b$ が切片でした。ただ、線形単回帰モデルでは式が少しだけ異なっていて、こちらは誤差項を考慮しなければなりません。したがいまして、下式のようになります。

$y = ax + b + u$

また、線形単回帰モデルでは一般的に $a$ や $b$ ではなく、 $\alpha$ と $\beta$ を用います。くわえて、全部で $n$ 個のデータ（ $x_1$ のとき $y_1$ 、 $x_2$ のとき $y_2$ 、...、 $x_n$ のとき $y_n$ ）が観測されたとすると、線形単回帰モデルの数式は

$y_i = \alpha + \beta x_i + u_i$

と表されます。それぞれの文字には以下のように名前がついています。

$y_i$ : 応答変数(目的変数, 従属変数)(response variable)
$x_i$ : 説明変数(独立変数)(explanatory variable)
$\alpha$ : 切片項(intercept term)
$\beta$ : 回帰係数(regression coefficient)
$u_i$ : 誤差項(error term)

線形単回帰モデルの仮定

意外と知られていないかもしれないのですが、この線形単回帰モデルには以下の仮定が置かれています。

説明変数 $x_i$ は確率変数ではなく、定数である
誤差項 $u_i$ の期待値は0である( $E[u_i]=0$ )
誤差項同士は無相関である( $E[u_i, u_j]=0, i \neq j$ )（つまり自己相関は0である）
誤差項の分散は均一である( $Var(u_i)=\sigma^2$ )

1は目的変数の $y_i$ のみを確率変数として扱うということ、3は時系列性を持ってしまっているということです。これは線形単回帰モデルが時系列データを扱えないというわけはなく、誤差項に時系列性を持たせないように説明変数を選択すれば可能になります。これらの仮定から外れたデータに対して線形単回帰モデルを適用したときに、どのようなことが起こるのかは今度実験的に解説したいと思います。

回帰係数を求める方法（回帰係数の推定）

さて、得られたデータをこのモデルに当てはめるにはどうすればいいかを考えてみます。これはつまり、線形単回帰モデルのパラメータである $\alpha$ と $\beta$ を得られたデータから求めることになります。機械学習で言えば訓練（学習）です。誤差項の無い1次関数であれば2点（2データ）あれば、連列方程式を用いて切片と傾きを計算できました。しかし誤差項がある線形単回帰モデルではそう単純な話でもありません。まずは下の2つの直線を見比べてみてください。

これは3データに対して線形単回帰モデルを適用したものとお考えください。どちらの方がより適切と言えるでしょうか。直感的には右側の方が良さそうです。良さそうでは評価が曖昧ですので、その理由を定量的に評価できるように数式化してみましょう。

数式化するにあたって、まずなぜ直感的に右側の方が良さそうなのかを考えてみます。それは各観測データがモデルの直線に近いからであると思います。それを図示してみます。

緑色の丸が観測データ、青色の直線がモデル、赤色の矢印がモデルと観測データとの誤差を表しています。こう考えると、直線を引く上で（モデルを推定する上で）大事なことは誤差をなるべく小さくすることだとわかります。誤差は「観測値」と「モデルによる予測値」の差分なので数式で表すと、

$y_i - (\alpha + \beta x_i)$

ここで注意しなければならないのが、観測値が直線の上にあるか下にあるかで誤差の符号が変わってしまうので、絶対値を取るか2乗してプラスに統一する必要があります。ここでは、後の計算が楽になるので2乗することにします。すべての観測値の分だけ誤差を2乗し、それらを足し合わせます。

$\displaystyle{ h(\alpha, \beta) = \sum_{i=1}^{n} (y_i - \alpha - \beta x_i)^2 }$

これを最小化するような $\alpha$ と $\beta$ を求めれば良いことになります。このアプローチを最小二乗法と呼びます。 $h(\alpha, \beta)$ を $\alpha$ と $\beta$ について最小化するので、それぞれについて偏微分して0とおいた連立方程式を解けば求まります。

最小二乗法の計算

実際に $\alpha$ と $\beta$ を求めてみます。もっと効率が良い方法があるとは思うのですが、愚直に計算していってみます。

$\displaystyle{ h(\alpha, \beta) = \sum_{i=1}^{n} (y_i - \alpha - \beta x_i)^2 \tag{1} }$

この関数を $\alpha$ と $\beta$ について偏微分します。

$\displaystyle{ \begin{cases} \frac{\partial h(\alpha, \beta)}{\partial \alpha} &= -2 \sum (y_i - \alpha - \beta x_i) \\ \frac{\partial h(\alpha, \beta)}{\partial \beta} &= -2 \sum x_i(y_i - \alpha - \beta x_i) \end{cases} }$

イコール0と置きます。

$\displaystyle{ \begin{cases} -2 \sum (y_i - \alpha - \beta x_i) = 0 \\ -2 \sum x_i(y_i - \alpha - \beta x_i) = 0 \end{cases} }$

式を変形して $\alpha$ と $\beta$ を求めていきます。

$\displaystyle{ \begin{cases} \sum (y_i - \alpha - \beta x_i) = 0 \\ \sum x_i(y_i - \alpha - \beta x_i) = 0 \end{cases} }$

$\sum$ を外します（ $\sum^n_{I=1} \alpha = \alpha + \alpha +...+ \alpha = n \alpha; \alphaをn回足す$ ）。

$\displaystyle{ \begin{cases} \sum y_i - n \alpha - \beta \sum x_i) = 0 \\ \sum x_i y_i - \alpha \sum x_i - \beta \sum x_i^2 = 0 \end{cases} }$

$n$ はデータ数なのでゼロはありえないため( $n\ne0$ )、 $n$ で割ります。

$\displaystyle{ \begin{cases} \frac{1}{n} \sum y_i - \alpha - \beta \frac{1}{n} \sum x_i = 0 \\ \frac{1}{n} \sum x_i y_i - \alpha \frac{1}{n} \sum x_i - \beta \frac{1}{n} \sum x_i^2 = 0 \end{cases} }$

ここで $x_i$ と $y_i$ の平均値は以下のように計算ができるので、

$\displaystyle{ \frac{1}{n}\sum^n_{i=1} x_i = \bar{x},\ \frac{1}{n}\sum^n_{i=1} y_i = \bar{y} \tag{2} }$

$\bar{x}$ 、 $\bar{y}$ を用いると、

$\displaystyle{ \begin{cases} \bar{y} - \alpha - \beta \bar{x} = 0 \\ \alpha \bar{x} + \frac{1}{n} \beta \sum x_i^2 = \frac{1}{n} \sum x_i y_i \tag{3} \end{cases} }$

したがいまして、 $\alpha$ の推定量 $\hat{\alpha}$ は

$\displaystyle{ \hat{\alpha} = \bar{y} - \hat{\beta} \bar{x} \tag{4} }$

ここで $\hat{\beta}$ は $\beta$ の推定量です。では $\beta$ を求めていきます。(3)式に(4)式を代入します。

$\displaystyle{ \begin{aligned} \bar{x} \bar{y} - \beta \bar{x}^2 + \frac{1}{n} \beta \sum x_i^2 &= \frac{1}{n} \sum x_i y_i \\ \cr -n \bar{x}^2 \beta + \sum x_i^2 \beta &= \sum x_i y_i - n \bar{x} \bar{y} \\ \cr \left ( \sum x_i^2 - n \bar{x}^2 \right ) \beta &= \sum x_i y_i - n \bar{x} \bar{y} \\ \end{aligned} }$

$\sum$ でくくります。添字 $i$ が付いていない文字は $n$ が外に出ます。

$\displaystyle{ \begin{eqnarray} \sum \left (x_i^2 - \bar{x}^2 \right ) \beta &= \sum \left ( x_i y_i - \bar{x} \bar{y} \right ) \\ \cr \beta &= \frac{\sum \left ( x_i y_i - \bar{x} \bar{y} \right )}{\sum \left ( x_i^2 - \bar{x}^2 \right )} \tag{5} \end{eqnarray} }$

だいぶそれらしくなってきました。続いて分子と分母をそれぞれ変形します（天下り式ではありますが・・・）。まず分母です。

$\displaystyle{ \begin{aligned} \sum \left ( x_i^2 - \bar{x}^2 \right ) &= \sum \left ( x_i^2 - 2x_i \bar{x} + \bar{x}^2 + 2 x_i \bar{x} - 2 \bar{x}^2 \right ) \\ \cr &= \sum \left ( x_i^2 - 2x_i\bar{x}+\bar{x}^2 \right ) + 2\bar{x} \sum x_i - 2n \bar{x}^2 \\ \cr &= \sum \left ( x_i^2 - 2x_i\bar{x}+\bar{x}^2 \right ) + 2n \bar{x}^2 - 2n \bar{x}^2 \\ \cr &= \sum \left ( x_i - \bar{x} \right )^2 \end{aligned} }$

なお、2行目から3行目への変形は平均値の式である(2)式を用いています。次は分子です。

$\displaystyle{ \begin{aligned} \sum \left ( x_i y_i - \bar{x} \bar{y} \right ) &= \sum \left ( x_i y_i - \bar{y} x_i - \bar{x} y_i + \bar{x} \bar{y} + \bar{y}x_i + \bar{x} y_i -2 \bar{x} \bar{y} \right ) \\ \cr &= \sum \left ( x_i y_i - \bar{y} x_i - \bar{x} y_i + \bar{x} \bar{y} \right ) +n \bar{x} \bar{y} +n \bar{x} \bar{y} -2n \bar{x} \bar{y} \\ \cr &= \sum \left ( x_i - \bar{x} \right ) \left ( y_i - \bar{y} \right ) \end{aligned} }$

したがいまして、(5)式は、

$\displaystyle{ \beta = \frac{\sum \left ( x_i - \bar{x} \right ) \left ( y_i - \bar{y} \right )}{\sum \left ( x_i - \bar{x} \right )^2} \tag{6} }$

とてもキレイになりました！さらに、以下の記号を用意します。

$\displaystyle{ \begin{aligned} Q_{xx} &= \sum \left ( x_i - \bar{x} \right )^2 \\ \cr Q_{xy} &= \sum \left ( x_i - \bar{x} \right ) \left ( y_i - \bar{y} \right ) \end{aligned} }$

すると、 $\beta$ の推定量は

$\displaystyle{ \hat{\beta} = \frac{Q_{xy}}{Q_{xx}} \tag{7} }$

という形になります。 $\alpha$ の推定量と並べて書いてみます。

$\displaystyle{ \hat{\alpha} = \bar{y} - \hat{\beta} \bar{x}, \ \hat{\beta} = \frac{Q_{xy}}{Q_{xx}} }$

推定量の考察

求めた推定量の数式について少し考えてみます。まず $\hat{\alpha}$ です。

$\displaystyle{ \hat{\alpha} = \bar{y} - \hat{\beta} \bar{x} }$

$\beta$ の推定量が含まれているため切片の値は $\hat{\beta}$ に依存していますが、 $x$ の平均値がかけられています。説明変数である $x$ を中心化もしくは標準化すると $\bar{x}=0$ となるため切片は $\bar{y}$ つまり目的変数の平均値です。さらに $y$ を中心化もしくは標準化すると回帰の直線は原点 $(0, 0)$ を通ります。

$\displaystyle{ \hat{\beta} = \frac{\sum \left ( x_i - \bar{x} \right ) \left ( y_i - \bar{y} \right )}{\sum \left ( x_i - \bar{x} \right )^2} }$

つづいて $\hat{\beta}$ です。こちらは、分母が $x$ の分散、分子が $x$ と $y$ の共分散になっています。共分散は2変数の連動性を表しているもので、データが標準化されている場合は相関係数と一致します。よって $x$ と $y$ とで無関係（無相関）であれば分子が０になりますので、回帰係数つまり直線の傾きは0になります。

参考

現代数理統計学の基礎

アラフォーからの機械学習

アラフォーです。機械学習や統計学の解説記事を書いていきます

線形判別分析を理解する

線形判別分析

線形判別関数のパラメータ

対数比による分類

2次判別分析

線形判別分析による次元削減

参考

バイアスとバリアンスを理解する

参考

ロジスティック回帰をはじめから

ロジスティック回帰について

ロジスティック回帰モデルを導出する

ロジスティック回帰モデルの数式を考察する

ロジスティック回帰モデルの解釈

回帰係数の推定

参考

リッジ回帰の考察（主成分分析との関係）

リッジ回帰の考察

リッジ回帰と主成分分析の関係

参考

固有値分解と特異値分解の説明

固有値分解(Eigenvalue Decompositions)

特異値分解(Singular Value Decompositions ; SVD)

参考

リッジ回帰の数式を導出する

最小二乗法による線形回帰モデル

最小二乗法とバイアスについて

参考

線形単回帰モデルを数式で理解する

線形単回帰モデル

線形単回帰モデルの概観

線形単回帰モデルの仮定

回帰係数を求める方法（回帰係数の推定）

最小二乗法の計算

推定量の考察

参考