決定係数をPythonで計算する

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

計算の内容

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

決定係数の計算
決定係数の計算

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

 積和:$T_{xy}/x$の偏差平方和: $T_{xx}$で傾き:slope(=2)を計算
 $y$の平均$-$傾き×$x$の平均で切片:intercept(=-1)を計算
 数式$y=2x-1$に$x$の各データを代入して$x$の予測値:predict_xを計算
 測定値:$y-$予測値で残差:resudial_yを計算
 1$-$残差の2乗の合計(=10)/$y$の偏差平方和(=50)で決定係数を計算
 xに[1,7]を代入して予測値を計算

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

用語 変数 意      味
xのリスト $x$ 説明変数$x$の配列NumPyで計算するのでndarray形式で定義
yのリスト $y$ 目的変数yの配列、ndarray形式で定義
xの平均 mean_x $x$の平均。平均はaverageよりmeanを使うことが多い
yの平均 mean_y $y$の平均
xの偏差平方和 t_xx 偏差平方和 回帰分析ではt_xxと表現することが多い
yの偏差平方和 t_yy 偏差平方和 回帰分析ではt_yyと表現することが多い
積和 t_xy 積和 回帰分析ではt_xyと表現することが多い
傾き slope 回帰直線の式での傾き 
切片 intercept 回帰直線の式での切片 
予測値 predict_x 求めた回帰直線の式に$x$を代入した予測値
残差 resudial_y $y$の測定値-予測値
決定係数 r2 相関係数の2乗

計算の前に、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については、さらに高度な分析を想定しているので、決定係数を求める方に重きが置かれているようです。