★データ解析備忘録★

ゆる〜い技術メモ

重回帰のパラメーター推定理論とSAS/IMLによる行列実装

統計学の「基本の「き」」である重回帰を復習しながらSAS/IMLで行列で実装していこうと思います。

理論

重回帰を式で表すと、
{ \displaystyle
y=\beta_0+\beta_1 x_1 +\cdots+\beta_p x_p +\varepsilon
}
となります。
これを行列で表すと
{ \displaystyle
\boldsymbol{Y}=\boldsymbol{X\beta}+\boldsymbol{\varepsilon}
}
要素を書き下すと
{ \displaystyle
\begin{bmatrix}
y_1 \\ \vdots \\ y_n
\end{bmatrix}
=
\begin{bmatrix}
1 & x_{11} & \cdots & x_{1p} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n1} & \cdots & x_{np}
\end{bmatrix}
\begin{bmatrix}
\beta_0 \\ \vdots \\ \beta_p
\end{bmatrix}
+
\begin{bmatrix}
\varepsilon_1 \\ \vdots \\ \varepsilon_n
\end{bmatrix}
}
となります。

最小二乗法による\boldsymbol{\beta}の推定

最小二乗法では、誤差の二乗和を最小化します。つまり、i番目の要素に対して、
{ \displaystyle
S(\boldsymbol{\beta})=\sum_{i=1}^{n}
\varepsilon^2_i
}
を最小化します。
ベクトル表記ではS(\boldsymbol{\beta})=\boldsymbol{\varepsilon}^T\boldsymbol{\varepsilon}を最小化することになるので
{ \displaystyle
\begin{eqnarray}
S(\boldsymbol{\beta})&=& \boldsymbol{\varepsilon}^T\boldsymbol{\varepsilon} \\
&=& (\boldsymbol{Y}-\boldsymbol{X\beta})^T (\boldsymbol{Y}-\boldsymbol{X\beta}) \\
&=& (\boldsymbol{Y}^T-\boldsymbol{\beta}^T \boldsymbol{X}^T)(\boldsymbol{Y}-\boldsymbol{X\beta}) \\
&=& \boldsymbol{Y}^T \boldsymbol{Y} - \boldsymbol{\beta}^T \boldsymbol{X}^T \boldsymbol{Y} - \boldsymbol{Y}^T \boldsymbol{X\beta} + \boldsymbol{\beta}^T \boldsymbol{X}^T \boldsymbol{X\beta}
\end{eqnarray}
}
\boldsymbol{\beta}に関して最小化するための一階の条件は
{ \displaystyle
\begin{eqnarray}
\frac{S(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = -2\boldsymbol{X}^T \boldsymbol{Y} + 2\boldsymbol{X}^T \boldsymbol{X\beta} = 0 \\
\boldsymbol{X}^T \boldsymbol{Y} - \boldsymbol{X}^T \boldsymbol{X\beta} =0 \\
\boldsymbol{X}^T \boldsymbol{Y} = \boldsymbol{X}^T \boldsymbol{X\beta} 
\end{eqnarray}
}
より、求める推定量は
{ \displaystyle
\boldsymbol{\beta} = (\boldsymbol{X}^T \boldsymbol{X})^{-1}\boldsymbol{X}^T \boldsymbol{Y}
}
となります。

最尤法による\boldsymbol{\beta}の推定

回帰式を
{ \displaystyle
\boldsymbol{Y}=\boldsymbol{X\beta}+\boldsymbol{\varepsilon}, \;\; \boldsymbol{\varepsilon} \sim MVN(\boldsymbol{0},\sigma^2\boldsymbol{I})
}
とすると、(\sigma^2は分散、\boldsymbol{I}単位行列)
確率密度関数
f(\boldsymbol{Y}| \boldsymbol{X}, \boldsymbol{\beta}, \sigma^2) = \frac{1}{(2\pi\sigma^2\boldsymbol{I})^{\frac{n}{2}}} \exp \left[-\frac{1}{2\sigma^2}(\boldsymbol{Y}-\boldsymbol{X\beta})^T (\boldsymbol{Y}-\boldsymbol{X\beta}) \right]
なので、尤度は
{ \displaystyle
L(\boldsymbol{\beta},\sigma^2|\boldsymbol{X},\boldsymbol{Y})=f(\boldsymbol{Y}| \boldsymbol{X}, \boldsymbol{\beta}, \sigma^2)
}
であるから、このモデルの対数尤度は
\log L(\boldsymbol{\beta},\sigma^2|\boldsymbol{X},\boldsymbol{Y})= -\frac{n}{2}\log{2\pi\sigma^2} - \frac{1}{2\sigma^2} (\boldsymbol{Y}-\boldsymbol{X\beta})^T (\boldsymbol{Y}-\boldsymbol{X\beta})
となるから、最尤推定量は、
{ \displaystyle
\boldsymbol{\beta} = (\boldsymbol{X}^T \boldsymbol{X})^{-1}\boldsymbol{X}^T \boldsymbol{Y}
}
最小二乗法と同じ結果になりました。

SAS/IMLによる実装

データ

Rにデフォルトで入っている attitude データを使います。
RのデータをSASで使う方法は、以前書いた記事をご覧ください。
y-mattu.hatenablog.com


attitude データは、「管理者に対する態度」というデータセットで、以下のものが int 型で含まれています。

  1. 全般的格付け
  2. 雇用者の苦情の処理
  3. 特別な特権の付与
  4. 学習の機会
  5. 昇進機会
  6. 批判的すぎる
  7. 進歩

f:id:songcunyouzai:20160128140848p:plain

今回は、「全般的格付け」を目的変数に、それ以外を説明変数にしてみます。

SASコードは以下の通りです。

推定した\boldsymbol{\beta}は以下のようになりました。
f:id:songcunyouzai:20160128140935p:plain

これが正しいのか、regプロシージャの推定で確かめます。
f:id:songcunyouzai:20160128135447p:plain
全然有意になってないww
まあでも値はあってました。
今回は単なる勉強のためのコードなのでこれで良しとします。

決定係数とかの行列実装は気が向いたらやります。

それにしてもimlregの行数の差ときたら...
SASのプロシージャにしろRのパッケージにしろ便利ですね。

おまけ:Rによる実装

行列で

data(attitude)
X <- as.matrix(attitude[,-1])
Y <- as.vector(attitude[,1])
X <- cbind(rep(1, nrow(X)), X)
Beta <- solve(t(X)%*%X) %*% t(X) %*% Y
print(Beta)

結果

                  [,1]
           10.78707639
complaints  0.61318761
privileges -0.07305014
learning    0.32033212
raises      0.08173213
critical    0.03838145
advance    -0.21705668

lm関数でやると

> summary(lm(rating~.,data=attitude))

Call:
lm(formula = rating ~ ., data = attitude)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.9418  -4.3555   0.3158   5.5425  11.5990 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 10.78708   11.58926   0.931 0.361634    
complaints   0.61319    0.16098   3.809 0.000903 ***
privileges  -0.07305    0.13572  -0.538 0.595594    
learning     0.32033    0.16852   1.901 0.069925 .  
raises       0.08173    0.22148   0.369 0.715480    
critical     0.03838    0.14700   0.261 0.796334    
advance     -0.21706    0.17821  -1.218 0.235577    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.068 on 23 degrees of freedom
Multiple R-squared:  0.7326,    Adjusted R-squared:  0.6628 
F-statistic:  10.5 on 6 and 23 DF,  p-value: 1.24e-05

当たり前ですが同じ結果です。