ラビットチャレンジレポート #8 ロジスティック回帰

# ロジスティック回帰モデル
* 分類問題を解くための教師あり機械学習モデル（教師データから学習）
    * 入力とm次元パラメータの整形結合をシグモイド関数に入力
    * 出力は$y=1$になる確率の値になる
* 扱うデータ
    * 入力（説明変数、特徴量）
        * m次元のベクトル(m=1の場合はスカラ)
    * 出力（目的変数）
        * 0 or 1
        * 例）
            * 合格 or 不合格
            * 勝ち or 負け
            * オス or メス
            * 生 or 死

### シグモイド関数
* 入力は実数、出力は必ず0~1の値
* 確率を表現
* 単調増加関数

```math
\sigma(x) = \frac{1}{1+exp(-ax)}
```

#### パラメータが変わるとシグモイド関数の形が変わる
* aを増加させると、x=0付近での曲線の勾配が増加
* aを極めて大きくすると、単位ステップ関数に近づく（x<0でf(x)=0、x>0でf(x)=1となるような関数）

####　シグモイド関数の性質
* シグモイド関数の微分は、シグモイド関数自身で表現することが可能
* 尤度関数の微分を行う際にこの事実を利用すると計算が容易

```math
\frac{d\sigma(x)}{dx}= \frac{d}{dx}(\frac{1}{1+exp(-ax)}) \\
= (-1)・(1+exp(-ax))^{-2}・exp(-ax)・(-a) \\
= \frac{a・exp(-ax)}{(1+exp(-ax))^2} \\
= \frac{a}{1+exp(-ax)}・\frac{1+exp(-ax)-1}{1+exp(-ax)} \\
= a \sigma(x)(1- \sigma(x))
```

#### シグモイド関数の出力をY=1になる確率に対応させる
* データの線形結合を計算
* シグモイド関数で算出した値が0.5以上であれば1、未満なら0と予測する

```math
P(Y=1|x) = \sigma(w_0+w_1x_1+・・+w_mx_m)
```

* 左辺
    * 説明変数の実現値が与えられた際にY=1になる確率
    * データが与えられたときにY=1となる確率

* 右辺
    * データのパラメータに対する線形結合
    * $w_0$が切片
    * $w_1$が回帰係数
    * $x_1$が説明変数

## 最尤推定
### 世の中の様々な確率分布
* 正規分布、t分布、ガンマ分布・・・
* ロジスティック回帰モデルではベルヌーイ分布を利用する

### ベルヌーイ分布
* 数学において確率$p$で1、確率$1-p$で0をとる、離散確立分布（例：コイン投げ）
* 「生成されるデータ」は分布のパラメータによって異なる

### ベルヌーイ分布のパラメータ推定
* ある分布を考えたとき、そのパラメータによって生成されるデータは変化する
* データからそのデータを生成したであろう尤もらしい分布を推定したい→最尤推定

##### 同時確率
* あるデータが得られたとき、それが同時に得られる確率
* 確率変数は独立であることを仮定すると、それぞれの確立の掛け算となる

##### 尤度関数
* データは仮定し、パラメータを変更する
* 尤度関数を最大化するようなパラメータを選ぶ方法を最尤推定という

#### ロジスティック回帰モデルの最尤推定
* 確率pはシグモイド関数となるため、推定するパラメータは重みパラメータとなる
* $(x_1, y_1), (x_2, y_2).・・・.(x_n, y_n)$を精製するに至ったもっともらしいパラメータを探す

##### 尤度関数を最大とするパラメータを探す（最尤推定）
* 対数をとると微分の計算が簡単
    * 同時確立の積が和に変換可能
    * 指数が積の演算に変換可能
* 対数尤度関数が最大になるテント尤度関数が最大になる点は同じ！
    * 対数関数は単調増加（ある尤度の値が$x_1 < x_2$の時必ず、$logx_1 < logx_2$となる）
* 尤度関数にマイナスをかけたものを最小化し、最小二乗法の最小化と合わせる

## 勾配降下法
* 反復学習によりパラメータを逐次的に更新するアプロサーチの一つ
* 学習率と呼ばれるハイパーパラメータでモデルのパラメータの収束しやすさを調整

```math
w_{k+1} = w_k - \eta \frac{dE(w)}{dw}
```

### 確率的勾配降下法(SGD)
* データを1つずつランダムに選んでパラメータを更新
* 勾配降下法でパラメータを1回更新するのと同じ計算量でパラメータをn回更新できるので効率よく最適な解を探索可能

```math
w_{k+1} = w_k - \eta(y_i-p_i)x_i
```

### 分類の評価方法
##### 再現率
* 本当にPositiveなものの中からPositiveと予測できる割合
* 誤りが多少多くても抜け漏れは少ない予測をしたい場合に利用

##### 適合率
* モデルがPositiveと予測したものの中で本当にPositiveである割合
* 見逃しが多くてもより正確な予測をしたい歳に利用

##### F値
* 再現率と適合率の調和平均

## 実装演習
* タイタニックのデータを使用して演習を行う！
* データは下記のような感じ

```python
titanic_df.head(5)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2183979/31e043c5-3be6-722b-891b-bc890856d1ce.png)

* 今回、使わないデータを削除

```python
titanic_df.drop(['PassengerId', 'Name', 'Ticket', 'Cabin'], axis=1, inplace=True)
```

* 年齢（Age）が欠損しているところを年齢の中央値にする

```python
titanic_df['Age'].fillna(titanic_df['Age'].mean(), inplace=True)
```

#### チケット価格から生死を予測
* data1が運賃のみのリスト
* lanel1が生死フラグのみのリスト

```python
data1 = titanic_df.loc[:, ["Fare"]].values
label1 =  titanic_df.loc[:,["Survived"]].values

from sklearn.linear_model import LogisticRegression
model=LogisticRegression()
model.fit(data1, label1)
```

* 最大値、平均、最小値から予測してみる

```python
data1.mean()
# 32.204207968574636

data1.max()
# 512.3292
```

```python
model.predict([[30]])
# array([0])

model.predict([[600]])
# array([1])
```

```python
w_0 = model.intercept_[0]
w_1 = model.coef_[0,0]

def sigmoid(x):
    return 1 / (1+np.exp(-(w_1*x+w_0)))

x_range = np.linspace(-1, 500, 3000)

plt.figure(figsize=(9,5))
plt.legend(loc=2)

plt.plot(data1,np.zeros(len(data1)), 'o')
plt.plot(data1, model.predict_proba(data1), 'o')
plt.plot(x_range, sigmoid(x_range), '-')
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2183979/67587e5b-0df3-58f5-7c3c-4ebcbee0fd91.png)

#### 次は2変数で生死を予測
* 使う変数は先ほど中央値で補完したAgeと下記で作成したPclass_Gender

```python
titanic_df['Gender'] = titanic_df['Sex'].map({'female': 0, 'male': 1}).astype(int)
titanic_df['Pclass_Gender'] = titanic_df['Pclass'] + titanic_df['Gender']

```

```python
np.random.seed = 0

xmin, xmax = -5, 85
ymin, ymax = 0.5, 4.5

index_survived = titanic_df[titanic_df["Survived"]==0].index
index_notsurvived = titanic_df[titanic_df["Survived"]==1].index

from matplotlib.colors import ListedColormap
fig, ax = plt.subplots()
cm = plt.cm.RdBu
cm_bright = ListedColormap(['#FF0000', '#0000FF'])
sc = ax.scatter(titanic_df.loc[index_survived, 'Age'],
                titanic_df.loc[index_survived, 'Pclass_Gender']+(np.random.rand(len(index_survived))-0.5)*0.1,
                color='r', label='Not Survived', alpha=0.3)
sc = ax.scatter(titanic_df.loc[index_notsurvived, 'Age'],
                titanic_df.loc[index_notsurvived, 'Pclass_Gender']+(np.random.rand(len(index_notsurvived))-0.5)*0.1,
                color='b', label='Survived', alpha=0.3)
ax.set_xlabel('Age')
ax.set_ylabel('Pclass_Gender')
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
ax.legend(bbox_to_anchor=(1.4, 1.03))
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2183979/6dc99bfc-001b-26a7-8b00-0adfad66a758.png)


```python
data2 = titanic_df.loc[:, ["Age", "Pclass_Gender"]].values

h = 0.02
xmin, xmax = -5, 85
ymin, ymax = 0.5, 4.5
xx, yy = np.meshgrid(np.arange(xmin, xmax, h), np.arange(ymin, ymax, h))
Z = model2.predict_proba(np.c_[xx.ravel(), yy.ravel()])[:, 1]
Z = Z.reshape(xx.shape)

fig, ax = plt.subplots()
levels = np.linspace(0, 1.0)
cm = plt.cm.RdBu
cm_bright = ListedColormap(['#FF0000', '#0000FF'])
#contour = ax.contourf(xx, yy, Z, cmap=cm, levels=levels, alpha=0.5)

sc = ax.scatter(titanic_df.loc[index_survived, 'Age'],
                titanic_df.loc[index_survived, 'Pclass_Gender']+(np.random.rand(len(index_survived))-0.5)*0.1,
                color='r', label='Not Survived', alpha=0.3)
sc = ax.scatter(titanic_df.loc[index_notsurvived, 'Age'],
                titanic_df.loc[index_notsurvived, 'Pclass_Gender']+(np.random.rand(len(index_notsurvived))-0.5)*0.1,
                color='b', label='Survived', alpha=0.3)

ax.set_xlabel('Age')
ax.set_ylabel('Pclass_Gender')
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
#fig.colorbar(contour)

x1 = xmin
x2 = xmax
y1 = -1*(model2.intercept_[0]+model2.coef_[0][0]*xmin)/model2.coef_[0][1]
y2 = -1*(model2.intercept_[0]+model2.coef_[0][0]*xmax)/model2.coef_[0][1]
ax.plot([x1, x2] ,[y1, y2], 'k--')
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2183979/40b16d38-acc1-4b8a-3540-78e9015a367e.png)


### モデル評価(クロスバリデーション)

```python
from sklearn.model_selection import train_test_split

traindata1, testdata1, trainlabel1, testlabel1 = train_test_split(data1, label1, test_size=0.2)
traindata1.shape
trainlabel1.shape

traindata2, testdata2, trainlabel2, testlabel2 = train_test_split(data2, label2, test_size=0.2)
traindata2.shape
trainlabel2.shape

data = titanic_df.loc[:, ].values
label =  titanic_df.loc[:,["Survived"]].values
traindata, testdata, trainlabel, testlabel = train_test_split(data, label, test_size=0.2)
traindata.shape
trainlabel.shape

eval_model1=LogisticRegression()
eval_model2=LogisticRegression()

predictor_eval1=eval_model1.fit(traindata1, trainlabel1).predict(testdata1)
predictor_eval2=eval_model2.fit(traindata2, trainlabel2).predict(testdata2)

eval_model1.score(traindata1, trainlabel1)
# 0.6573033707865169

eval_model1.score(testdata1,testlabel1)
# 0.6927374301675978

eval_model2.score(traindata2, trainlabel2)
# 0.776685393258427

eval_model2.score(testdata2,testlabel2)
# 0.776536312849162
```


### モデル評価(混合行列)

```python
from sklearn.metrics import confusion_matrix
confusion_matrix1=confusion_matrix(testlabel1, predictor_eval1)
confusion_matrix2=confusion_matrix(testlabel2, predictor_eval2)
```

```python
confusion_matrix1
# array([[110,   9],
#        [ 46,  14]])

confusion_matrix2
# array([[99,  9],
#        [31, 40]])
```

```python
fig = plt.figure(figsize = (7,7))
#plt.title(title)
sns.heatmap(
    confusion_matrix1,
    vmin=None,
    vmax=None,
    cmap="Blues",
    center=None,
    robust=False,
    annot=True, fmt='.2g',
    annot_kws=None,
    linewidths=0,
    linecolor='white',
    cbar=True,
    cbar_kws=None,
    cbar_ax=None,
    square=True, ax=None, 
    #xticklabels=columns,
    #yticklabels=columns,
    mask=None)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2183979/c29cb82a-2583-311f-74cf-bcef3dbf66b2.png)

```python
fig = plt.figure(figsize = (7,7))
#plt.title(title)
sns.heatmap(
    confusion_matrix2,
    vmin=None,
    vmax=None,
    cmap="Blues",
    center=None,
    robust=False,
    annot=True, fmt='.2g',
    annot_kws=None,
    linewidths=0,
    linecolor='white',
    cbar=True,
    cbar_kws=None,
    cbar_ax=None,
    square=True, ax=None, 
    #xticklabels=columns,
    #yticklabels=columns,
    mask=None)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2183979/afc1772f-7625-4025-8761-3fef54c21b39.png)