今回と次回でpythonで重回帰分析を実行する方法を二つ紹介します。
今回はscikit-learnのLinearRegressionを使う方法です。
ドキュメントはこちら。
sklearn.linear_model.LinearRegression
最初に検証用のダミーデータを作ります。
$x_{i,j}$を -10 ~ 10の一様分布からサンプリングし、次の式で$y_i$を作ります。
$x_{i,3}$の係数が0になっていることから分かる通り、$x_{i,3}$は$y_i$には何の関係もない値です。
また、ノイズとして正規分布に従う乱数加えておきます。
$$
y_i = 5 + 2x_{i,0} -3x_{i,1} + 4x_{i,2} + 0x_{i,3} + \varepsilon_i,\\
\varepsilon_i \sim N(0,1) , \ \ \ i = 0,1,\cdots, 99
$$
サンプルデータを作って保存するコードがこちら。
import numpy as np
X = np.random.uniform(-10, 10, size=(100, 4))
y = 5 + X@[2, -3, 4, 0] + np.random.normal(0, 1, 100)
np.savez("sample_data.npz", X=X, y=y)
早速、回帰分析して回帰係数と定数項ついでに決定係数を求めてみましょう。
(回帰分析の目的が予測モデルを作ることであれば、データを訓練用と評価用に分けるべきなのですが、
今回は回帰係数を求める方法の紹介が主目的なので分けていません。)
import numpy as np
from sklearn.linear_model import LinearRegression
# データの読み込み
npzfile = np.load("sample_data.npz")
X = npzfile["X"]
y = npzfile["y"]
# モデルのインスタンス生成
model = LinearRegression()
#学習
model.fit(X, y)
# LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)
# 回帰係数
print(model.coef_)
# [ 2.00044153 -2.99255801 3.98231315 0.01708309]
# 切片
print(model.intercept_)
# 4.943464324307898
# 決定係数 R^2
model.score(X, y)
# 0.9988893054003364
各コードの後ろにつけているコメントが僕の環境で実行した時の結果です。
回帰係数も切片もそれぞれほぼ正解に近い値が算出されていますね。
注:サンプルデータを乱数で生成しているので、データ生成からやり直せば結果は変わります。
scikit-learnを使うと、非常に手軽に回帰分析ができることがわかりました。
次回はstatsmodelsで同じことをやってみます。