メモランダム!!

自分用の端書のため,他の人が読めるようには書いていません.悪しからずm(_ _)m

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

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の場合には毎回結果が異なり,正直者と判定される場合や嘘つきと判定される場合に分かれました.その様子を表したのが下の図です.

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

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」)をダウンロードしてインストールしなしました・・・

Pandas:グループ毎に括って最大の値を含む列を抜き出す

PythonのライブラリーであるPandasを使って,「グループ毎に括って最大の値を含む列を抜き出す」方法のメモです.

対象とするのはこんなデータ 

Sensor Time Value
0 T-A 10:00:00 25
1 T-B 10:00:01 30
2 T-C 10:00:02 104
3 T-B 10:00:03 52
4 T-C 10:00:04 41
5 T-A 10:00:05 91
6 T-C 10:00:06 102
7 T-B 10:00:07 40
8 T-B 10:00:08 101
9 T-C 10:00:09 97

 
3種類のSensorから時間と検知された値(Value)が返ってきます.
Sensorの種類ごとに括って,その中で最大の値を出す時間(と言うか行)を抜き出す,というのが今回の狙いです.
 
コードは下記の通り.
特定の列に含まれる最大値を含む行を抜き出す(日本語が複雑・・)方法が分からなかったのでやや力技を使いました.

import pandas as pd

# データ生成
df = pd.DataFrame(
        {'Time': ['10:00:00', '10:00:01', '10:00:02', '10:00:03', '10:00:04','10:00:05', '10:00:06', '10:00:07', '10:00:08', '10:00:09'],
         'Sensor': ['T-A', 'T-B', 'T-C', 'T-B', 'T-C', 'T-A', 'T-C', 'T-B' ,'T-B' ,'T-C'],
           'Value': [25, 30, 104, 52, 41, 91, 102, 40 ,101 ,97]})

# datetime型に変換
df['Time'] = df['Time'].apply(lambda dd: pd.to_datetime(dd))

# 'Sensor'で括る    
df_g = df.groupby('Sensor')

# 'Sensor'でくくられたデータフレームの中で'Value'に最大値を含む行を抜き出す
def select(xx):
    # 'Value'に最大値を含む行を抜き出す(そういうメソッドがあるのかもしれないけど分からなかった)
    val_r = xx[xx['Value'] == max(xx['Value']) ]
    
    # 全く同じ行があった場合は削除(このデータの場合は無いですけど)
    val_r = val_r.drop_duplicates()
    
    return val_r

df_new = df_g.apply(select)

ご覧の通り,Valueに含まれる最大値を含む行は無理やり書きました.
結果は下記の通り
 

Sensor Time Value
Sensor 
T-A 5 T-A 2017-10-08 10:00:05 91
T-B 8 T-B 2017-10-08 10:00:08 101
T-C 2 T-C 2017-10-08 10:00:02 104

Anacondaをインストールしようとすると「Failed to create Anaconda menus」と言われる

Anacondaをアップデートしたらおかしくなったので再インストールした.

そしたら「Failed to create Anaconda menus」と言われてインストールできない.

shirabeta.net

アンインストールしてもフォルダが残っていたので削除したらインストールできたっぽい

Python:秒数が小数点以下のパース

Pandasでデータフレームに格納するときに,時間にナノ秒まで含まれている時のパースのやり方についてのメモ.

対象とするのはこんなデータ.

//file.csv
17:22:59.703371360,10
17:22:59.788956621,20
17:22:59.790719017,30
17:22:59.813919277,20
17:22:59.891942610,10
17:22:59.898820371,20
17:22:59.919604329,30

こんな感じで,小数点以下がすごい長いデータをパースしたいわけです.
一例としてこんな感じでしょうか.

import pandas as pd

# データ読み込み
df = pd.read_csv('file.csv', names=('time', 'val'))

# datetime型に変換
from dateutil import parser
df.time = df.time.apply(lambda x: parser.parse(x))

# 表示
print(df.dtypes)

出力すると

time    datetime64[ns]
val              int64
dtype: object

parserを使って無理やりdatetime型に変えています.
あとはapplyで一行づつ適用してます.
若干力技な気がしますが...

追記

このフォーマットの場合はparse_datesで指定するだけでも行けるっぽい.

import pandas as pd

# データ読み込み
df = pd.read_csv('time.csv', names=('time', 'val'),parse_dates=['time'])

# 表示
print(df.dtypes)

ランダム関数について

numpyランダム関数を色々使ってみました.

まずはインポート

In [1]: import numpy as np


random()は0から1までの値が得られる

In [10]: np.random.random()
Out[10]: 0.45226808024834264


引数を入れると,入れた数字分の乱数が得られる

In [12]: np.random.random(5)
Out[12]: array([ 0.93385679,  0.6383795 ,  0.20373784,  0.76812241,  0.42309771])


範囲指定は出来ないみたいです

In [13]: random.uniform(1,100)
Traceback (most recent call last):


random.uniformで値を指定できるみたいです

In [14]: np.random.uniform(1,100)
Out[14]: 8.543481477277378

In [17]: np.random.uniform(2,5)
Out[17]: 2.2877891814898934


引数を入れないとrandom()と同じなんですかね?

In [15]: np.random.uniform()
Out[15]: 0.4615576165499794


範囲を指定して,個数も指定できるみたいです

In [18]: np.random.uniform(2,5,3)
Out[18]: array([ 4.67283745,  2.26006292,  4.18645001])