StatsModelsとscikit-learnのElastic Net(Lasso回帰)を比較してみた

 最近スパース推定の勉強会に参加しております.使っているテキストはこちら.

スパース推定法による統計モデリング (統計学One Point)

スパース推定法による統計モデリング (統計学One Point)

この本の第2章のElastic Netについて勉強がてらにプログラムを書いてみて,振る舞いを観察してみました.

 今回使うデータはこちら.リンク先のヘッダの「Earn」が観測日の売上.「MaxTemp」が観測日の最高気温で「MinTemp」が最低気温.「Humi」は観測日の湿度.「Wind」は風の強さ.「Sun」は日照時間です.この「Earn」売上に直結する項目は何かを分析することが目的です.
ちなみに,このデータは昔お世話になったアイスクリーム屋さん・ハンバーガ屋さんの統計をイメージして作成したものです.なので実際の現象とは異なりますのであしからず.

各項目の相関係数を出してみるとこんな感じです.「Wind」がわざとらしいほどに無相関ですね.
f:id:shu10038:20180713095238p:plain
各項目の関係をグラフにするとこんな感じ.
f:id:shu10038:20180712223511j:plain
明らかに「Wind」がやる気がない数字が入ってますね.(まぁ,だってそういう風に作ったんだもん)

StatsModelsのエラスティックネットを用いてPythonスクリプトを書いてみた

 統計を用いた分析をする時に「StatsModels」を使うか「scikit-learn」を使うか悩ましかったりします.StatsModelsはPandasのDataFrameを直接使うことができるので,私はどちらかというとStatsModelsを使うことが多いですね.現時点ではscikit-learnはDataFrameを受け付けないですが,そのうち変更があったらどうなるか分かりません.といってもscikit-learn昔書いたこちらの記事を参照にDataFrameからListに変換してあげればいいだけなんですけどね.

StatsModelsのエラスティックネットのドキュメントはこちらこちら

さて,下記のように書いてみました.

この書籍では罰則項を増加させたときの回帰係数を描画した解パス図というのを紹介していました.これに習って同じものを描いてみるとこんな感じ.
f:id:shu10038:20180714090330j:plain
あれ?なんか変だぞ??この場合は一番売上に一番関係なさそうな「Wind」が先に消えるんじゃないの??
なんか変なのでscikit-learnではどういう振る舞いを見せるのかやってみました.

scikit-learnによるエラスティックネット

scikit-learnを使ったエラスティックネットはネット上にサンプルコードが山ほどあるのでここでは特に記載しません.Gitには上げておいたので,興味がある方は下記を参照ください.

さて,同じように解パス図を書くとこんな感じ.
f:id:shu10038:20180714090331p:plain

あれ?StatsModelsとsklearnで結果が違う・・・・
そんな馬鹿な,,という感じですが原因不明です.これでは余りにも内容がない記事になってしまいますので,あとで加筆します.


ではそんな感じで.ではでは.

データフレームの各行の差を計算してグラフ化する:Pandas

やりたいことしては,タイトル通り『データフレームの各行の差を計算してグラフ化する』ことです.Pandasのデータフレームを用いて各列の時間差を計算してグラフに出力するということをやります.

今回はTimedelta型を使っているのですが,データフレームの差を計算して描画しようとすると「'Timedelta' object has no attribute 'plot'」というエラーメッセージ出てきます.無理やりint型になおして出力するとナノ秒になってしまうので,これを上手く工夫して出力しなくてはなりません.

ここからが本題

さて,こんなデータを使います.各駅の名前とそこを通過する時間が記載されているデータがあるとします.

//time.csv
hakata,10:00:05.156
yoshiduka,10:29:54.731
kadomatu,11:03:44.195
katuragawa,11:38:32.395
tendou,12:09:21.151
iiduka,12:46:10.364
kokura,13:21:59.996

このとき前の駅から今の駅まで到着するのにかかった時間を計算してグラフに出力します.

Python3のコード

import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab

# CSVデータの読み込み
df = pd.read_csv('data/time.csv', header=None, names=['sta', 'time'])

# 時間型に変換
df['time'] = df['time'].apply(lambda xx: pd.to_datetime(xx))

# 時間差の計算.各行の差を計算する
df_e =df['time'].diff(1).fillna(0)

# 描画
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(df['sta'], df_e)

# ナノ秒から変換
def toChange(y, i):
    h = int(y / 1000000000 / 3600)
    m = int(y / 1000000000 / 60)
    s = int(y / 1000000000 % 60)
    return '%(h)02d:%(m)02d:%(s)02d' % {'h': h, 'm': m, 's': s}

from matplotlib.ticker import FuncFormatter as ff

ax.yaxis.set_major_formatter(ff(toChange))
plt.grid()

plt.show()

結果

するとこんなふうに出力されます
f:id:shu10038:20180610204723j:plain

Pandasデータフレーム内の複数の文字を「.str.split」を使って区切る

Pandas(@Python)のデータフレーム内のstringを複数の文字で区切りたい場合があります.
「.str.split」を使うのですが少し工夫が必要です.
まぁ,簡単に言うと縦棒「|」で区切るわけですよ.

import pandas as pd
df = pd.DataFrame(
        {'Time': ['jt=10:00:00:-p', 'jt=10:01:01:-p', 'jt=10:10:02:-p', 'jt=10:50:03:-p', 'jt=11:00:04:-p'],
           'Value': [25, 30, 104, 52, 41]})

df_spr = df['Time'].str.split('=|:-', expand=True)

print(df_spr)

出力すると・・

    0         1  2
0  jt  10:00:00  p
1  jt  10:01:01  p
2  jt  10:10:02  p
3  jt  10:50:03  p
4  jt  11:00:04  p

ベイズ更新をやってみた!

3月の末にベイズ推定の勉強会に参加してきました!
math-unknown.connpass.com
この勉強会では単に講師の話を聞くだけではなく,みんなでワイワイディスカッションしながら進めるような形式で,グループ内の誰かが疑問に思ったことを後で調べることでより理解が深まったように思います.
そして講師の方の説明の流暢さに脱帽でした!

その中でもベイズ更新の例としてオオカミ少年の例が興味深いと思いました.勉強会の案内のところに次のリンクが貼ってあったので,共有しておきます.
okandayo.hatenablog.com
ベイズ更新を使ってオオカミ少年の例をシミュレーションしてみようと思いましたが,なかなか時間がなく出来ませんでした.ようやく取り掛かることが出来たのでちょこっとやってみたわけです.

問題設定

次のような風に問題文を作ってみます.「オオカミ少年」のお話は昔NHKの番組で見たことがある程度なのでうろ覚えですが,少年が嘘をついて「オオカミが来たぞ!」と叫んでいるうちに村人が誰も少年を信用しなくなってしまい,本当にオオカミが来た時に村人が食べられてしまった話だったと思います.

問題文;オオカミ少年
ある村では嘘つきは75%の確率で嘘をつき,正直者は85%の確率で本当のことを証言することが知られている.
この村のある少年は事あるごとに「オオカミが来たぞ!」と叫んでいるが,本当に狼が来るときと来ない時があるようだ.この少年は嘘をついている,もしくは時々間違えるようである.この少年が100回叫んだ結果から,この少年は正直者か嘘つきかを答えよ.

問題文から仮定と各値を書き出す

この問題文から,次の様のような仮定を置きます.

  • H_{1}:少年は嘘つき
  • H_{2}:少年は正直者

「オオカミが来たぞ!」と叫ぶイベントに対して,本当にオオカミが来る事象と来ない事象があります.これを次のように置きます.

  • A:少年は嘘をつく(「オオカミが来たぞ!」と叫んだときにオオカミが来ない)
  • \bar { A } :少年は本当のことを言う(「オオカミが来たぞ!」と叫んだときにオオカミが来る)

A\bar { A } は「オオカミが来たぞ!」と叫んで来ないときと実際に来たときに相当します.これらには次の図のような関係があります.
]

さて,次に問題文からそれぞれの値を次のように置きます.

  • P\left( A|{ H_{ 1 } } \right) =0.15:ある村で嘘つきが本当の証言をする確率
    • →嘘つきが「オオカミが来たぞ!」と叫んで,オオカミが来る確率
  • P\left( A|{ H_{ 2 } } \right) =0.75:ある村で正直者が本当の証言をつく確率
    • →正直者が「オオカミが来たぞ!」と叫んでオオカミが来る確率
  • P\left( \overline { A } |{ H_{ 1 } } \right)=0.85 :ある村で嘘つきが嘘の証言をする確率
    • →嘘つきが「オオカミが来たぞ!」と叫んでもオオカミが来ない確率
  • P\left( \overline { A } |{ H_{ 2 } } \right) =0.25:ある村で正直な人が間違えて嘘の証言をする確率
    • →正直者が間違えて「オオカミが来たぞ!」と叫んでオオカミが来ない確率

これを表にするとこんな感じです.

  オオカミが来る A オオカミが来ない  \bar { A }
嘘つき H_{1} 0.15 0.85
正直 H_{2} 0.75 0.25

更新則を求める

さて,次にベイズの定理と全確率の公式から

  1. 「オオカミが来たぞ!」と叫んで本当に狼が来た場合,少年が嘘つきである確率
  2. 「オオカミが来たぞ!」と叫んで狼が来なかった場合,少年が嘘つきである確率

をそれぞれ求めます.ベイズの定理と全確率は下記を参照ください.

「オオカミが来たぞ!」と叫んで本当に狼が来た場合,少年が嘘つきである確率


\begin{eqnarray*} P\left( { H_{ 1 } }|{ A } \right)  & = & \frac { P\left( { A }|{ H_{ 1 } } \right) P\left( { H_{ 1 } } \right)  }{ P\left( { A } \right)  }  \\  & = & \frac { P\left( { A }|{ H_{ 1 } } \right) P\left( { H_{ 1 } } \right)  }{ \sum _{ H }^{  }{ P\left( { A }|{ H } \right) P\left( { H } \right)  }  }  \\  & = & \frac { P\left( { A }|{ H_{ 1 } } \right) P\left( { H_{ 1 } } \right)  }{ P\left( { A }|{ H_{ 1 } } \right) P\left( { H_{ 1 } } \right) +P\left( { A }|{ H_{ 2 } } \right) P\left( { H_{ 2 } } \right)  }  \end{eqnarray*}

「オオカミが来たぞ!」と叫んで狼が来なかった場合,少年が嘘つきである確率


\begin{eqnarray*} P({ H_{ 2 } }|\overline { A } ) & = & \frac { P({ { \overline { A }  } }|{ H_{ 2 } })P\left( { H_{ 2 } } \right)  }{ P({ { \overline { A }  } }) }  \\  & = & \frac { P({ { \overline { A }  } }|{ H_{ 2 } })P({ H_{ 2 } }) }{ \sum _{ H }^{  }{ P({ { \overline { A }  } }|{ H })P\left( { H } \right)  }  }  \\  & = & \frac { P({ { \overline { A }  } }|{ H_{ 2 } })P\left( { H_{ 2 } } \right)  }{ P({ { \overline { A }  } }|{ H_{ 1 } })P\left( { H_{ 1 } } \right) +P({ { \overline { A }  } }|{ H_{ 2 } })P\left( { H_{ 2 } } \right)  }  \end{eqnarray*}

シミュレーション

次に実際にシミュレーションしてみます.Pythonでこんな感じに書きました.

import numpy as np
import matplotlib.pyplot as plt

# 確率の設定
p_AH1 = 0.15 # 嘘つきがが嘘をついてオオカミが来ない確率
p_AH2 = 0.75 # 正直者が本当のことを言ってオオカミが来る確率
p_barAH1 = 0.85 # 嘘つきが間違えて正しいことを言う確率
p_barAH2 = 0.25 # 正直者が間違える確率

# 嘘つきである確率の初期値
p_l = [0.5]

# ループ計算
for i in range(99):
    rand = np.random.random()

    #「オオカミが来たぞ!」と叫んだときに狼が来る
    if rand >0.5:
        p_l.append(p_l[i]*p_barAH1/(p_barAH1*p_l[i]+p_barAH2*(1-p_l[i])))

    #「オオカミが来たぞ!」と叫んだときに狼が来ない
    else:
        p_l.append(p_l[i]*p_AH1/(p_AH1*p_l[i]+p_AH2*(1-p_l[i])))

# 描画        
plt.plot(range(100),p_l)
plt.xlabel('times')
plt.ylabel('H1')
plt.grid()
plt.show()
シミュレーション結果

狼が来る確率を0.4.0.5,0.6とそれぞれ変化させてみたのが下の図です.(ちなみにx軸のレンジが100だと見づらいので80にしています)

今回は問題設定上,正直者と判定されやすいのでこうなりました.
ただし,狼が来る確率を0.4の場合には毎回結果が異なり,正直者と判定される場合や嘘つきと判定される場合に分かれました.その様子を表したのが下の図です.

というわけで今日はこんな感じで.

ATOMをCドライブ以外にインストールする方法

ATOMってエディタがあって,これはでは「こいつはCドライブ以外にインストールできないのか?」と思っていましたが,ひょんなことからDドライブとかにもインストールできました.


ある時ATOMを再インストールしようとすると「Installation has failed」と表示されました.なので下記を参考にさせていただきました.
takumi9942.net

これによれば

  1. AtomSetup.exeあるいはAtomSetup Beta exeを7zipソフトで開く
  2. 中のファイルを短いパスのところに展開する(例:C:\Temp)
  3. コマンドプロンプトを起動し、展開したフォルダまで移動する。(例:コマンドプロンプトを起動して、「cd c:\Temp」と入力する)
  4. 次のコマンドを実行する「Update.exe –install=.」

ということでしたが,この3番めのファイルを展開する際にD:\Tempとかで展開すればDドライブにもインストールできるというわけです.

おしまい.

PandasとStatsModelsを使って重回帰分析をしてみた

Pythonで重回帰分析をやってみました

多変量解析法入門 (ライブラリ新数学大系)

多変量解析法入門 (ライブラリ新数学大系)

この文献のP43~86を参考を参考にしています.

マンガでわかる統計学 回帰分析編

マンガでわかる統計学 回帰分析編

また,回帰分析の全体の流れは上の書籍を参考にしています.

今回はStatsModelsというライブラリを使いました.

データを可視化する

データファイルを読み込んでどんなデータかを観察します.
今回は3次元データなので3次元プロットをしてみます.

# データの読み込み
import pandas as pd

file = 'reg_test_data.csv'
df = pd.read_csv(file)
df = df.set_index('val')
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 初期化
fig = plt.figure()
ax = Axes3D(fig)

# 軸ラベルを設定する
ax.set_xlabel("X-axis")
ax.set_ylabel("Y-axis")
ax.set_zlabel("Z-axis")

# データを無理やりarrayに変換する
xs = df.T.ix[:,0].values
ys = df.T.ix[:,1].values
zs = df.T.ix[:,2].values

# 描画
ax.scatter3D(xs, ys, zs)
plt.show()

f:id:shu10038:20180128193842p:plain
3次元プロットをしてみるとこんな感じです.

f:id:shu10038:20180128193845p:plain
横からみるとこんな感じです.x1とx2には相関はなさそうですが,x1とy,x2とyには相関がありそうです.

散布図行列を描く

描画をしただけだと感覚的にしかわからないので,定量的に評価します.
今回は散布図行列という可視化をしてみました.

from pandas.tools.plotting import scatter_matrix
scatter_matrix(df.T)
plt.show()

f:id:shu10038:20180128193843p:plain

相関係数を計算する

同じく,各変数の相関を定量評価します.

# 相関係数を表示
df.T.corr()
print(df.T.corr())
val        x1        x2         y
val                              
x1   1.000000 -0.170384  0.675085
x2  -0.170384  1.000000  0.603907
y    0.675085  0.603907  1.000000

この通りx1とx2には相関はなさそうですが,x1とy,x2とyには相関がありそうというのがわかります.


重回帰分析の計算をする

ここからが実際の分析です.
今回使ったStatsModelsのドキュメントは下記のとおりです

import statsmodels.formula.api as sm
reg = "y ~ x1 + x2"
model = sm.ols(formula=reg, data=df.T)

# 回帰分析を実行する
result = model.fit()
==============================================================================
Dep. Variable:                      y   R-squared:                       0.988
Model:                            OLS   Adj. R-squared:                  0.987
Method:                 Least Squares   F-statistic:                     703.1
Date:                Sun, 28 Jan 2018   Prob (F-statistic):           4.53e-17
Time:                        19:34:45   Log-Likelihood:                -25.297
No. Observations:                  20   AIC:                             56.59
Df Residuals:                      17   BIC:                             59.58
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
Intercept      0.7622      0.829      0.919      0.371        -0.988     2.512
x1             1.0248      0.034     29.785      0.000         0.952     1.097
x2             1.0218      0.037     27.524      0.000         0.943     1.100
==============================================================================
Omnibus:                        4.291   Durbin-Watson:                   1.543
Prob(Omnibus):                  0.117   Jarque-Bera (JB):                1.439
Skew:                           0.024   Prob(JB):                        0.487
Kurtosis:                       1.687   Cond. No.                         86.1
==============================================================================

回帰式を使ってもう一回描画する

# 先程の結果から値を獲得
b0,b1,b2= result.params

# もう一回描画する
import numpy as np

# 変数の区間の指定
x = np.arange(0, 21, 3)
y = np.arange(0, 21, 3)

# メッシュ表示
X, Y = np.meshgrid(x, y)

# 回帰式を代入する
Z = b0 + b1*X +b2*Y

# 初期化
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_wireframe(X,Y,Z)

# 軸名を設定する
ax.set_xlabel("X-axis")
ax.set_ylabel("Y-axis")
ax.set_zlabel("Z-axis")

# 表示範囲を指定する
ax.set_xlim(-5, 25)
ax.set_ylim(0, 25)
ax.set_zlim(0, 35)

# データを無理やりarrayに変換する
xs = df.T.ix[:,0].values
ys = df.T.ix[:,1].values
zs = df.T.ix[:,2].values

結構上手くフィティングされているように見えます.
f:id:shu10038:20180128193844p:plain

真横から見ると,ピタッと合っている気がします
f:id:shu10038:20180128193846p:plain

PyCharm:interpolaterの設定にまつわるエラー

年末年始にかけてPyCharmをインストールし直したら,エラーが出まくりました・・・
自分が勉強不足のところもあるわけですが,直すのに時間が掛かりました.
色々試行錯誤したのでメモっておきます.
ちなみにインストールしているバージョンは「Anaconda3-5.0.1-Windows-x86_64」「pycharm-community-2017.3.2」です.

そもそもRun出来ない場合

症状

f:id:shu10038:20180104105644j:plain
Runしようとすると

Please Select a valid Python interpolater

と出る.

ちなみにバージョンはCommunity版の2017.3

解決法

stackoverflowに書いてあったのでメモ

  1. 左上の「File」
  2. 「Settings」
  3. 「Project:○○(プロジェクトの名前)」
  4. Project interpolater右のギアマーク
  5. Add Local

ライブラリが読み込めない,importしない

症状

f:id:shu10038:20180113153744p:plain
f:id:shu10038:20180113153745p:plain

  • Pandasとかを読み込もうとすると,「ModuleNotFoundError: No module named ~~」とか出る
  • improt OSなどはできる
  • 上の写真のようにPandasのところだけ「no module name」と出る
  • Jupyter Notebookとかではimportできる

解決法

これは正しい解決法か知りませんが,Projectを作成する際に,「Existing interplater」の値にAnaconda直下のPython3のEXEファイルを代入しました
f:id:shu10038:20180113155841j:plain

結局どうしたかというと・・・

その後一週間くらいしたら,python.exeの実行ファイルがなくなっていたりと,不具合続出でした.
だんだん面倒になってきたので古いバージョン(「Anaconda3-4.2.0-Windows-x86_64」「pycharm-community-2016.1.5」)をダウンロードしてインストールしなしました・・・