ノート/本

回帰分析

たとえば、図XXは、1930年に取られたアメリカの自動車の速度(speed、単位はmph)と 停止距離(dist、単位はフィート)のデータを、散布図に示したものです。 この図を見ると、速度が速いと停止距離も長いという関係がありそうです。 もしその関係が成り立っているのなら、速度XXmphの時には停止距離がYYフィートだろう という予想ができます。このような予測をするために速度と停止距離の関係を求める、というのが 回帰分析でやろうとしていることです。

cars.png

<多変量解析という言葉? いつ説明?>

単回帰分析

単回帰分析は、1つの量的変数の値から他の1つの量的変数の値を予測する手法です。 予測の元としてこちらから値を与える変数を独立変数または説明変数と呼び、 予測の結果値が得られる変数を従属変数または目的変数と呼びます。上記の予測の例では、 速度を説明変数に、停止距離を目的変数に取ります。(注:この例では、速度と停止距離の 関係は対称で、停止距離から速度を推測するモデルも作れます。)

この問題を解くときに、よくグラフの上で「それらしい」直線を引いて、引いた直線上の 欲しい速度(横軸)に当たる点の停止距離(縦軸)の値を読むことがあります。単回帰分析は 基本的にその作業を行います。ただし、「それらしい」直線をきちんと理論的に求めます。

それぞれの説明変数(速度)をxi、その時の観測値(停止距離)をyiとします。つまり、観測 したデータは (x0, y0), (x1, y1), (x2, y2) ... (xn, yn) のn個のペアからなっているわけです。 モデルは直線なので、 y = ax + b で表すことができます。観測したデータにいちばんよくあてはまる、それらしい直線 y = ax + b の係数aとbを求めたいわけです。これが求まれば、ある速度xtの時の 停止距離yt は yt = at x + bによって求まります。

このaとbを測定値 (x0, y0), (x1, y1), (x2, y2) ... (xn, yn) から求めるには、次のように 考えます。速度xiの時の、式から予測される停止距離 yhat は yhat = axi + b なので、 誤差を、誤差=実際に掛かった停止距離ー予測される停止距離= yi - (axi + b) となります。 全てのデータの誤差の絶対値を小さくするために、それぞれの誤差の二乗の和 <式> sum ei2 = sum ( yi - (axi + b) )2 を最小にするような a と b を計算します。この方法は「最小二乗法」と呼ばれています。

<コラム:計算法>〜〜後から埋める。aとbでそれぞれ偏微分して=0として整理する。

pandasの環境で単回帰分析を行うには

〜〜〜

# -*- coding: utf-8 -*-
import pandas as pd        # pandasライブラリを読み込み、それをpdと名付ける
import statsmodels.api as sm
from rpy2.robjects import r, pandas2ri
pandas2ri.activate()
cars = r['cars']
# 単回帰分析
results = sm.OLS(cars['dist'], sm.add_constant(cars['speed'])).fit()
   # add_constantについては本文中で説明
print(results.summary())    # 回帰分析の結果の概要の印刷
b, a = results.params         # 回帰分析で得られた傾きaと切片bの抽出 
# 結果のグラフ表示
plt.plot(cars['speed'], cars['dist'], 'o')   # 元データの散布図描画
x = cars['speed']
plt.plot(x, a*x+b)    # xをspeedとしたときの、直線y=ax+bの描画
plt.show()

~~
http://wcs.hatenablog.com/entry/2016/11/08/231703
add_constantという関数は1という値だけの列を追加していて、これが切片(説明変数にかかわらずオフセットされる量)となる。 Statsmodelsの流儀で、モデルに切片を使う場合はこのようにしなければならない。

~~
http://qiita.com/yubais/items/24f49df99b487fdc374d
ここで ϕ0(x)≡1ϕ0(x)≡1 なので、この行列の1列目はすべて 1 となる。先のコードで「おまじない」と称して np.repeat(1, nsample) という列を加えたのはこのような理由による。

~~

                            OLS Regression Results
==============================================================================
Dep. Variable:                  speed   R-squared:                       0.896
Model:                            OLS   Adj. R-squared:                  0.894
Method:                 Least Squares   F-statistic:                     423.5
Date:                  月, 19  6月 2017   Prob (F-statistic):           9.23e-26
Time:                        14:57:35   Log-Likelihood:                -153.74
No. Observations:                  50   AIC:                             309.5
Df Residuals:                      49   BIC:                             311.4
Df Model:                           1
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
dist           0.3081      0.015     20.578      0.000       0.278       0.338
==============================================================================
Omnibus:                       13.356   Durbin-Watson:                   1.280
Prob(Omnibus):                  0.001   Jarque-Bera (JB):               14.098
Skew:                          -1.137   Prob(JB):                     0.000868
Kurtosis:                       4.262   Cond. No.                         1.00
==============================================================================

~~

得られた回帰直線 y=ax+b の傾きaとbの値は

print('a=', a, 'b=', b)

のようにして表示すると

a= 3.93240875912 b= -17.5790948905

となっていることがわかります。

~~

cars-regression.png

更に、x (speed) のそれぞれの値に対する残差(実測値ー回帰直線での予想値)は、 results.residとして得られる。results.residはxの個数分だけ残差が並んだデータフレーム になっています。それを用いて、xに対する散布図と、残差の値のヒストグラムを、 次のようにして図示できます。

# residualsの散布図を表示
plt.plot(cars['speed'], results.resid, 'o')   # speedと残差results.residの散布図
plt.show()

# residualsのヒストグラムを表示
results.resid.hist(bins=20)    # results.residのヒストグラム
plt.show()

cars-resid-scatter.png

cars-resid-hist.png

また、決定係数(寄与率)は

~~
Rsquared-definition.png

~~
のように定義されるものですが、results.rsuared として求められています。

print('Rsquared', results.rsquared)

の出力は

Rsquared 0.651079380758

のように得られることがわかります。決定率は*****

重回帰分析

Multiple Regression using Statsmodels


添付ファイル: fileRsquared-definition.png 1件 [詳細] filecars-resid-scatter.png 1件 [詳細] filecars-resid-hist.png 1件 [詳細] filecars-regression.png 1件 [詳細] filecars.png 3件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2017-06-20 (火) 09:13:23 (6d)