【python】正規・非正規分布ナイーブベイズでデータをモデル化する【kaggle,naive bayes,Gaussian,Kernel】

[Zhuang Jiaさんのカーネル](https://www.kaggle.com/jiazhuang)を使って、ナイーブベイズとガウシアンナイーブベイズ、非ガウシアンナイーブベイズを勉強する。

#早速本編

相関が少ない・無いような特徴量では、ナイーブベイズが単純かつ強力な手法であることが実証された。

ナイーブベイズは単純ベイズなどとも呼ばれる。

特徴量Xがあり、目的のYが-1か1で与えられていたとする。
予測元となるデータX0からyを予測したいとする。
条件付き確率で表現すると、

$$
P(y = 1 | X = X_{0})
$$

$$
P(y = -1 | X = X_{0})
$$

を計算して、どちらが確率的に大きい値となるかでyを予測することが出来る。


ではどう計算するのか。


$$
P(y = 1 | X = X_{0}) = \frac{P(X = X_{0}|y = 1) \times P(y = 1) }{ P(X = X_{0})}
$$

を計算する。

(Qiitaのマークダウンで数式書くのしんどい！)


#ナイーブベイズ

特徴量が多数あった時、以下のようになる。(分母は省略)
それぞれの条件付き確率は独立にして考えることができる(このへんが"ナイーブ"と呼ばれる理由らしい)

$$
P(y = 1 | X = X_{0}) = P(X^{(1)} = X_{0}^{(1)}|y = 1) \times P(X^{(2)} = X_{0}^{(2)}|y = 1)・・・・
$$

こんなデータで考える。

```py
import pandas as pd
import numpy as np

data = pd.DataFrame({
    'X1': [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3],
    'X2': ['S', 'M', 'M', 'S', 'S', 'S', 'M', 'M', 'L', 'L', 'L', 'M', 'M', 'L', 'L'],
    'y': [-1, -1, 1, 1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, -1]
})
```

$$
X_{0} = (2,S)
$$

として与えられたとき、y= 1 or -1　を予測する計算をしていく。

まず1と-1の個数をカウントすると、1が9個、-1が6個。
よって、

$$
P(y = 1)=\frac{9}{15} 
$$

$$
P(y = -1)=\frac{6}{15}
$$

y=1であった組み合わせの中を細かく見ていくと、X1=1は2回、X1=1は3回、X1=1は4回登場している。

さらにX2は、X2=Sは1回、X2=Mは4回、X2=Lは4回登場している。

$$
P(y=1|X=X_{0})=P(y=1)⋅P(X1=2|y=1)⋅P(X2=S|y=1)
$$

$$
\frac{9}{15} ・\frac{3}{9}・\frac{1}{9} = \frac{1}{45}
$$

これはy=-1の時でも同じように全体からの割合で計算できる。

$$
P(y=-1|X=X_{0})=P(y=-1)⋅P(X1=2|y=-1)⋅P(X2=S|y=-1)
$$

$$
\frac{6}{15} ・\frac{2}{6}・\frac{3}{6} = \frac{1}{15}
$$

よって、
(2,S)の組み合わせは、y=-1となる可能性が高いことが計算できた。

#ガウシアンナイーブベイズ

離散値の予測じゃなくて連続値の予測をしたい時ってどうするか？
一点の値(たとえばX=1)を推定しようとすると確率は0になってしまうので、X±Δの微小区間で確率を計算するように考えます。


```py
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
from scipy.stats import norm

x = np.linspace(-5, 5)
y = norm.pdf(x)
plt.plot(x, y)
plt.vlines(ymin=0, ymax=0.4, x=1, colors=['red'])
```

![01.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/15958f60-1aad-aa64-1bb5-1796a276f053.png)

X=1の時を正規分布に当てはめて、この時の確率密度を求めることで、その値を求める確率の値とする。

データ2万行200列のデータから、各列の平均と各列の標準偏差を求める。
trainとtestについて計算するので、計200*4のデータフレームになる。


```py
train = pd.read_csv('C:/train.csv', index_col=0)
test = pd.read_csv('C:/test.csv', index_col=0)
target = train.target.values#教師データが入っている
train.drop('target', axis=1, inplace=True)
train.shape, target.shape, test.shape
```


```py
pos_idx = (target == 1)
neg_idx = (target == 0)
stats = []
for col in train.columns:
    stats.append([
        train.loc[pos_idx, col].mean(),
        train.loc[pos_idx, col].std(),
        train.loc[neg_idx, col].mean(),
        train.loc[neg_idx, col].std()
    ])
    
stats_df = pd.DataFrame(stats, columns=['pos_mean', 'pos_sd', 'neg_mean', 'neg_sd'])
stats_df.head()
```

```
 	pos_mean 	pos_sd 	neg_mean 	neg_sd
0 	11.156418 	3.270293 	10.626681 	3.008564
1 	-1.017613 	4.220638 	-1.695770 	4.024813
2 	11.156633 	2.841075 	10.665876 	2.612961
3 	6.864113 	2.070898 	6.788979 	2.040082
4 	11.131337 	1.649266 	11.072412 	1.620103
```

```py
# priori probability 事前確率
ppos = pos_idx.sum() / len(pos_idx) #y=1である確率
pneg = neg_idx.sum() / len(neg_idx) #y=0である確率

def get_proba(x):
    # we use odds P(target=1|X=x)/P(target=0|X=x)　オッズ比を計算する。
    return (ppos * norm.pdf(x, loc=stats_df.pos_mean, scale=stats_df.pos_sd).prod()) /\
           (pneg * norm.pdf(x, loc=stats_df.neg_mean, scale=stats_df.neg_sd).prod())

#norm.pdfで正規分布の確率密度関数に当てはめて、出てきた値をppos,pnegにかける

tr_pred = train.apply(get_proba, axis=1)
#trainのデータに定義した関数をそれぞれ処理させる。
#正規分布なので、列平均に近いほど大きい値になる?
```

こうして計算した予測値と元データを使ってaucを計算する。

```py
from sklearn.metrics import roc_auc_score
roc_auc_score(target, tr_pred)
```

```
0.8899703016513566
```

aucは0.89と高い値になり、正解率も高いモデルとなった。

#非ガウシアンナイーブベイズにしてみる

正規分布を仮定して行ったが、もしデータが非正規であったらどうしたらよいのか？
たとえば1列目のvar_0をつかってplotしてみる。
そもそもこのtrainは正規分布になっているのか？


```py
plt.figure(figsize=(20, 10))
plt.subplot(1, 2, 1)
sns.distplot(train.loc[pos_idx, 'var_0'])
plt.plot(np.linspace(0, 20), norm.pdf(np.linspace(0, 20), loc=stats_df.loc[0, 'pos_mean'], scale=stats_df.loc[0, 'pos_sd']))
plt.title('target==1')
plt.subplot(1, 2, 2)
sns.distplot(train.loc[neg_idx, 'var_0'])
plt.plot(np.linspace(0, 20), norm.pdf(np.linspace(0, 20), loc=stats_df.loc[0, 'neg_mean'], scale=stats_df.loc[0, 'neg_sd']))
plt.title('target==0')
```

![02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/64c2144b-5fd0-8cdf-d6d9-ef7f48bcc077.png)

正規分布とは呼べなさそうな分布をしている。
カーネル密度を使ってフィットさせてみる。


```py
from scipy.stats.kde import gaussian_kde

plt.figure(figsize=(20, 10))
plt.subplot(1, 2, 1)
sns.distplot(train.loc[pos_idx, 'var_0'])
kde = gaussian_kde(train.loc[pos_idx, 'var_0'].values)
plt.plot(np.linspace(0, 20), kde(np.linspace(0, 20)))
plt.title('target==1')
plt.subplot(1, 2, 2)
sns.distplot(train.loc[neg_idx, 'var_0'])
kde = gaussian_kde(train.loc[neg_idx, 'var_0'].values)
plt.plot(np.linspace(0, 20), kde(np.linspace(0, 20)))
plt.title('target==0')
```

![03.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/c4dddd58-4a48-1cb3-35eb-7fc94bac94f6.png)

カーネル関数がかなりフィットさせている。

```py
stats_df['pos_kde'] = None
stats_df['neg_kde'] = None
for i, col in enumerate(train.columns):
    stats_df.loc[i, 'pos_kde'] = gaussian_kde(train.loc[pos_idx, col].values)
    stats_df.loc[i, 'neg_kde'] = gaussian_kde(train.loc[neg_idx, col].values)

def get_proba2(x):
    proba = ppos / pneg
    for i in range(200):
        proba *= stats_df.loc[i, 'pos_kde'](x[i]) / stats_df.loc[i, 'neg_kde'](x[i])
    return proba
```

```py
%%time
get_proba2(train.iloc[0].values)
```

```py
def get_col_prob(df, coli, bin_num=100):
    bins = pd.cut(df.iloc[:, coli].values, bins=bin_num)
    uniq = bins.unique()
    uniq_mid = uniq.map(lambda x: (x.left + x.right) / 2)
    dense = pd.DataFrame({
        'pos': stats_df.loc[coli, 'pos_kde'](uniq_mid),
        'neg': stats_df.loc[coli, 'neg_kde'](uniq_mid)
    }, index=uniq)
    return bins.map(dense.pos).astype(float) / bins.map(dense.neg).astype(float)

tr_pred = ppos / pneg
for i in range(200):
    tr_pred *= get_col_prob(train, i)

roc_auc_score(target, tr_pred)
```

```
0.9087237544204513
```

0.909という高いAUCを出すモデルが出来た。
オーバーフィットしている感はあるが、まあ大丈夫でしょう。
これをtestデータに対して使って終わりとする。


```py
te_pred = ppos / pneg
for i in range(200):
    te_pred *= get_col_prob(test, i)

pd.DataFrame({
    'ID_code': test.index,
    'target': te_pred
}).to_csv('sub.csv', index=False)
```


#まとめ
ナイーブベイズ、ガウシアンナイーブベイズ、非正規のとき
と実践していきました。

#以上

今回はZhuang JiaさんがChris Deotte さんのカーネルを見て実践してみようと思ったものらしいので、Zhuangさんが感銘を受けたというソースも当たってみるつもりです。