最小二乗法

$x_i$は$p-1$個の要素を持つ説明変数ベクトルとする。
$y_i$は目的変数とすると線形モデルは下記の様に表される

```math
y_i=\beta_1+\sum_{j=2}^{p}x_{ij}\beta_j
```

この時、$\beta$を切片項を含むp要素のベクトルとし
$X_1$ ~ $X_N$を切片項の1を含むN行p列の行列としてまとめると

```math
\hat{y}=X^T\hat{\beta}
```

とかける。この時、

```math
RSS(\beta)=\sum_{i=1}^N(y_i -x_i^T\beta)^2
```

を最小とする$\beta$を求めたい。($X\beta)^Tyとy^T(X\beta)$はベクトルであることからかける順番は気にしなくてよく

```math
\begin{align}
RSS(\beta)&=(y-X\beta)^T(y-X\beta)\\
&=y^Ty-(X\beta)^Ty-y^TX\beta+(X\beta)^TX\beta\\
&=y^Ty-2(X\beta)^Ty+(X\beta)^TX\beta\\
&=y^Ty-2(X\beta)^Ty+\beta^TX^TX\beta
\end{align}
```

この値を最小化する$\beta$を求めたい。$f(\beta)$とみて$\beta$について偏微分し0となる$\beta$を求める。まずは第二項の$P(\beta)=2(X\beta)^Ty$について考える


```math
\qquad P(\beta) = 2
\begin{bmatrix}

\begin{pmatrix}
  x_{11} & x_{12} & \cdots & x_{1p}\\
  x_{21} & x_{22} & \cdots & x_{2p}\\
  \cdots & \cdots & \cdots & \cdots\\
  x_{N1} & x_{N2} & \cdots & x_{Np}
\end{pmatrix}

\begin{pmatrix}
  \beta_{1}\\
  \beta_{2}\\
  \cdots\\
  \beta_{p}\\
\end{pmatrix}

\end{bmatrix}^T

  \begin{pmatrix}
    y_{1}\\
    y_{2}\\
    \cdots\\
    y_{N}\\
  \end{pmatrix}\\

  P(\beta)=2

  \begin{bmatrix}

  \begin{pmatrix}
    x_{11}\beta_{1} + \cdots + x_{1p}\beta_{p}\\
    x_{21}\beta_{1} + \cdots + x_{2p}\beta_{p}\\
    \cdots + \cdots + \cdots\\
    x_{N1}\beta_{1} + \cdots + x_{Np}\beta_{p}\\
  \end{pmatrix}

  \end{bmatrix}^T

  \begin{pmatrix}
    y_{1}\\
    y_{2}\\
    \cdots\\
    y_{N}\\
  \end{pmatrix}

```

```math
\begin{align}

  P(\beta) &= 2(x_{11}\beta_{1}+\cdots+x_{1p}\beta_{p})y_1\\
   &+ 2(x_{21}\beta_{1}+\cdots+x_{2p}\beta_{p})y_2\\
   &+\cdots\\
   &+2 (x_{N1}\beta_{1}+\cdots+x_{Np}\beta_{p})y_N\\

  P(\beta) &= \sum_{j=1}^Ny_j\sum_{i=1}^px_{ji}\beta_i

\end{align}
```

ここで$\beta_1$から$\beta_p$についてそれぞれ微分すると

```math
\frac{\delta P}{\delta \beta_1} = 2(x_{11}y_1+x_{21}y_2+\cdots+x_{N1}y_N) \\
\frac{\delta P}{\delta \beta_2} = 2(x_{12}y_1+x_{22}y_2+\cdots+x_{N2}y_N) \\
\cdots \\
\frac{\delta P}{\delta \beta_p} = 2(x_{1p}y_1+x_{2p}y_2+\cdots+x_{Np}y_N)
```

ここで、XはN行p列、YはN要素ベクトルであり、$\frac{\delta P}{\delta \beta}$をp要素のベクトルと考えると

```math
\frac{\delta P}{\delta \beta}=2X^Ty\\
```

次に$Q(\beta)=\beta^TX^TX\beta$を\betaについて微分する

```math
\begin{align}

  Q(\beta) &=

  \begin{pmatrix}
    \beta_{1}\\
    \beta_{2}\\
    \cdots\\
    \beta_{p}\\
  \end{pmatrix}^T

  \begin{pmatrix}
    x_{11} & x_{12} & \cdots & x_{1p}\\
    x_{21} & x_{22} & \cdots & x_{2p}\\
    \cdots & \cdots & \cdots & \cdots\\
    x_{N1} & x_{N2} & \cdots & x_{Np}
  \end{pmatrix}^T

  \begin{pmatrix}
    x_{11} & x_{12} & \cdots & x_{1p}\\
    x_{21} & x_{22} & \cdots & x_{2p}\\
    \cdots & \cdots & \cdots & \cdots\\
    x_{N1} & x_{N2} & \cdots & x_{Np}
  \end{pmatrix}

  \begin{pmatrix}
    \beta_{1}\\
    \beta_{2}\\
    \cdots\\
    \beta_{p}\\
  \end{pmatrix}\\

&=

  \begin{pmatrix}
    \beta_{1} & \cdots & \beta_{p}\\
  \end{pmatrix}

  \begin{pmatrix}
    x_{11} & x_{21} & \cdots & x_{N1}\\
    x_{12} & x_{22} & \cdots & x_{N2}\\
    \cdots & \cdots & \cdots & \cdots\\
    x_{1p} & x_{2p} & \cdots & x_{Np}
  \end{pmatrix}

  \begin{pmatrix}
    x_{11} & x_{12} & \cdots & x_{1p}\\
    x_{21} & x_{22} & \cdots & x_{2p}\\
    \cdots & \cdots & \cdots & \cdots\\
    x_{N1} & x_{N2} & \cdots & x_{Np}
  \end{pmatrix}

  \begin{pmatrix}
    \beta_{1}\\
    \beta_{2}\\
    \cdots\\
    \beta_{p}\\
  \end{pmatrix}

\end{align}

```

$X^TX$について考えると

```math
\begin{align}
X^TX
&=
  \begin{pmatrix}
    x_{11} & x_{21} & \cdots & x_{N1}\\
    x_{12} & x_{22} & \cdots & x_{N2}\\
    \cdots & \cdots & \cdots & \cdots\\
    x_{1p} & x_{2p} & \cdots & x_{Np}
  \end{pmatrix}

  \begin{pmatrix}
    x_{11} & x_{12} & \cdots & x_{1p}\\
    x_{21} & x_{22} & \cdots & x_{2p}\\
    \cdots & \cdots & \cdots & \cdots\\
    x_{N1} & x_{N2} & \cdots & x_{Np}
  \end{pmatrix}\\
&=
  \begin{pmatrix}

    x_{11}x_{11} + x_{21}x_{21} + \cdots + x_{N1}x_{N1}
      & x_{11}x_{12} + x_{21}x_{22} + \cdots + x_{N1}x_{N2}
      & \cdots
      & x_{11}x_{1p} + x_{21}x_{2p} + \cdots + x_{N1}x_{Np}\\

    x_{12}x_{11} + x_{22}x_{21} + \cdots + x_{N2}x_{N1}
      & \cdots
      & \cdots
      & \cdots\\

    \cdots
      & \cdots
      & \cdots
      & \cdots\\

    x_{1p}x_{11} + x_{2p}x_{21} + \cdots + x_{Np}x_{N1}
      & \cdots
      & \cdots
      & x_{1p}x_{1p} + x_{2p}x_{2p} + \cdots + x_{Np}x_{Np}\\
  \end{pmatrix}\\
&=
  \begin{pmatrix}
    \sum_i^N x_{i1}x_{i1} & \sum_i^N x_{i1}x_{i2} & \cdots & \sum_i^N x_{i1}x_{ip}\\
    \sum_i^N x_{i2}x_{i1} & \sum_i^N x_{i2}x_{i2} & \cdots & \sum_i^N x_{i2}x_{ip}\\
    \cdots & \cdots & \cdots & \cdots & \\ 
    \sum_i^N x_{ip}x_{i1} & \sum_i^N x_{ip}x_{i2} & \cdots & \sum_i^N x_{ip}x_{ip}\\
  \end{pmatrix}
\end{align}
```

各項は行番号をr、列番号をcとして、$\sum_i^N x_{ir}x_{ic}$と一般化できる。これを$\sum_i^N x_{ir}x_{ic} = X^2_{rc}$と表すとして

```math
X^TX \beta
=
\begin{pmatrix}
  X^2_{11}\beta_1 + X^2_{12}\beta_2 + \cdots + X^2_{1p}\beta_p\\
  X^2_{21}\beta_1 + X^2_{22}\beta_2 + \cdots + X^2_{2p}\beta_p\\
  \cdots\\
  X^2_{p1}\beta_1 + X^2_{p2}\beta_2 + \cdots + X^2_{pp}\beta_p
\end{pmatrix}
```

```{math}
\begin{align}
Q(\beta)=
\beta^T X^TX \beta
=& (X^2_{11}\beta_1 + X^2_{12}\beta_2 + \cdots + X^2_{1p}\beta_p)\beta_1\\
&+ (X^2_{21}\beta_1 + X^2_{22}\beta_2 + \cdots + X^2_{2p}\beta_p)\beta_2\\
&\cdots\\
&+ (X^2_{p1}\beta_1 + X^2_{p2}\beta_2 + \cdots + X^2_{pp}\beta_p)\beta_p\\

\frac{\delta Q}{\delta \beta_1}
=& 2X^2_{11}\beta_1 + X^2_{12}\beta_2 + X^2_{13}\beta_3 +\cdots + X^2_{1p}\beta_p\\
&+ X^2_{21}\beta_2 + X^2_{31}\beta_3 + \cdots + X^2_{p1}\beta_p
\end{align}
```

$X^2_{rc}=X^2_{cr}$であることから

```math
\frac{\delta Q}{\delta \beta_1}
= 2(X^2_{11}\beta_1 + X^2_{12}\beta_2 + X^2_{13}\beta_3 +\cdots + X^2_{1p}\beta_p)\\

\frac{\delta Q}{\delta \beta_2}
= 2(X^2_{21}\beta_1 + X^2_{22}\beta_2 + X^2_{23}\beta_3 +\cdots + X^2_{2p}\beta_p)\\

\frac{\delta Q}{\delta \beta_p}
= 2(X^2_{p1}\beta_1 + X^2_{p2}\beta_2 + X^2_{p3}\beta_3 +\cdots + X^2_{pp}\beta_p)\\
```

これは$X^2\beta=X^TX\beta$である

つまり

```math
\frac{\delta RSS(\beta)}{\delta \beta} = 2X^TX\beta-2X^Ty=0\\
X^TX\beta = X^Ty
```

となればよく、$X^TX$が正則であるならば、両辺に左からから$(X^TX)^{-1}$をかけて

```math
\beta=(X^TX)^{-1}X^Ty
```

Rで実装すると

```{R}
y <- c(80,90,120,110,130,170) #output vector

x0 <- c(1,1,1,1,1,1) #intersept
x1 <- c(5,5,7,5,8,2)
x2 <- c(6,8,10,12,12,12)
X <- cbind(x0, x1, x2) #input matrix

b <- solve(t(X) %*% X) %*% t(X) %*% y

b

lm(y~x1+x2)
```