［R］ゼロから作る最小二乗法１：単回帰

「自分で数式をコードに落としていって動かす」という作業は非常に勉強になると思ったので，いろんなアルゴリズムをゼロから作っていきたいと思います。第一弾として単回帰からはじめていきます。

理論の要点を整理してから実装する構成で述べていきます。

最小二乗法だけでなく統計的推論についても触れていきます。

モデル

線形回帰（linear regression）は，予測の目的変数\(Y_i\)と説明変数\(X_{i1}, X_{i2}, ..., X_{ik}\)の間に次のような線形関係を仮定したモデルを置いて予測する手法です。

\[ Y_i = \beta_0 + \beta_1X_{i1} + \beta_2 X_{i2} + \cdots + \beta_k X_{ik} + u_i \]

今回は線形回帰のなかでも説明変数が１つの単回帰モデル（simple regression model） \[ Y_i = \beta_0 + \beta_1X_{i} +u_i \] での場合で考えます。

誤差関数

パラメータ\(\beta_0, \beta_1\)を推定するために多く使われる方法は，実測値\(Y_i\)と予測値\(\hat{Y}_i\)の誤差2乗和（Sum of Squared Error）

\[ SSE = \sum_{i=1}^n (Y_i - \hat{Y}_i)^2 \]

あるいは平均2乗誤差（Mean of Squared Error: MSE）

\[ MSE = \frac{1}{n} \sum_{i=1}^n (Y_i - \hat{Y}_i)^2 \]

を最小にするパラメータを採用するという最小二乗法（least squares method）です。

どちらの誤差関数も同様にU字型（下に凸）の関数になり，関数の接線の傾きが0になる地点（関数を微分して0になる点）で誤差が最小になるパラメータが得られます¹。

どちらの誤差関数でも，最小二乗法でのパラメータ推定式は同じものになります。

パラメータの推定

回帰係数のパラメータは「誤差関数（例えばSSE）を微分してゼロになる点」という次の条件 \[ \begin{cases}\frac{\partial SSE}{\partial \beta_0} = -2 \sum_{i=1}^n (Y_i - \beta_0 - \beta_1 X_i) = 0\\\frac{\partial SSE}{\partial \beta_1} = -2 \sum_{i=1}^n (Y_i - \beta_0 - \beta_1 X_i)X_i = 0\end{cases} \] の解となるものであり， \[ \begin{align} \hat{\beta}_0 &= \bar{Y} - \beta_1 \bar{X}\\ \hat{\beta}_1 &= \frac{\sum_{i=1}^n(X_i-\bar{X})(Y_i-\bar{Y_i})} {\sum_{i=1}^n (X_i - \bar{X})^2} \end{align} \] になります（\(\bar{Y},\bar{X}\)はそれぞれ\(Y_i,X_i\)の平均値）。

Rで実装

式に入れるだけなのでシンプルですね

# 関数を定義
SimpleRegression_OLS <- function(x, y) {
  # β_1の推定
  beta_1 = sum((x - mean(x))*(y - mean(y))) / sum((x - mean(x))^2)
  # β_0の推定
  beta_0 = mean(y) - beta_1 * mean(x)
  
  # 結果を格納
  result = c("beta_0" = beta_0, "beta_1" = beta_1)
  return(result)
}

# womenデータセットで試す
SimpleRegression_OLS(x = women$height, y = women$weight)

##    beta_0    beta_1 
## -87.51667   3.45000

# lm()と比較
res <- lm(weight ~ height, women)
res

## 
## Call:
## lm(formula = weight ~ height, data = women)
## 
## Coefficients:
## (Intercept)       height  
##      -87.52         3.45

同じ値が出ました。

統計的検定

$u_i \sim N(0, \sigma^2)$：誤差項$u_i$は相互に独立に平均0，分散$\sigma^2$の正規分布$N(0, \sigma^2)$に従う

という仮定を置いた単回帰モデルをsimple normal regression modelと呼び，そこでは次のような方法で統計的推論を行うことができます。

標準誤差

推定量\(\hat{\beta}_0, \hat{\beta}_1\)の標準誤差は $$SE(\hat{\beta}_0)=\sqrt{\left(\frac{1}{n}+\frac{\bar{X}^2}{\sum_{i=1}^n(X_i - \bar{X})^2}\right)\hat{\sigma}^2} \\ SE(\hat{\beta}_1)=\sqrt{\frac{\hat{\sigma}^2}{\sum_{i=1}^n(X_i-\bar{X})^2}}$$ で，ここで$\hat{\sigma}^2$は誤差項の分散の推定量（残差分散）であり，残差平方和を\(n-2\)で割ったもの $$\hat{\sigma}^2 = \frac{\sum_{i=1}^n \hat{u}_i^2}{n-2}$$ です。

検定統計量

検定統計量は次のようになります²。 \[ t = \frac{\hat{\beta}_1 - \beta_1}{SE(\hat{\beta}_1)} \]

これは自由度\(n-2\)の\(t\)分布に従うことが知られています。

仮説体系は多くの場合 \[ H_0: \beta_1 = 0, \hspace{2em} H_1: \beta_1 \neq 0 \] なので，実際に使う検定統計量は \[ t = \frac{\hat{\beta}_1}{SE(\hat{\beta}_1)} \] になります。

そして，有意水準\(\alpha\)のもとで次の判定を行います。 $$\begin{cases} |t| > t_{\alpha / 2}(n-2)\text{のとき，帰無仮説を棄却する}\\ |t| \leq t_{\alpha / 2}(n-2)\text{のとき，帰無仮説を棄却しない} \end{cases}$$

あるいは\(p\)値の算出をします。

Rで実装

まず標準誤差を出します

# 標準誤差
StandardError <- function(x, y){
  # 残差の取得
  coef = SimpleRegression_OLS(x,y)
  y_hat = coef[1] + coef[2] * x
  u = y - y_hat
  
  # 誤差項の推定
  n = length(x)
  sigma_hat = sum(u^2)/(n-2)
  
  # 標準誤差(β_0)の推定
  A = (1/n)+(mean(x)^2/sum((x - mean(x))^2))
  SE_0 = sqrt(A * sigma_hat)
  
  # 標準誤差(β_1)の推定
  SE_1 = sqrt(sigma_hat / sum((x - mean(x))^2))
  
  # 結果の格納
  result = c("SE_0" = SE_0, "SE_1" = SE_1)
  return(result)
}

# womenデータセットで試す
StandardError(x = women$height, y = women$weight)

次に検定を処理します

tTest <- function(x, y){
  # 係数を算出
  beta = SimpleRegression_OLS(x, y)
  # 標準誤差を算出
  SE = StandardError(x, y)
  # 検定統計量を算出
  t_0 = beta[1] / SE[1]
  t_1 = beta[2] / SE[2]
  # p値
  n = length(x)
  p_0 = pt(q = -abs(t_0), df = (n-2))*2 # 両側検定なので2倍する
  p_1 = pt(q = -abs(t_1), df = (n-2))*2 # 両側検定なので2倍する
  # 結果
  result = matrix(c(t_0, t_1, p_0, p_1), ncol = 2,
                  dimnames = list(c("beta_0","beta_1"),c("t_value","p_value")))
  return(result)
}

# womenデータセットで試す
tTest(x = women$height, y = women$weight)

##          t_value      p_value
## beta_0 -14.74103 1.711082e-09
## beta_1  37.85531 1.090973e-14

lmのほうと比較してみます

summary(res)

## 
## Call:
## lm(formula = weight ~ height, data = women)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7333 -1.1333 -0.3833  0.7417  3.1167 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -87.51667    5.93694  -14.74 1.71e-09 ***
## height        3.45000    0.09114   37.85 1.09e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.525 on 13 degrees of freedom
## Multiple R-squared:  0.991,  Adjusted R-squared:  0.9903 
## F-statistic:  1433 on 1 and 13 DF,  p-value: 1.091e-14

ちゃんと同じ値になっています。

おまけ

OLS推定量の導出

SSEを誤差関数にした最小二乗推定量の導出についてメモ。

まず， \[ SSE = \sum_{i=1}^n (Y_i - \hat{Y}_i)^2 = \sum_{i=1}^n (Y_i - \beta_0 - \beta_1 X_i)^2 \] の微分は，合成微分律から， \[ \frac{\partial SSE}{\partial \beta_0} = -2 \sum_{i=1}^n (Y_i - \beta_0 - \beta_1 X_i) \\ \frac{\partial SSE}{\partial \beta_1} = -2 \sum_{i=1}^n (Y_i - \beta_0 - \beta_1 X_i)X_i \] となる。

誤差を最小化するパラメータは次の条件の解である \[ \begin{cases} \frac{\partial SSE}{\partial \beta_0} = -2 \sum_{i=1}^n (Y_i - \beta_0 - \beta_1 X_i) = 0 \\ \frac{\partial SSE}{\partial \beta_1} = -2 \sum_{i=1}^n (Y_i - \beta_0 - \beta_1 X_i)X_i = 0 \end{cases} \] これを整理すると正規方程式（normal equations） \[ \begin{cases} \begin{align} n \beta_0 + \sum_{i=1}^n \beta_1 X_i &= \sum_{i=1}^n Y_i \\ \sum_{i=1}^n \beta_0 X_i + \sum_{i=1}^n \beta_1 X_i^2 &= \sum_{i=1}^n X_iY_i \end{align} \end{cases} \] が得られ，このうち上の式に\((1/n)\sum_{i=1}^n X_i\)を乗じて \[ \begin{cases} \begin{align} \sum_{i=1}^n \beta_0 X_i + \frac{1}{n}\beta_1 (\sum_{i=1}^n X_i)^2 &= \frac{1}{n}(\sum_{i=1}^n X_i)(\sum_{i=1}^n Y_i) \\ \sum_{i=1}^n \beta_0 X_i + \sum_{i=1}^n \beta_1 X_i^2 &= \sum_{i=1}^n X_iY_i \end{align} \end{cases} \] 下の式から引けば \[ \begin{align} \sum_{i=1}^n \beta_1 X_i^2 - \frac{1}{n}\beta_1 (\sum_{i=1}^n X_i)^2 &= \sum_{i=1}^n X_iY_i - \frac{1}{n}(\sum_{i=1}^n X_i)(\sum_{i=1}^n Y_i) \\ (\sum_{i=1}^n X_i^2 - \frac{1}{n}(\sum_{i=1}^n X_i)^2)\beta_1 &= \sum_{i=1}^n X_iY_i - \frac{1}{n}(\sum_{i=1}^n X_i)(\sum_{i=1}^n Y_i) \end{align} \] となるから，\(\beta_1\)は \[ \begin{align} \beta_1 &= \frac{\sum_{i=1}^n X_iY_i - \frac{1}{n}(\sum_{i=1}^n X_i)(\sum_{i=1}^n Y_i)} {\sum_{i=1}^n X_i^2 - \frac{1}{n}(\sum_{i=1}^n X_i)^2} \\ &= \frac{\sum_{i=1}^n(X_i-\frac{1}{n}\sum_{i=1}^n X_i)(Y_i-\frac{1}{n}\sum_{i=1}^n Y_i)} {\sum_{i=1}^n (X_i - \frac{1}{n}\sum_{i=1}^n X_i)^2} \\ &= \frac{\sum_{i=1}^n(X_i-\bar{X})(Y_i-\bar{Y_i})} {\sum_{i=1}^n (X_i - \bar{X})^2} \end{align} \] となる。

一方，\(\beta_0\)の値は正規方程式の上の式を\(n\)で割れば \[ \beta_0 + \frac{1}{n}\sum_{i=1}^n \beta_1 X_i = \frac{1}{n}\sum_{i=1}^n Y_i \] から \[ \begin{align} \beta_0 &= \frac{1}{n}\sum_{i=1}^n Y_i - \frac{1}{n}\sum_{i=1}^n \beta_1 X_i \\ &= \bar{Y} - \beta_1 \bar{X} \end{align} \] として得られる。

参考

経済分析のための統計的方法

---

1. こういう形の誤差関数は少しずつ誤差を下げていく勾配降下法のような数値解でも，いきなり誤差の最小値を算出して使う最小二乗法のような解析解でも解くことができます。↩

2. 切片の係数\(\hat{\beta}_0\)のほうも同じ形の式で検定統計量が算出されます↩