決定係数をPythonで計算する

Pythonを使ってサンプルデータから決定係数を計算します。はじめに、NumPyを使って積和や偏差平方和を手順通りに計算します。次に機械学習で使うscikit-learnや統計解析で使うstatsmodelsによる方法をご紹介します。

計算の内容

次の表を測定値のデータとして使います。説明変数を$x$、目的変数を$y$として、決定係数を計算します。

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

なお、プログラムで使う変数は次の通りです。

計算の前に、NumPyを使い、$x$、$y$およびdのデータをndarray形式の配列で定義します。

配列の定義

1. import numpy as np
2. arr_x=np.array([2,3,4,5,6])
3. arr_y=np.array([1,7,8,9,10])
4. arr_d=np.array([1,7])
5. n=len(arr_x)

NumPyで計算

手順通りの計算

はじめに特殊な関数を使わずNumPyで手順通り計算します。

NumPyで相関係数を計算する

1. mean_x = sum(arr_x)/n
2. mean_y = sum(arr_y)/n
3. t_xx = sum((arr_x-mean_x)**2)
4. t_yy = sum((arr_y-mean_y)**2)
5. t_xy = sum((arr_x-mean_x)*(arr_y-mean_y))
6. slope = t_xy/t_xx
7. intercept = (1/n)*sum(arr_y)-(1/n)*slope*sum(arr_x)
8. print('slope=', slope, 'intercept=', intercept)
9. predict_x = intercept+slope*arr_x
10. print('predict=', predict_x)
11. resudial_y = arr_y-predict_x
12. print('resudial=', resudial_y)
13. r2 = 1-(sum(e_y**2))/t_yy
14. print('r_square=', r2)
15. predict_d = intercept+slope*arr_d
16. print('x=[1,7] --> y=', predict_d)

結果は次の通りです。

slope= 2.0 intercept= -1.0
predict= [ 3.  5.  7.  9. 11.]
resudial= [-2.  2.  1.  0. -1.]
r_square= 0.8
x=[1,7] --> y= [ 1. 13.]

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]]）に変換します。もともとscikit-learnは説明変数が複数あることを前提としているためです。ここでは、X（大文字）として$x$と区別します。Irオブジェクトのscoreメソッドにより決定係数を計算します。相関係数は決定係数の平方根で求めます。

scikit-learnによる回帰式の計算

1. from sklearn.linear_model import LinearRegression
2. lr=LinearRegression()
3. X=arr_x.reshape(n,1)
4. lr.fit(X,arr_y)
5. print('slope=',lr.coef_,'intercept=',lr.intercept_)
6. predict_x=lr.predict(X)
7. print('predict=',predict_x)
8. resudial_y=arr_y-predict_x
9. print('resudial=',resudial_y)
10. r2=lr.score(X,arr_y)
11. print('r_square=',r2)
12. D=arr_d.reshape((len(arr_d),1))
13. predict_d=lr.predict(D)
14. print('x=[1,7] --> y=',predict_d)

結果は次の通りです。

slope= [2.] intercept= -1.0000000000000018
predict= [ 3.  5.  7.  9. 11.]
resudial= [-2.  2.  1.  0. -1.]
r_square= 0.8
x=[1,7] --> y= [ 1. 13.]

Xはreshape関数で2次元とする必要があり、次のようなリストになっています。

X= [[2]
 [3]
 [4]
 [5]
 [6]]

さらに、傾きslopeがリストになっています。これらは説明変数が2つ以上あることを想定しています。

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

statsmodelsによる方法

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

statsmodelsのための配列の定義

1. X= [[1. 2.]
2. [1. 3.]
3. [1. 4.]
4. [1. 5.]
5. [1. 6.]]

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

Scikit-learn と同じようにOLSクラス（Ordinary Least Squares regressionの略）でインスタンス化し、rsquaredプロパティで決定係数を計算します。相関係数は決定係数の平方根で求めます。

statsmodelによる回帰式の計算

1. import statsmodels.api as sm
2. X = sm.add_constant(arr_x)
3. model = sm.OLS(arr_y,X)
4. results = model.fit()
5. print('[intercept,slope]=',results.params)
6. predict_x=results.fittedvalues #results.predict(X)でも可
7. print('predict=',predict_x)
8. resudial_y=results.resid
9. print('resudial=',resudial_y)
10. D = sm.add_constant(arr_d)
11. r2=results.rsquared
12. print('r_square=',r2)
13. predict_d=results.predict(D)
14. print('x=[1,7] --> y=',predict_d)

結果は次の通りです。

[intercept,slope]= [-1.  2.]
predict= [ 3.  5.  7.  9. 11.]
resudial= [-2.00000000e+00  2.00000000e+00  1.00000000e+00 -3.55271368e-15
 -1.00000000e+00]
r_square= 0.8
x=[1,7] --> y= [ 1. 13.]

まとめ

統計学のテキストでは相関係数がよくつかわれますが、scikit-learn、statsmodelsについては、さらに高度な分析を想定しているので、決定係数を求める方に重きが置かれているようです。