この記事ではpythonを使って回帰分析する際のコードをまとめていく。いくつか方法がありそうなので、ライブラリごとに書いていければと思う。
今回はstatsmodels.api編。
コードはこちら。詳細は参考サイトを見て頂ければ。今回もメモ書き程度に残していく
%matplotlib notebook import numpy as np import statsmodels.api as sm import matplotlib.pyplot as plt import seaborn as sns #データを作成 x = np.array([1, 2, 3, 4, 5]) y = np.array([2, 6, 6, 9, 6]) #フィッティングモデル import statsmodels.api as sm X = sm.add_constant(x) #説明変数に定数項を追加（切片ありの際必要） re = sm.OLS(y, X).fit() #OLS（最小二乗法）による最適化 #結果 print("回帰係数(b)=", re.params[1]) print("切片(a)：",re.params[0]) #予想値 x_ex = np.arange(-1, 7, 0.1) y_ex = re.params[1] * x_ex + re.params[0] #プロット sns.set_style("whitegrid", {'grid.linestyle':'--'}) plt.rcParams["font.family"] = 'meiryo' fig1 = plt.figure(figsize=(3,3), facecolor='w') ax1 = fig1.add_subplot(111) ax1.set_xlim(-1, 7) ax1.set_ylim(-2, 10) ax1.set_xlabel("X",fontsize=10,fontname='meiryo') ax1.set_ylabel("Y",fontsize=10,fontname='meiryo') ax1.scatter(x, y, marker="o",alpha=1,edgecolors="#08699E",facecolor="#08699E", s=40) #実データ ax1.plot(x_ex, y_ex ,'r--',color="#000000") #回帰の予想値
出力はこんな感じ。
回帰係数(b)= 1.0999999999999999 切片(a)： 2.5000000000000004
予想値を出したいときは、statsmodels.stats.outliers_influence.summary_tableという関数を使うとよい。
from statsmodels.stats.outliers_influence import summary_table st, data, ss2 = summary_table(re, alpha=0.05) y_pred = data[:,2]
このstには下記のようなデータが格納されている。dataには数値データ、ss2には表のラベルが格納されている。
Predicted Valueという列には回帰直線にxを入れた時の予想値が返ってくる。プロットするとこの通り。この出力はデータ順に並んでいのだが、xの値が入ってない。
#プロット ax1.scatter(x, y_pred, marker="o",alpha=1,edgecolors="#08699E",facecolor="#08699E", s=40) #回帰の予想値 ax1.plot(x_ex, y_ex ,'r--',color="#000000") #回帰直線
これが便利なのは、信頼区間も出してくれることである。"alpha"という引数が棄却域の広さを指定している。（0.05であれば0.95、つまり95%が信頼区間になる）信頼区間の詳細な記事はこちらを参考頂きたい。
mean ciが回帰直線の信頼区間である。
#回帰直線の信頼区間をプロット ax1.plot(x, data[:,5], color='#dc0000',alpha=0.5)#信頼区間の上限 ax1.plot(x, data[:,4], color='#dc0000',alpha=0.5)#信頼区間の下限 ax1.fill_between(x, data[:,5], data[:,4],facecolor='#dc0000' ,alpha=0.1)#信頼区間の範囲 ax1.scatter(x,y, marker="o",alpha=1,edgecolors="#999999",facecolor="#999999", s=40)#元データ ax1.plot(x_ex, y_ex ,color="#dc0000") #回帰直線
Predicted ciはデータ点の信頼区間である。
#デ-タ点の信頼区間をプロット ax1.plot(x, data[:,7], color="#019fe6",alpha=0.5)#信頼区間の上限 ax1.plot(x, data[:,6], color="#019fe6",alpha=0.5)#信頼区間の下限 ax1.fill_between(x, data[:,6], data[:,7],facecolor="#019fe6" ,alpha=0.1)#信頼区間の範囲 ax1.scatter(x,y, marker="o",alpha=1,edgecolors="#6276b5",facecolor="#019fe6", s=40)#元データ ax1.plot(x_ex, y_ex , color="#999999",alpha=1)#回帰直線
以前計算した値とちゃんと同じになっていて安心した。
もちろんちゃんと重回帰分析もできる。
%matplotlib notebook import numpy as np import statsmodels.api as sm from statsmodels.sandbox.regression.predstd import wls_prediction_std import matplotlib.pyplot as plt import seaborn as sns import pandas as pd cols = ['X', 'Y', 'Z'] data = pd.DataFrame(index=[], columns=cols) x = np.array([1, 2, 3, 4, 5]) y = np.array([2, 6, 6, 9, 6]) z = np.array([1, 3, 4, 7, 9]) data['X'] = x data['Y'] = y data['Z'] = z X_Y = data.drop('Z', axis=1) Z = data['Z'] #直線回帰モデルを設定 X_Y_c = sm.add_constant(X_Y) #説明変数に定数項を追加（切片ありの際必要） re = sm.OLS(Z, X_Y_c).fit() #OLS（最小二乗法）による最適化 #結果 from statsmodels.stats.outliers_influence import summary_table st, data, ss2 = summary_table(re, alpha=0.05) print(re.params)
こんな感じで切片と回帰係数が得られる。
const -1.160630 X 2.017323 Y -0.015748
ちなみに信頼区間とかも出せるっぽい。このへんはちゃんと勉強してからにしたいので、とりあえず出力だけ載せておく。
from statsmodels.stats.outliers_influence import summary_table st, data, ss2 = summary_table(re, alpha=0.05) print(st)
入力の値であるX,Yの値はこの表に含まれず、あくまでもX,Yのリストの先頭から順に計算して出力されてくる。
参考
Statsmodelの使い方: PythonのStatsmodelsを使って回帰分析を行う - wonderful cool something
