---
title: "Least Squares with Matrix 行列で最小二乗法"
author: "ryamada"
date: "2016年12月23日"
output:
html_document:
toc: true
toc_depth: 6
number_section: true
---
```{r setup, include=FALSE}
library(knitr)
opts_chunk$set(echo = TRUE)
library(rgl)
knit_hooks$set(rgl = hook_webgl)
```
１個の説明変数、１個の被説明変数があって、線形回帰を行うとき。
$$
Y \sim aX + b
$$
$$
\sum_{i=1}^n (y-\hat{y})^2 = \sum_{i=1}^n (y_i-(ax_i+b))^2
$$
を最小にするような、$a,b$の値を求める。
複数の説明変数があって、１個の被説明変数があるときには、説明変数の数だけの係数からなるベクトル$\mathbf{a}$を用いて
$$
y \sim X \mathbf{a} + b
$$
というモデルになる。
$$
\sum_{i=1}^n (y_i-\hat{y_i})^2 = \sum_{i=1}^n (y_i-(x_i^T a +b))^2 = ||y-(Xa+b)||^2
$$
を最小にするような$\mathbf{a},b$を求める。
今、$b=0$と固定したモデルの場合には
$$
y \sim Xa
$$
$$
||y-X\mathbf{a}||^2
$$
を最小にするような$\mathbf{a}$を求めることが課題である。
$b=0$でない場合には、$X$に１列足して、$X'$とする。
加えた列は、すべて値を1とする。
また、$\mathbf{a}$にも、１つ要素を増やしてやる。増やした分の推定値が$b$になる。
$$
a = (X^TX)^{-1}X^T y
$$
とにかく、計算してみる。
```{r}
d <- 1
n <- 1000
X <- matrix(rnorm(d*n),ncol=d)
a <- rnorm(d)
b <- rnorm(1)
a
y <- X %*% a + b + rnorm(n)*0.1
plot(X,y)
abline(b,a,col=2)
```
```{r}
X. <- cbind(X,rep(1,n))
a. <- solve(t(X.)%*%X.)%*%t(X.) %*% y
a.
print(c(a,b))
```
確かにうまく行っている。
$$
a = (X^T X)^{-1} X^T y
$$
どうしてこの式でよいのか、というのは、ふつうの統計学の教科書に書いてある。
それよりは、この式の行列のサイズを確認することに時間を使ってみよう。
+ $X$は$n\times d$行列。
+ $X^T$は$d \times n$ 行列。
+ $X^TX$ は$d \times d$正方行列。
+ $(X^TX)^{-1}$も$d \times d$ 正方行列。
+ $(X^TX)^{-1}X^T$は$d\times n$ 行列。
+ $(X^TX)^{-1}X^T y$ は$d \times 1$ 行列。
求めたい$\mathbf{a}$のそれと一致している。
```{r}
dim(X.)
dim(t(X.))
dim(t(X.)%*%X.)
dim(solve(t(X.)%*%X.))
dim(solve(t(X.)%*%X.)%*%t(X.))
dim(solve(t(X.)%*%X.)%*%t(X.)%*%y)
```
この要素の意味合いはなんだろうか？
+ $X$は$d$個の変数の値を列ベクトルとして持つ行列。
+ $X^T$は$n$サンプルごとの変数の値を列ベクトルとして持つ行列。
+ $X^TX$ は変数のペアワイズな内積を要素とする行列。異なる二つのベクトルの内積が大きいということは、それらが近い関係にあることを意味し、内積が0に近いということは、相互に独立性が強いことを意味する。
+ $(X^TX)^{-1}$は、ペアワイズな内積の逆数のような行列。逆にしたので、相互に独立な間柄を大きくとり扱い、相互に近い関係は軽く扱おうとする行列
+ $(X^TX)^{-1}X^T$は変数を重く扱うか軽く扱うかを考慮して、各サンプルの変数値を加減しなおした値を格納した行列
+ $(X^TX)^{-1}X^T y$ は変数の軽重を考慮した上で、それが$y$の値を大きくする方に働いているか、小さい方に働いているかで重みづけをして足し合わせた値。
+ それが係数。
もし、複数の説明変数があり、相互に相関が強いとすると、変数の重みづけの際に、両者の重みを小さ目にすることになり、結果として、個々の説明変数の係数は小さ目にすることになる。
以下のプロットを見ると、説明変数の値と、重みづけを勘案しなおした値とは線形な関係にあることが解る。
```{r}
plot(X.[,1],(solve(t(X.)%*%X.)%*%t(X.))[1,])
#plot(X.[,2],(solve(t(X.)%*%X.)%*%t(X.))[2,])
```
```{r}
library(mvtnorm)
d <- 3
n <- 1000
X <- rmvnorm(n,mean=rep(0,d),sigma=diag(rep(1,d)))
a <- rnorm(d)
b <- rnorm(1)
a
y <- X %*% a + b + rnorm(n)*0.1
X. <- cbind(X,rep(1,n))
a. <- solve(t(X.)%*%X.)%*%t(X.) %*% y
a.
print(c(a,b))
```
重み勘案前と後の値の関係をプロットしてみる。
これは３つの説明変数が相互に独立な場合である。
なぜなら、行列$X$を作ったときに３変数の分散共分散行列を単位行列で与えたからである。
```{r}
plot(X.[,1],(solve(t(X.)%*%X.)%*%t(X.))[1,])
plot(X.[,2],(solve(t(X.)%*%X.)%*%t(X.))[2,])
plot(X.[,3],(solve(t(X.)%*%X.)%*%t(X.))[3,])
```
説明変数３個の場合で、その３変数の分散共分散行列が単位行列でないように作成せよ。
正定値行列の作成の仕方を思い出すこと。
そのうえで、勘案前後のプロットをせよ。
各変数の間に相関が生じているので、勘案のされ方は、各サンプルの３変数の値の取り方によって差が生じるため、勘案前後の関係にばらつきが生じることを確認せよ。
```{r echo=FALSE}
d <- 3
library(GPArotation)
R <- Random.Start(d)
lambdas <- runif(d)
M <- R %*% diag(lambdas) %*% t(R)
X <- rmvnorm(n,mean=rep(0,d),sigma=M)
a <- rnorm(d)
b <- rnorm(1)
a
y <- X %*% a + b + rnorm(n)*0.1
X. <- cbind(X,rep(1,n))
a. <- solve(t(X.)%*%X.)%*%t(X.) %*% y
a.
print(c(a,b))
plot(X.[,1],(solve(t(X.)%*%X.)%*%t(X.))[1,])
plot(X.[,2],(solve(t(X.)%*%X.)%*%t(X.))[2,])
plot(X.[,3],(solve(t(X.)%*%X.)%*%t(X.))[3,])
```
実際に
$$
\mathbf{a} = (X^TX)^{-1}X^T
$$
をこのまま解こうとすると、$X^TX$の逆行列を計算する必要が出る。
逆行列の計算は負荷が大きく誤差が出やすいという性質もあるので、逆行列を算出せずに、$\mathbf{a}$を計算する方法がいくつか知られている。
いずれも、対角化・三角化行列を活用する方法である。
Rの線形回帰関数lm()では、デフォルトの計算方法としてQR法を取っていることが、lm()関数のソースを読むことで確認できる。
Rの関数は、関数名をプロンプトに書き込むことで、そのコードが読めることがある。lm()関数もそのようにしてコードが読める関数である。lm()関数のコードを表示させ、QR法がデフォルトであることを確認し、どのようにしてそれを確認したかを説明せよ。
```{r echo=TRUE}
lm
```
