線形回帰モデル　vol.3

# 前回記事
https://qiita.com/yudaiyamashita/items/73c5573fb174036d7c7d

# 中心極限定理

標本サイズが大きければ大きいほど分布が正規分布に近づくという定理。

# 正規分布
説明変数が次元2のベクトルの正規分布の例は下の通りです。高さ固定でどの角度から見ても、こんな感じです。
(えっ、バーバパパだって？)

![normalization.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/698527/e009e5e0-d1e7-d72d-0725-6a16c56164ce.png)


目的変数ベクトル(つまり回帰式$"y= f(x)"$の$y$の部分)$Y$は正規分布により観測されたと仮定する. 

こうすると誤差ベクトル$\boldsymbol \epsilon$はここまでの仮定に加え,　各成分は互いに独立で正規分布$N(0, \sigma^2)$に従う。具体的には以下の通りの式で同値となる. 

掛け算の省略記号$\Pi$(要は$\Sigma$記号の積バージョン)を用いると
確率密度関数$f(\epsilon)$は

$$f(\epsilon)= \prod_{i} \frac{1}{\sqrt{2\pi\sigma^2}} \exp (-\frac{1}{2\sigma^2} \epsilon_i^2)= \frac{1}{(2\pi\sigma^2)^{n/2}}\exp (-\frac{1}{2\sigma^2} \boldsymbol \epsilon^T\boldsymbol \epsilon)$$

そして誤差ベクトル$\epsilon$は平均ベクトル$\boldsymbol 0$, $\rm{cov} (\boldsymbol \epsilon)= \sigma^2 I_n$と表される。

参考

$\frac{1}{\sqrt{2\pi\sigma^2}} \exp (-\frac{1}{2\sigma^2} x^2)$を$-\infty$から$\infty$まで積分すると

$$\int_{-\infty}^{\infty} \frac{1}{\sqrt{2\pi\sigma^2}} \exp (-\frac{1}{2\sigma^2} x^2)dx= \frac{1}{\sqrt{2\pi\sigma^2}}\sqrt{\frac{\pi}{\alpha}} \hspace{1em}(\alpha= \frac{1}{2\sigma^2})= 1$$

(導出方法について簡単に)
$$\int_{-\infty}^{\infty}  \exp (-\beta x^2)dx$$
($\beta$は何らかの定数)
を2乗したものを計算した後、平方根をとることを考える。被積分関数は非負であることを考慮すると、正の平方根を考えればよい。

定義上、積分変数はそろえる必要もないので、
$$\int_{-\infty}^{\infty}  \exp (-\beta x^2)dx\int_{-\infty}^{\infty}  \exp (-\beta y^2)dy$$
無限小(余接バンドル)$dx, dy$は計算上微小量$\Delta x$などとして考えればよく(それこそ数学科でもなければ)、普通の数通り、交換法則、結合法則が成り立つから(・・・・本当はフビニの定理云々の話を挟む必要がある)
$$\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}  \exp[ (-\beta) (x^2+y^2)]dxdy$$
として置換積分を行う。多変数の置換積分はヤコビアン(ヤコビ行列)の知識が必要となる. 
結果だけ書けば
$x = r\cos\theta$, $y= r\sin \theta$, $dxdy=rdr d\theta$
となりあとは、簡単に計算できることがわかるだろう。

こういったものをガウス積分という

# ガウス型線形回帰モデル

誤差の出方に正規分布を仮定した線形回帰モデル

$\boldsymbol y= X\boldsymbol \beta+ \boldsymbol \epsilon$, $\boldsymbol \epsilon \sim N_n(\boldsymbol 0, \sigma^2 I_n)$

をガウス型線形回帰モデルという。
観測値ベクトル$\boldsymbol y$をガウス型線形回帰モデルとして
$E[\boldsymbol y]= E[X\boldsymbol \beta]+ E[\boldsymbol \epsilon]= X\boldsymbol \beta$, $\rm{cov}(\boldsymbol y)= E[(\boldsymbol y - \bar{\boldsymbol y})(\boldsymbol y - \bar{\boldsymbol y})^T] =  E[(\boldsymbol y - X\boldsymbol \beta)(\boldsymbol y - X\boldsymbol \beta)^T]= E[\boldsymbol \epsilon \boldsymbol \epsilon^T]= \sigma^2I$

($\because $ 今までの話どおり、仮定、平均は期待値、平均の平均は平均(元の値)、誤差平均は0(仮定)、など)

こうなると$N(X\boldsymbol \beta , \sigma^2I)$と記号では表示できる. 
尤度という用語を導入する。尤度(ゆうど)とは、確からしさのことで高ければ高いほど都合が良い。上限は下記の通り高々1。
変数$\boldsymbol \beta, \sigma^2$の関数とみた各尤度$f$の積、尤度関数$L$を次の様に定義する. ($L$はlikelihoodの頭文字。)
$$(\prod_i f(\epsilon_i, \sigma^2)=) L(\boldsymbol \beta, \sigma^2)= \frac{1}{(2\pi\sigma^2)^{n/2}}\exp (-\frac{1}{2\sigma^2} (\boldsymbol y - X\boldsymbol \beta)(\boldsymbol y - X\boldsymbol \beta)^T]$$

これを対数をとって, 対数尤度関数$l$とする。
$$l(\boldsymbol \beta, \sigma^2)=\log L(\boldsymbol \beta, \sigma^2)= -\frac{n}{2} \log (2\pi\sigma^2) - \frac{1}{2\sigma^2} (\boldsymbol y - X\boldsymbol \beta)(\boldsymbol y - X\boldsymbol \beta)^T$$
と表す。

これの最大値(最小値もだが)を求めれば元の最大値がわかるはず、
最大値を求めるのは2次関数でもなければ、微分が鉄則。

$$\frac{\partial l}{\partial \boldsymbol \beta} = \frac{1}{\sigma^2} (X^T \boldsymbol y -X^TX\boldsymbol \beta) = 0$$
$$\frac{\partial l}{\partial \sigma^2} = -\frac{n}{2\sigma^2} + -\frac{1}{2\sigma^4}(\boldsymbol y -X\boldsymbol \beta)^T(\boldsymbol y -X\boldsymbol \beta) = 0$$
この尤度方程式を解くことによって、$\boldsymbol \beta$と$\sigma^2$の最尤推定値
$$\hat{\boldsymbol \beta} =\frac{1}{\sigma^2} (X^T \boldsymbol y -X^TX\boldsymbol \beta) = 0$$

以下を$\beta$, $\sigma^2$の最尤推定値と呼ぶ. 
$\hat{\boldsymbol \beta}=(X^TX)^{-1}X^T \boldsymbol y$, $\hat{\sigma}^2= \frac{1}{n}(\boldsymbol y - X \hat{\beta})^T (\boldsymbol y - X \hat{\beta})$
最尤法は尤度をmaxにする方法である。
以上より最尤法によって推定した線形回帰モデルが求まる。
$$\boldsymbol y = \hat{\boldsymbol \beta}_0+ \hat{\boldsymbol \beta}_1x_1+ \cdots + \hat{\boldsymbol \beta}_p x_p= \hat{\boldsymbol \beta}\boldsymbol x$$


よって予測値ベクトル、残差ベクトルおよび誤差分散の推定値は以下の通りの式で与えられる.
これらは後述のvol.5(予定)での幾何的なイメージによってより明瞭に理解できるであろう。

$\hat{\boldsymbol y} = X \hat{\boldsymbol \beta}=(X^TX)^{-1}X^T \boldsymbol y$, $\boldsymbol e = \boldsymbol y - \hat{\boldsymbol y}= [I_n - X(X^TX)^{-1}X^T]\boldsymbol y$
$$\sigma^2 = \frac{1}{n}(\boldsymbol y - X \hat{\boldsymbol \beta})^T (\boldsymbol y - X \hat{\boldsymbol \beta}) \frac{1}{n}(\boldsymbol y -  \hat{\boldsymbol y})^T (\boldsymbol y -  \hat{\boldsymbol y})$$

以下一行、再掲

$$\boldsymbol y= X\boldsymbol \beta+ \boldsymbol \epsilon,　 \boldsymbol \epsilon \sim N_n(\boldsymbol 0, \sigma^2 I_n)$$

誤差ベクトル$\boldsymbol \epsilon$に対して$n$次元正規分布$N_n(\boldsymbol 0, \sigma^2 I_n)$を仮定したガウス型線形回帰モデルであった。これは, 観測値ベクトル$\boldsymbol y$が$n$次元正規分布$N_n(X\boldsymbol \beta, \sigma^2 I_n)$に従って観測されることを示す. (かっこの中身は平均と誤差の2乗平均である)

このとき,　回帰係数ベクトル$\boldsymbol \beta$の最尤推定量$\hat{\boldsymbol \beta}$ は, 正規分布に従う確率変数の線形結合であることから以下の論理より, 平均ベクトル$\boldsymbol \beta$, 平均共分散行列$\sigma^2(X^TX)^{-1}$の$(p+1)$次元正規分布$N_{p+1}(\boldsymbol \beta, \sigma^2(X^TX)^{-1})$に従うことがわかる. 

$n$次元確率ベクトル$Y$は, 平均ベクトル$\boldsymbol \mu$, 分散共分散行列$\Sigma$の
$n$次元正規分布$N_n(\boldsymbol \mu, \Sigma)$に従うものとすると、$m$次元確率ベクトル$Z=c+A\boldsymbol Y$は、平均ベクトル$c+A\boldsymbol \mu$, 分散共分散行列$A\Sigma A^T$の$m$次元正規分布$N_m(\boldsymbol c + A \boldsymbol \mu, A\Sigma A^T)$に従う.

よって$A=(X^TX)^{-1}X^T$とすると　$\boldsymbol \beta$の最尤推定量$\hat{\boldsymbol \beta}= (X^TX)^{-1}X^T \boldsymbol y$は$N_{p+1}(\boldsymbol \beta, \sigma^2(X^TX)^{-1})$であることが示された.

https://qiita.com/yudaiyamashita/items/b5d09edb377a232fbba6

続く