回帰直線の式をPythonで計算する

Pythonを使ってサンプルデータから回帰直線の式を作り、残差や予測値を計算します。はじめに、NumPyを使って積和や偏差平方和を手順通りに計算します。次の特殊な関数や機械学習で使うscikit-learnさらには統計解析で使うstatsmodelsによる方法をご紹介します。

計算の内容

表のサンプルの測定値のデータとして使います。はじめに説明変数を$x$、目的変数を$y$として回帰直線の傾きと切片を求め、ここから元のデータ$x$のそれぞれの値について予測値を計算します。次に、測定値と予測値の残差を計算し、最後にデータ$x$が0および6だった場合の予測値を計算します。

詳しい手順は次の通りです。

なお、全体を通して使う変数は次の通りです。

NumPyで計算

準備としてサンプルの$x$、$y$および最後に代入するデータをndarray形式で定義します。

NumPyによるデータの定義

1. import numpy as np
2. x = np.array([1, 4, 3, 5, 2])
3. y = np.array([4, 10, 6, 11, 4])
4. data = np.array([0, 6])

手順通りの計算

はじめに関数を使わず手順通り計算します。NumPyで使うためで$x,y$を定義し、簡便法を使い積和、偏差平方和を求めたのち、回帰直線の式の傾き、切片を計算します予測値や、残差はNumPyのブロードキャストが使えるのですっきりと表現できます。

最後に0と6の場合を予測する時も、predict_dを求める式に配列を代入すると、配列で答えが返ります。

NumPyを使って簡易法により計算

1. n=len(x)
2. t_xy=sum(x*y)-(1/n)*sum(x)*sum(y)
3. t_xx=sum(x**2)-(1/n)*sum(x)**2
4. slope=t_xy/t_xx
5. intercept=(1/n)*sum(y)-(1/n)*slope*sum(x)
6. print('傾き=',slope,'切片=',intercept)
7. predict_x=intercept+slope*x
8. print('予測値=',predict_x)
9. resudial_y=y-predict_x
10. print('残差=',resudial_y)
11. predict_d=intercept+slope*data
12. print('x=[0,6] --> y=',predict_d)

結果は次の通りです。

傾き= 2.0 切片= 1.0
予測値= [ 3.  9.  7. 11.  5.]
残差= [ 1.  1. -1.  0. -1.]
x=[0,6] --> y= [ 1. 13.]

Polyfit関数とpoly1d関数による方法

次にNumPyのpolyfit関数とpoly1dメソッドを使います。これらの概要は次の通りです。polyfit関数の引数に$x$、$y$の配列を指定することで、変数fitに傾きと切片の配列が返ります。

次に、poly1dメソッドを使って回帰直線の式を求めます。といっても、poly1dメソッドは回帰分析に限ったものではなく、引数で指定した多項式のオブジェクトを構築するものです。

ここではpoly1dにより、funcオブジェクトに$y=2x+1$の式がセットされます。そこでfuncの引数として$x$を指定すると$y$の予測値が計算されます。以降は、配列の引き算で残差を求め、最後にdataとして配列[0.6]をfuncに代入し予測値を計算します。

NumPyのpolyfitによる回帰式の計算

1. fit = np.polyfit(x, y, 1)
2. print('[傾き,切片]=', fit)
3. func = np.poly1d(fit)
4. predict_x = func(x)
5. print('予測値=', predict_x)
6. resudial_y = y-predict_x
7. print('残差=', resudial_y)
8. predict_d = func(data)
9. print('x=[0,6] --> y=', predict_d)

これらの関数は、もっと奥深いものがあります。詳しくは次を参照してください。

Numpy Document polyfit関数の説明

Numpy Document poly1d関数の説明

scikit-learn

機械学習でよく使われるサイキットラーン（scikit-learn）による方法です。はじめにlr＝LinearRegression()でscikit-learnにある回帰分析の機能を使うことを宣言します（LinearRegressionクラスのインスタンス化といいます）。次にfitメソッドでxとyの組み合わせについて学習をさせます。このときに、配列x：[1 4 3 5 2]をreshape関数で2次元の配列X[[1] [4] [3] [5] [2]]に変換します。もともとscikitlearnは説明変数が複数あることを前提としているためです。ここでは、X（大文字）としてxと区別します。以降は、predictメソッドで予測していきます。x=0，6の予測値も、[[0] [6]]のように変換してから計算します。

sklearnによる回帰式の計算

1. from sklearn.linear_model import LinearRegression
2. lr = LinearRegression()
3. X = x.reshape((len(x), 1))
4. lr.fit(X, y)
5. print('傾き=',lr.coef_, '切片=', lr.intercept_)
6. predict_x = lr.predict(X)
7. print('予測値=', predict_x)
8. resudial_y = y-predict_x
9. print('残差=', resudial_y)
10. D = data.reshape((len(data), 1))
11. predict_d = lr.predict(D)
12. print('x=[0,6] --> y=', predict_d)

傾き= [2.] 切片= 0.9999999999999982
予測値= [ 3.  9.  7. 11.  5.]
残差= [ 1.  1. -1.  0. -1.]
x=[0,6] --> y= [ 1. 13.]

これまでに比べると少し複雑ですが、scikitlearnは、注目の機械学習が詰め込まれたモジュールなので、機械学習の考え方を垣間見ることができます。

statsmodelsによる方法

最後にstatsmodelsです。まず、add_constant関数で、配列x：[1 4 3 5 2]をX[[1. 1.] [1. 4.] [1. 3.] [1. 5.] [1. 2.]]に変換します。

これは重回帰分析を想定しているためです。

Scikitlearn と同じようにOLSクラス（Ordinary Least Squares regressionの略）でインスタンス化と、predictメソッドで予測します。ただし、predict_xを予測するときに、esults.predict(X)としても問題ありませんが、測定値についてはすでに学習済みなのでXに対する予測値にはfittedvaluesプロパテティ、残差についてはresidプロパティで計算することができます。

statsmodelsによる回帰式の計算

1. import statsmodels.api as sm
2. X = sm.add_constant(x)
3. model = sm.OLS(y, X)
4. results = model.fit()
5. print('[切片,傾き]=', results.params)
6. predict_x = results.fittedvalues # results.predict(X)でも可
7. print('予測値=', predict_x)
8. resudial_y = results.resid
9. print('残差=', resudial_y)
10. D = sm.add_constant(data)
11. predict_d = results.predict(D)
12. print('x=[0,6] --> y=', predict_d)
13. print(results.summary())

[切片,傾き]= [1. 2.]
予測値= [ 3.  9.  7. 11.  5.]
残差= [ 1.00000000e+00  1.00000000e+00 -1.00000000e+00  1.77635684e-15
 -1.00000000e+00]
x=[0,6] --> y= [ 1. 13.]

statsModels は統計モデル用のライブラリで、このほかにも非常に多彩な機能があります。

まとめ

回帰直線の式を求める方法はさまざまですが、scikit-learn、statsmodelsについては、さらに高度な分析を想定しているので、少し手間が煩雑になっています。