「自分で数式をコードに落としていって動かす」という作業は非常に勉強になると思ったので，いろんなアルゴリズムをゼロから作っていきたいと思います。第一弾として単回帰からはじめていきます。
理論の要点を整理してから実装する構成で述べていきます。
最小二乗法だけでなく統計的推論についても触れていきます。
モデル
線形回帰（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になる点）で誤差が最小になるパラメータが得られます1。
どちらの誤差関数でも，最小二乗法でのパラメータ推定式は同じものになります。
パラメータの推定
回帰係数のパラメータは「誤差関数（例えば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
同じ値が出ました。
統計的検定
：誤差項は相互に独立に平均0，分散の正規分布に従う
という仮定を置いた単回帰モデルをsimple normal regression modelと呼び，そこでは次のような方法で統計的推論を行うことができます。
標準誤差
推定量\(\hat{\beta}_0, \hat{\beta}_1\)の標準誤差は で，ここでは誤差項の分散の推定量（残差分散）であり，残差平方和を\(n-2\)で割ったもの です。
検定統計量
検定統計量は次のようになります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\)のもとで次の判定を行います。
あるいは\(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-14lmのほうと比較してみます
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} \] として得られる。
参考
- 作者: 岩田暁一
- 出版社/メーカー: 東洋経済新報社
- 発売日: 1983/04/01
- メディア: 単行本
- 購入: 2人 クリック: 25回
- この商品を含むブログ (4件) を見る
