NumPyで行列　回帰分析 その6 最小2乗法｜Pythonで数学を学ぼう！ 第27回

線形回帰モデルの設定

回帰分析における回帰直線とは、多数のデータを近似する直線ともいえます。多数のデータから回帰係数を算出する合理的考え方が最小2乗法です。

詳しくその仕組みをみていきます。途中過程の理解が難しいと思われるかもしれませんが、最終的な結果はスッキリした公式に帰着する様子がご覧いただけます。

まずこのあとの重回帰分析の解説のことも踏まえて、変数ならびにインデックスを決めます。データ数をn個、説明変数をx（m個）、目的変数をyとします。

したがって、m個の説明変数x₁、x₂、…、x_mと目的変数yのデータ組がn個あることになります。i番目のデータを変数の上付きとしてxⁱ、yⁱとします。

説明変数xについては、データのi番目を上付き、説明変数のj番目を下付きとして、x_jⁱと表すことにします。

図1：データ表

これにより、目的変数と説明変数を次のように定義できます。回帰式は説明変数の1次式、すなわち線形であるとします。

y=θ₀＋θ₁x₁＋θ₂x₂＋…＋θ_mx_m

ここで、θ₀、θ₁、θ₂、…、θ_mは回帰係数です。

最後にもう一つ変数を設定します。当たり前のことですがデータ表のデータはこの回帰式を満たしません。回帰式はデータの近似式です。そこで両者の差を誤差eとおけばn個のデータについて次のn本の関係式が成り立ちます。

y¹=θ₀＋θ₁x₁¹＋θ₂x₂¹＋…＋θ_jx_j¹＋…＋θ_mx_m¹＋e¹
y²=θ₀＋θ₁x₁²＋θ₂x₂²＋…＋θ_jx_j²＋…＋θ_mx_m²＋e²
⁝
yⁱ=θ₀＋θ₁x₁ⁱ＋θ₂x2ⁱ＋…＋θ_jx_jⁱ＋…＋θ_mx_mⁱ＋eⁱ
⁝
yⁿ=θ₀＋θ₁x₁ⁿ＋θ₂x₂ⁿ＋…＋θ_jx_jⁿ＋…＋θ_mx_mⁿ＋eⁿ

回帰分析における最小2乗法

ここからが最小2乗法の出番です。このまま続けていくと重回帰分析の計算となるので、まずは回帰分析における最小2乗法をみていきます。その後で、再びここから重回帰分析における最小2乗法の計算をみていきます。

説明変数が1つ（m=1）の場合が回帰分析（単回帰分析）です。n個のデータのi番目のデータ（xⁱ, yⁱ）は次の関係式を満たします。

yⁱ=θ₀＋θ₁x₁ⁱ＋eⁱ

この式における差eⁱができるだけ小さくなるような直線が、最もあてはまりのよい（良い近似）直線と考えることができます。

そこで、各データでの差の2乗和(e¹)²＋(e²)²＋…＋(eⁿ)²を最小になるように回帰係数θ₀、θ₁の値を決定します。このような方法を最小2乗法といいます。eⁱ＝yⁱ-(θ₀＋θ₁x₁ⁱ) なので、差の2乗和を

Q＝{y¹-(θ₀＋θ₁x₁¹)}2＋{y²-(θ₀＋θ₁x₁²)}2＋…＋{yⁿ-(θ₀＋θ₁x₁ⁿ)}²

とおきます。

図2：最小2乗法　差の2乗和Qが最小となるように回帰係数を決定する

2乗和Qはθ₀とθ₁の関数とみることができます。そこで、Qの微分をすることでQを最小化するθ₀とθ₁を求めることができます。

回帰直線の傾きθ₁＝〔x₁とyの共分散〕/〔x₁の分散〕＝Cov(x_1,y)/Var(x₁)
回帰直線の切片θ₀＝〔yの平均〕- θ₁×〔x₁の平均〕＝mean(y)-θ₁×mean(x₁)

図3：2乗和Qをθ₀とθ₁で微分する

平均・分散・共分散

2つの回帰係数には代表的統計値である平均mean、分散Var、そして共分散Covがシンプルに現れています。そこで、これらの統計値の計算公式をpythonのコードとともに確認してみます。

図4：平均の公式

簡単のために5個の数値データを用います。np.array()を用いてベクトル・データとして定義します。データの総和にはsum()関数が便利です。引数に配列を渡せば全要素の総和を返します。

データ数はlen(x)として取得できるので、sum(x)/len(x)によりxの平均が得られます。データyの平均についても同様です。NumPyライブラリにはベクトル（配列）の平均を返すnp.mean()関数があります。引数に配列を渡せばいいので便利です。

>>>  import numpy as np
>>>
>>>   x = np.array([1, 2, 3, 4, 5])
>>>   y = np.array([76, 60, 83, 90, 32])
>>>
>>>   print('xの平均')
>>>   mean_x = sum(x)/len(x)
>>>   print(mean_x)
>>>   print(np.mean(x))
>>>
>>>   print('yの平均')
>>>   mean_y = sum(y)/len(y)
>>>   print(mean_y)
>>>   print(np.mean(y))

これを実行すると次のような結果が得られます。
xの平均
3.0
3.0
yの平均
68.2
68.2

続いて分散と共分散です。これらは平均と比べると一般的ではありませんが、統計値として基本的なものです。平均ではデータの散らばり具合はわかりません。

分散はデータの散らばり具合を表す数値です。分散は0以上の値をとります。分散が小さいほどデータは平均値に集まっていることを、逆に大きいほどデータが平均値からばらついていることを表します。

分散の公式（その1）は分散の定義式「分散＝偏差（データと平均の差）の2乗平均」です。pythonのコードでは、xの平均np.mean(x)を用いて「sum((x-np.mean(x))**2)/len(x)」とできます。

もう一つよく使われる分散の公式（その2）「分散＝（2乗平均）－（平均の2乗）」があります。2乗平均はnp.mean(x**2)、平均の2乗はnp.mean(x)**2となります。

そしてNumPyライブラリにはベクトル（配列）の分散を返すnp.var()関数があります。引数に配列を渡せばいいので便利です。これら3つのコードを同時に表示してみます。コードの下に実行結果を示します。

>>>   print('xの分散')
>>>   var1_x = sum((x - np.mean(x))**2)/len(x)
>>>   print(var1_x)
>>>   var2_x = np.mean(x**2) - np.mean(x)**2
>>>   print(var2_x)
>>>   print(np.var(x))
>>>
>>>   print('yの分散')
>>>   var1_y = sum((y - np.mean(y))**2)/len(y)
>>>   print(var1_y)
>>>   var2_y = np.mean(y**2) - np.mean(y)**2
>>>   print(var2_y)
>>>   print(np.var(y))

これを実行すると次のような結果が得られます。
xの分散
2.0
2.0
2.0
yの分散
426.56000000000006
426.5599999999995
426.56000000000006

最後に共分散です。2種類（2変数）のデータの関係を表す数値です。例えば身長と体重の関係を考えてみます。一般に身長が大きいほど体重も大きくなる傾向があります。身長と体重の共分散は正となり、「正の相関がある」といいます。

一方、共分散が負である場合は、一方が増加するともう一方は減少する傾向にあるといえます。これを「負の相関がある」といいます。また、共分散の絶対値が大きいほど相関関係は強いことを表し、共分散が0に近いほど相関関係は小さいといえます。

共分散の公式（その1）は共分散の定義式「共分散＝（xの偏差 × yの偏差）の平均」です。偏差は分散の場合と同じです。pythonのコードでは、xの平均np.mean(x)とyの平均np.mean(y)を用いて「sum((x - np.mean(x))*(y - np.mean(y)))/len(x)」とできます。

もう一つよく使われる共分散の公式（その2）「共分散＝（xyの平均）－（xの平均）×（yの平均）」があります。xyの平均はnp.mean(x*y)とすればいいので「np.mean(x*y) - np.mean(x) * np.mean(y)」と表されます。

そして、NumPyライブラリにはベクトル（配列）の共分散を返すnp.cov()関数があります。引数に配列を2つ渡します。すると結果は2×2の行列（配列）が出力されます。それぞれ次のような成分を意味します。
array([[xの分散 , 共分散],
[共分散 , yの分散]])
回帰分析では共分散だけが必要になるので2行1列成分を次のようにして取り出します。
>>>  np.cov(x, y, bias=True)[0][1]

np.cov()関数を使う時の注意がパラメータbias（初期値False）です。Falseだと偏差の積を（データ数-1）で割ります。Trueにするとデータ数で割ります。Trueに指定する必要があります。これら3つのコードを同時に表示してみます。

>>>  print('xとyの共分散')
>>>  cov1_xy = sum((x - np.mean(x))*(y - np.mean(y)))/len(x)
>>>  print(cov1_xy)
>>>  cov2_xy = np.mean(x*y) - np.mean(x) * np.mean(y)
>>>  print(cov2_xy)
>>>  print(np.cov(x, y, bias=True))
>>>  print(np.cov(x, y, bias=True)[0][1])

これを実行すると次のような結果が得られます。
xとyの共分散
-11.6
-11.600000000000023
[[ 2. -11.6 ]
[-11.6 426.56]]
-11.600000000000001

図5：分散と共分散の公式

回帰係数の算出

以上で回帰家数を求める準備ができました。xとyにそれぞれ身長と体重を配列として定義します。2つの回帰係数の公式から次のようにtheta_1（傾き）とtheta_0（切片）を計算します。

>>>  theta_1 = np.cov(x, y, bias=True)[0][1]/np.var(x)
>>>  theta_0 = np.mean(y) - theta_1 * np.mean(x)

平均・分散・共分散はすべてNumPyの関数を使うことで、回帰係数の公式はそのままpythonで表すことができます。

>>>  import numpy as np
>>>
>>>  x = np.array([165, 170, 172, 175, 170, 172, 183, 187, 180, 185]) # 身長
>>>  y = np.array([52, 60, 63, 65, 78, 70, 72, 85, 90, 94]) # 体重
>>>
>>>  theta_1 = np.cov(x, y, bias=True)[0][1]/np.var(x)
>>>  theta_0 = np.mean(y) - theta_1 * np.mean(x)
>>>
>>>  print(f'回帰直線の傾きθ_1 = {theta_1}')
>>>  print(f'回帰直線の切片θ_0 = {theta_0}')

これを実行すると次のような結果が得られます。
回帰直線の傾きθ_1 = 1.4625684723067556
回帰直線の切片θ_0 = -184.36579427875833

以上のPythonコードは「Pythonで数学を学ぼう！第27回.ipynb」としてまとめてあります。

このプログラムファイルは次からダウンロードできます。
https://drive.google.com/file/d/1GKY45QRANTPobraRXs-VhNv7uC5MEbfQ/view?usp=share_link