最小二乗法による単回帰の回帰係数の推定

クリスカ^{∗}

2011 年 7 月 27 日

1 回帰分析とは

身長と体重の間には, 一方が増加すれば他方も増加するという関係があることが分かります。このような関係が成り立っているとき, この間の関係式を推定しようというのが回帰分析です. 回帰分析は複数種類のデータを分析するための 1 つであり, あるデータを他のデータの関数によって説明するものです.

2 データの統計量

2 種類のデータ x, y について n 個の観測値 (x_{i}, y_{i}) (i = 1, 2, …, n) が与えられていたとします.

データの平均 : \overline{x} = \frac{1}{n} ∑_{i=1}^{n} x_{i}

データの分散 : (s_{x})^{2} = \frac{1}{n} ∑_{i=1}^{n}(x_{i} − \overline{x})^{2} = \frac{1}{n} ∑_{i=1}^{n} x_{i}^{2} − (\overline{x})^{2} = \overline{x_{i}^{2}} − (\overline{x})^{2}

データの共分散 : s_{xy} = \frac{1}{n} ∑_{i=1}^{n}(x_{i} − \overline{x})(y_{i} − \overline{y}) = \frac{1}{n} ∑_{i=1}^{n} x_{i}y_{i} − \overline{x} \overline{y} = \overline{xy} − \overline{x} \overline{y}

データの相関係数 : r_{xy} = \frac{s_{xy}}{s_{x}s_{y}}

Remark データから平均を差し引き標準偏差で割って, 位置, 尺度の調整をすると, 平均は 0, 標準偏差は 1 に揃ったことになります. この操作をデータの標準化といいます. x_{i} を z_{i} = \frac{x_{i} − \overline{x}}{s_{x}}, y_{i} を w_{i} = \frac{y_{i} − \overline{y}}{s_{y}} と標準化したデータ z_{i}, w_{i} を考えると, 定義より,

r_{xy} = \frac{1}{n} ∑_{i=1}^{n} z_{i}w_{i}

となります. これは標準化したデータ z_{i}, w_{i} の共分散に等しいことが分かります.

∗ k-nakagawa@h6.dion.ne.jp

このことから, −1 ≤ r_{xy} ≤ 1 が言えます. 証明は以下の通りです.

\frac{1}{n} ∑_{i=1}^{n}(z_{i} ± w_{i})^{2} = \frac{1}{n} ∑_{i=1}^{n} (z_{i}^{2} ± 2z_{i}w_{i} + w_{i}^{2})

= \frac{1}{n} ∑_{i=1}^{n} z_{i}^{2} ± \frac{2}{n} ∑_{i=1}^{n} z_{i}w_{i} + \frac{1}{n} ∑_{i=1}^{n} w_{i}^{2}

= 1 ± 2r_{xy} + 1 = 2(1 ± r_{xy})

左辺は常に非負であるので 1 ± r_{xy} ≥ 0 より −1 ≤ r_{xy} ≤ 1 が言えます. また, x′ = ax + b, y′ = cx + d (ac > 0) のように線型変換した相関係数を r′_{xy} とすると線型不変性 r′_{xy} = r_{xy} も言えます.

3 回帰分析

相関係数が 1 (もしくは −1) に近く, 正の相関 (負の相関) が認められるとき, 一方の変数を他方の変数の一次関数として表すことが考えられます. つまり, y を x の関数として考えるならば, y = α + βx というモデルになりますが, このモデルを (線型) 回帰モデル, このモデルを用いた分析を回帰分析といいます. また, パラメータ α, β を回帰係数といいます. 上記のように y を x の関数として考えるとき, x を説明変数 (または独立変数), y を被説明変数 (または従属変数) といい, y を x によって説明することを「y を x に回帰する」といいます. ここでは説明変数が 1 つの場合を考えていますが, 観測データが (x_{1i}, x_{2i}, …, x_{mi}, y_{i}) (i = 1, 2, …, n) として m 種類あり, y = α + β_{1}x_{1} + β_{2}x_{2} + ⋯ + β_{m}x_{m} というモデル (y を x_{1}, x_{2}, …, x_{m} に回帰するモデル) を考えることもできます. このように説明変数が複数ある場合を重回帰と呼び, これに対して説明変数が 1 つである場合を単回帰と呼びます. α, β_{1}, β_{2}, … は観測値を用いて推定することになりますが, 推定された α̂, β̂_{1}, β̂_{2}, … を用いた関係式を回帰式 (もしくは回帰直線) といいます.

Remark 相関係数 r_{xy} は x と y の直線的関係を測る尺度であり, この値から非線形の関係を読みとることはできません.

4 最小二乗法による単回帰の回帰係数の推定

観測値 (x_{i}, y_{i}) (i = 1, 2, …, n) を用いた単回帰の回帰係数の推定について考えることにします. この場合, n 個の観測値に対して, y_{i} = α + βx_{i} (i = 1, 2, …, n) という回帰モデルをあてはめて回帰係数 α, β の推定を行うことになります. しかし, 観測値が 1 つの直線上に完全に乗ることはまずあり得ません. ですので, 回帰式の誤差である y_{i} − (α + βx_{i}) を小さくする α, β を推定値とすることになります. 最小二乗法では, 誤差の平方和である下式 S を最小にする α, β を推定値とします.

S = ∑_{i=1}^{n}{y_{i} − (α + βx_{i})}^{2}

x_{i}, y_{i} は変数ではなく, 既に確定した値をとっていることに注意すると, S は α, β の関数だということが分かります. S は非負値の二次式であるので, これを最小化する α, β は, α と β についての偏導関数が 0 となることが必要です. よって, 次の連立方程式の解が求める推定量 α̂, β̂ となります.

\{ \begin{matrix}\frac{∂S}{∂α} = 0 \\ \frac{∂S}{∂β} = 0 \end{matrix}

⟺ \{ \begin{matrix}−2 ∑_{i=1}^{n}\{y_{i} − (α + βx_{i})\} = 0 \\ −2 ∑_{i=1}^{n} x_{i}\{y_{i} − (α + βx_{i})\} = 0 \end{matrix}

⟺ \{ \begin{matrix}∑_{i=1}^{n} \{y_{i} − (α̂ + β̂x_{i})\} = 0 \\ ∑_{i=1}^{n} x_{i} \{y_{i} − (α̂ + β̂x_{i})\} = 0 \end{matrix}

⟺ \{ \begin{matrix}α̂n + β̂ ∑_{i=1}^{n} x_{i} & = & ∑_{i=1}^{n} y_{i} \\ α̂ ∑_{i=1}^{n} x_{i} + β̂ ∑_{i=1}^{n} x_{i}^{2} & = & ∑_{i=1}^{n} x_{i}y_{i} \end{matrix}

ここで, 式変形によって得られた最後の α̂, β̂ に関する連立方程式を正規方程式といいます. 正規方程式の上式の両辺を n で割って整理すると

α̂ + β̂\frac{1}{n} ∑_{i=1}^{n} x_{i} = \frac{1}{n} ∑_{i=1}^{n} y_{i} ⟺ α̂ + β̂ \overline{x} = \overline{y}
⟺ α̂ = \overline{y} − β̂ \overline{x}

となります. これを下式に代入して, 両辺を n で割って整理すると

(\overline{y} − β̂ \overline{x})\frac{1}{n} ∑_{i=1}^{n} x_{i} + β̂\frac{1}{n} ∑_{i=1}^{n} x_{i}^{2} = \frac{1}{n} ∑_{i=1}^{n} x_{i}y_{i} ⟺ (\overline{y} − β̂ \overline{x}) \overline{x} + β̂\overline{x^{2}} = \overline{xy}

⟺ β̂ (\overline{x^{2}} − (\overline{x})^{2}) = \overline{xy} − \overline{x} \overline{y} ⟺ β̂(s_{x})^{2} = s_{xy}

⟺ β̂ = \frac{s_{xy}}{(s_{x})^{2}}

⟺ β̂ = \frac{r_{xy}s_{y}}{s_{x}}

\overline{y} = α̂ + β̂ \overline{x} より, 回帰直線は観測データの平均 (\overline{x}, \overline{y}) を通ることが分かります. y = α̂ + β̂x を y を被説明変数 (x を説明変数) とする回帰直線といいます. y を x で回帰するともいいます.

また,

y = α̂ + β̂x ⟺ y = (\overline{y} − β̂ \overline{x}) + β̂x

⟺ y − \overline{y} = β̂(x − \overline{x})

⟺ y − \overline{y} = \frac{r_{xy}s_{y}}{s_{x}} (x − \overline{x})

⟺ \frac{y − \overline{y}}{s_{y}} = r_{xy} × \frac{x − \overline{x}}{s_{x}}

と変形できることより (y の標準化) = (相関係数) × (x の標準化) という関係があることも分かります. 次に, 回帰式に x_{i} を代入して得られる値 ŷ_{i} を y_{i} の内挿値といい, y_{i} と ŷ_{i} の差 e_{i} = y_{i} − ŷ_{i} を残差といいます. これを用いると

\{ \begin{matrix}∑_{i=1}^{n} \{y_{i} − (α̂ + β̂x_{i})\} = 0 \\ ∑_{i=1}^{n} x_{i} \{y_{i} − (α̂ + β̂x_{i})\} = 0 \end{matrix} ⟺ \{ \begin{matrix}∑_{i=1}^{n} e_{i} = 0 \\ ∑_{i=1}^{n} e_{i}x_{i} = 0 \end{matrix}

が分かります. これは, e_{i} の平均が 0 であり, また, {e_{i}} と {x_{i}} はベクトルとして直交していることを示しています.

5 多重線型回帰と正規方程式

単回帰のときと同様な方法で求めます. 観測値 (x_{1i}, x_{2i}, …, x_{mi}, y_{i}) (i = 1, 2, …, n) を用いた多重線型回帰の回帰係数の推定について考えることにします. この場合, n 個の観測値に対して, y_{i} = α + β_{1}x_{1i} + β_{2}x_{2i} + ⋯ + β_{m}x_{mi} (i = 1, 2, …, n) という回帰モデルをあてはめて回帰係数 α, β_{1}, β_{2}, …, β_{m} の推定を行うことになります. 回帰式の誤差である y_{i} − (α + β_{1}x_{1i} + β_{2}x_{2i} + ⋯ + β_{m}x_{mi}) を小さくする α, β_{1}, β_{2}, …, β_{m} を推定値とすることになります. 最小二乗法では, 誤差の平方和である下式 S を最小にする α, β_{1}, β_{2}, …, β_{m} を推定値とします.

S = ∑_{i=1}^{n}{y_{i} − (α + β_{1}x_{1i} + β_{2}x_{2i} + ⋯ + β_{m}x_{mi})}^{2}

S は非負値の二次式であるので, これを最小化する α, β_{1}, β_{2}, …, β_{m} は, 各々の変数についての偏導関数が 0 となることが必要です. よって, 次の連立方程式の解が求める推定量 α̂, β̂_{1}, β̂_{2}, …, β̂_{m} となります.

\frac{∂S}{∂α} = \frac{∂S}{∂β_{1}} = \frac{∂S}{∂β_{2}} = ⋯ = \frac{∂S}{∂β_{m}} = 0

⟺ { −2 ∑_{i=1}^{n}{y_{i} − (α + β_{1}x_{1i} + β_{2}x_{2i} + ⋯ + β_{m}x_{mi})} = 0 ; −2 ∑_{i=1}^{n} x_{1i}{y_{i} − (α + β_{1}x_{1i} + β_{2}x_{2i} + ⋯ + β_{m}x_{mi})} = 0 ; −2 ∑_{i=1}^{n} x_{2i}{y_{i} − (α + β_{1}x_{1i} + β_{2}x_{2i} + ⋯ + β_{m}x_{mi})} = 0 ; ⋮ ; −2 ∑_{i=1}^{n} x_{mi}{y_{i} − (α + β_{1}x_{1i} + β_{2}x_{2i} + ⋯ + β_{m}x_{mi})} = 0

⟺ { ∑_{i=1}^{n}{y_{i} − (α + β_{1}x_{1i} + β_{2}x_{2i} + ⋯ + β_{m}x_{mi})} = 0 ; ∑_{i=1}^{n} x_{1i}{y_{i} − (α + β_{1}x_{1i} + β_{2}x_{2i} + ⋯ + β_{m}x_{mi})} = 0 ; ∑_{i=1}^{n} x_{2i}{y_{i} − (α + β_{1}x_{1i} + β_{2}x_{2i} + ⋯ + β_{m}x_{mi})} = 0 ; ⋮ ; ∑_{i=1}^{n} x_{mi}{y_{i} − (α + β_{1}x_{1i} + β_{2}x_{2i} + ⋯ + β_{m}x_{mi})} = 0

⟺ { α̂n + β̂_{1} ∑_{i=1}^{n} x_{1i} + β̂_{2} ∑_{i=1}^{n} x_{2i} + ⋯ + β̂_{m} ∑_{i=1}^{n} x_{mi} = ∑_{i=1}^{n} y_{i} ; α̂ ∑_{i=1}^{n} x_{1i} + β̂_{1} ∑_{i=1}^{n} x_{1i}^{2} + β̂_{2} ∑_{i=1}^{n} x_{1i}x_{2i} + ⋯ + β̂_{m} ∑_{i=1}^{n} x_{1i}x_{mi} = ∑_{i=1}^{n} x_{1i}y_{i} ; α̂ ∑_{i=1}^{n} x_{2i} + β̂_{1} ∑_{i=1}^{n} x_{2i}x_{1i} + β̂_{2} ∑_{i=1}^{n} x_{2i}^{2} + ⋯ + β̂_{m} ∑_{i=1}^{n} x_{2i}x_{mi} = ∑_{i=1}^{n} x_{2i}y_{i} ; ⋮ ; α̂ ∑_{i=1}^{n} x_{mi} + β̂_{1} ∑_{i=1}^{n} x_{mi}x_{1i} + β̂_{2} ∑_{i=1}^{n} x_{mi}x_{2i} + ⋯ + β̂_{m} ∑_{i=1}^{n} x_{mi}^{2} = ∑_{i=1}^{n} x_{mi}y_{i}

ここで, 式変形によって得られた最後の α̂, β̂_{1}, β̂_{2}, …, β̂_{m} に関する連立方程式を正規方程式といいます. 正規方程式は式の数が少ないときは容易に解けますが, 式の数が多いときは電子計算機が必要になります.

Remark 多重線形回帰の正規方程式は行列を用いることにより以下のように書き直すことが出来ます.

{ α̂n + β̂_{1} ∑_{i=1}^{n} x_{1i} + β̂_{2} ∑_{i=1}^{n} x_{2i} + ⋯ + β̂_{m} ∑_{i=1}^{n} x_{mi}  = ∑_{i=1}^{n} y_{i} ; α̂ ∑_{i=1}^{n} x_{1i} + β̂_{1} ∑_{i=1}^{n} x_{1i}^{2} + β̂_{2} ∑_{i=1}^{n} x_{1i}x_{2i} + ⋯ + β̂_{m} ∑_{i=1}^{n} x_{1i}x_{mi} = ∑_{i=1}^{n} x_{1i}y_{i} ; α̂ ∑_{i=1}^{n} x_{2i} + β̂_{1} ∑_{i=1}^{n} x_{2i}x_{1i} + β̂_{2} ∑_{i=1}^{n} x_{2i}^{2} + ⋯ + β̂_{m} ∑_{i=1}^{n} x_{2i}x_{mi} = ∑_{i=1}^{n} x_{2i}y_{i} ; ⋮ ; α̂ ∑_{i=1}^{n} x_{mi} + β̂_{1} ∑_{i=1}^{n} x_{mi}x_{1i} + β̂_{2} ∑_{i=1}^{n} x_{mi}x_{2i} + ⋯ + β̂_{m} ∑_{i=1}^{n} x_{mi}^{2} = ∑_{i=1}^{n} x_{mi}y_{i}

⇕

(\begin{matrix}1 & 1 & ⋯ & 1 \\ x_{11} & x_{12} & ⋯ & x_{1n} \\ x_{21} & x_{22} & ⋯ & x_{2n} \\ ⋮ & ⋮ & ... & ⋮ \\ x_{m1} & x_{m2} & ⋯ & x_{mn} \end{matrix}) (\begin{matrix}1 & x_{11} & x_{21} & ⋯ & x_{m1} \\ 1 & x_{12} & x_{22} & ⋯ & x_{m2} \\ ⋮ & ⋮ & ... & ⋮ \\ 1 & x_{1n} & x_{2n} & ⋯ & x_{mn} \end{matrix}) ( α̂ ; β̂_{1} ; β̂_{2} ; ⋮ ; β̂_{m}) = (\begin{matrix}1 & 1 & ⋯ & 1 \\ x_{11} & x_{12} & ⋯ & x_{1n} \\ x_{21} & x_{22} & ⋯ & x_{2n} \\ ⋮ & ⋮ & ... & ⋮ \\ x_{m1} & x_{m2} & ⋯ & x_{mn} \end{matrix}) (y_{1} ; y_{2} ; ⋮ ; y_{n})

ここで, X = (\begin{matrix}1 & x_{11} & x_{21} & ⋯ & x_{m1} \\ 1 & x_{12} & x_{22} & ⋯ & x_{m2} \\ ⋮ & ⋮ & ... & ⋮ \\ 1 & x_{1n} & x_{2n} & ⋯ & x_{mn} \end{matrix}) とおくと上式は

X^{T} X ( α̂ ; β̂_{1} ; β̂_{2} ; ⋮ ; β̂_{m}) = X^{T} (y_{1} ; y_{2} ; ⋮ ; y_{n})

と変形できるので, 最小二乗法による推定値は

( α̂ ; β̂_{1} ; β̂_{2} ; ⋮ ; β̂_{m}) = (X^{T} X)^{−1}X^{T} (y_{1} ; y_{2} ; ⋮ ; y_{n})

となります. 単回帰のときのように, 各パラメータについて明示的な式を示すことはできませんが, 上式右辺の行列計算をした結果が最小二乗法による推定値となります.