#1.ロジスティック回帰モデル
$ x = (x_1,x_2,\cdots,x_m)^T \in \mathbb{R}^m $
$ y = (0,1) $
#2.シグモイド関数
●入力は実数・出力は必ず0~1の値
●(クラス1に分類される)確率を表現
●単調増加関数
公式
\sigma (x) = \frac{1}{1 + exp^{(-ax)}}
●aを増加させると，x=0付近での曲線の勾配が増加する
●aを極めて大きくすると，単位ステップ関数(x<0でf(x)=0，x>0でf(x)=1となるような関数)に近づく
●バイアス変化は段差の位置
●シグモイド関数の微分は、シグモイド関数自身で表現することが可能
●尤度関数の微分を行う際にこの事実を利用すると計算が容易
シグモイド関数の微分
\begin{eqnarray}
\frac{\partial \sigma (x)}{\partial \sigma} &=& \frac{\partial}{\partial x}
\Big(\frac{1}{1 + exp^{(-ax)}} \Big) \\
&=& (-1) \cdot \{ 1 + exp^{(-ax)} \}^{-2} \cdot exp(-ax) \cdot (-a) \\
&=& \frac{a \ exp(-ax)}{\{ 1 + exp(-ax) \}^2} \\
&=& \frac{a}{1+exp(-ax)} \cdot \frac{1+exp(-ax)-1}{1+exp(-ax)} \\
&=& a \sigma (x) \ (1 - \sigma (x))
\end{eqnarray}
●データの線形結合を計算
●シグモイド関数の入力→出力が確率に対応
●[表記]i番目データを与えた時のシグモイド関数の出力をi番目のデータがY=1になる確率とする
$$ P(Y = 1 |\ \boldsymbol{x}) = \sigma (w_0 + w_1 x_1 + \cdots + w_m x_m) $$
求めたい値P(Y = 1 | x)は説明変数の実現値が与えられた際にY=１になる確率
表記
$$ p_i = \sigma (w_0 + w_1 x_{i1} + \cdots + w_m x_{im}) $$
#3.最尤推定
$$ P(y) = P^y (1 - p)^{1-y} $$
・n回の試行でy1~ynが同時に起こる確率(p固定)
$$ P(y_1,y_2,\cdots,y_n;p) = \prod_{i=1}^{n} \color{blue}{p} ^ {y_i} (1-
\color{blue}{p})^{1-{y_i}} $$
・y1~ynのデータが得られた際の尤度関数
$$ P(y_1,y_2,\cdots,y_n;p) = \prod_{i=1}^{n} \color{red}{p} ^
\color{blue}{y_i} (1 - \color{red}{p})^{1 - \color{blue}{y_i}} $$
青字は既知のパラメータ、赤字は未知のパラメータ
\begin{eqnarray}
P(Y= y_1 |\ \boldsymbol{x}_1) =p_1^{y_1} (1 - p_1)^{1-y_1} &= \sigma(w^T
\boldsymbol{x}_1)^{y_1} (1 - \sigma (w^T \boldsymbol{x}_1))^{1-y_1} \\
P(Y= y_2 |\ \boldsymbol{x}_2) =p_2^{y_2} (1 - p_2)^{1-y_2} &= \sigma(w^T
\boldsymbol{x}_2)^{y_2} (1 - \sigma (w^T \boldsymbol{x}_2))^{1-y_2} \\
\vdots \\
P(Y= y_n |\ \boldsymbol{x}_n) =p_n^{y_n} (1 - p_n)^{1-y_n} &= \sigma(w^T
\boldsymbol{x}_n)^{y_n} (1 - \sigma (w^T \boldsymbol{x}_n))^{1-y_n} \\
\end{eqnarray}
\begin{eqnarray}
P(\color{blue}{y_1,y_2,\cdots,y_n} | \color{red}{w_0,w_1,\cdots,w_m}) &=&
\prod_{i=1}^{n} \color{red}{p_i}^{y_i} (1 - \color{red}{p_i})^{1 - {y_i}} \\
&=& \prod_{i=1}^{n} \sigma (\color{red}{w}^T \color{blue}
{\boldsymbol{x}_i})^\color{blue}{y_i} (1 - \sigma (\color{red}
{w}^T \color{blue}{\boldsymbol{x}_i}))^{1-\color{blue}{y_i}} \\
&=& L( \color{red}{w} )
\end{eqnarray}
青字は既知のパラメータ、赤字は未知のパラメータ
●尤度関数を最大とするパラメータを探索(推定)
○ 対数を取ると微分の計算が簡単
■ 同時確率の積が和に変換可能
■ 指数が積の演算に変換可能
○ 対数尤度関数が最大になる点と尤度関数が最大になる点は同じ
■ 対数関数は単調増加（ある尤度の値がx1 < x2の時、必ずlog(x1) < log(x2)となる)
○ 「尤度関数にマイナスをかけたものを最小化」し、「最小2乗法の最小化」と合わせる
\begin{eqnarray}
E(w_0,w_1,\cdots,w_m) &=& - Log(w_0,w_1,\cdots,w_m) \\
&=& - \sum_{i=1}^{n} \big\{y_i \ log p_i + (1 - y_i) \ log (1-p_i) \big\}
\end{eqnarray}
○ 反復学習によりパラメータを逐次的に更新するアプローチの一つ
○ ηは学習率と呼ばれるハイパーパラメータでモデルのパラメータの収束しやすさを調整
●なぜ必要か
○[線形回帰モデル(最小2乗法)] ▶MSEのパラメータに関する微分が0になる値を解析に求めることが可能
○[ロジスティック回帰モデル(最尤法)] ▶対数尤度関数をパラメータで微分して0になる値を求める必要があるのだが、解析的にこの値を求めることは困難である。
$$ w(k+1) = w^k - \eta \frac{\partial E(w)}{\partial w} $$
●対数尤度関数を係数とバイアスに関して微分
\begin{eqnarray}
\frac{\partial E(w)}{\partial w} &=& - \sum_{i=1}^{n} \frac{\partial
E_i(w)} {\partial p_i} \frac{\partial p_i}{\partial w} \\
&=& - \sum_{i=1}^{n} \big( \frac{y_i}{p_i} - \frac{1-y_i}{1-p_i} \big) \
p_i(1 - p) \ \boldsymbol{x}_i \\
&=& - \sum_{i=1}^{n} (y_i (1 - p_i) - p_i(1 - y_i) ) \boldsymbol{x}_i \\
&=& - \sum_{i=1}^{n} (y_i - p_i) \boldsymbol{x}_i
\end{eqnarray}
1行目：連鎖率により、分解
２行目：対数尤度関数のpに関する微分
３行目：シグモイド関数の微分
４行目：式の整理
○パラメータが更新されなくなった場合、それは勾配が0になったということ。
○少なくとも反復学習で探索した範囲では最適な解がもとめられたことになる。
$$ w^{(k+1)} = w^{(k)} + \eta \sum_{i=1}^{n} (y_i - p_i) \boldsymbol{x}_i $$
○勾配降下法では、パラメータを更新するのにN個全てのデータに対する和を求める必要がある
■ nが巨大になった時にデータをオンメモリに載せる容量が足りない、計算時間が莫大になる可能性がある
■ 確率的勾配降下法を利用して解決
$$ w^{(k+1)} = w^k + \eta (y_i - p_i) \boldsymbol{x}_i $$
#4.ロジスティック回帰モデルの性能指標
| 検証用データの結果 | |||
| positive | negative | ||
| モデルの予測結果 | positive | 真陽性(True Positive) ~正しくpositive と判別した個数 |
偽陰性(False Positive) ~間違えてpositiveと判別した個数 |
| negative | 偽陽性(False Positive) ~間違えて Negativeと判別した個数 |
真陰性(True Negative) ~正しくnegativeと判別した個数 |
|
$$ \frac{TP + TN}{TP + FN + FP + TN} $$
●再現率(Recall)
○「本当にPositiveなもの」の中からPositiveと予測できる割合
(NegativeなものをPositiveとしてしまう事象については考えていない)
○「誤り(False Positive)が多少多くても抜け漏れは少ない」予測をしたい際に利用
○使用例) 病気の検診で「陽性であるものを陰性と誤診(False Negative)」としてしまうのを避けたい。
陰性を陽性であると誤診(False Positive)」とするものが少し増えたとしても再検査すればいい。
$$ \frac{TP}{TP + FN} $$
●適合率(Precision)
○モデルが「Positiveと予測」したものの中で本当にPositiveである割合
(本当にPositiveなものをNegativeとしてしまう子については考えていない)
○見逃し(False Negative)が多くてもより正確な」予測をしたい際に利用
○重要なメールをスパムメールと誤判別」されるより「スパムと予測したものが確実にスパム」である方が便利。
スパムメールを検出できなくても(False Negative)、自分でやればいい。
$$ \frac{TP}{TP + FP} $$
●F値
○理想的にはどちらも高いモデルがいいモデルだが、両者はトレードオフの関係にあり、どちらかを小さくすると、もう片方の値が大きくなってしまう。
○PrecisionとRecallの調和平均
■ RecallとPrecisionのバランスを示している。高ければ高いほどRecallとPrecisionがともに高くなる
#5.ハンズオン
●設定
○タイタニックの乗客データを利用しロジスティック回帰モデルを作成
●目標
○年齢３０歳の男性は生き残れるか
#グーグルドライブをマウント
from google.colab import drive
drive.mount('/content/drive')
#パスの追加
import sys
sys.path.append('/content/drive/MyDrive/titanic/study_ai_ml_google')
# 必要なモジュールをインストール
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.linear_model import LogisticRegression
from pandas import DataFrame
import seaborn as sns
%matplotlib inline
# タイタニックデータの取得
titanic_df = pd.read_csv('/content/drive/MyDrive/titanic/study_ai_ml_google/data/titanic_train.csv')
| PassengerId | Survived | Pclass | Name | Sex | Age | SibSp | Parch | Ticket | Fare | Cabin | Embarked | |
| 0 | 1 | 0 | 3 | Braund, Mr. Owen Harris | male | 22.0 | 1 | 0 | A/5 21171 | 7.2500 | NaN | S |
| 1 | 2 | 1 | 1 | Cumings, Mrs. John Bradley (Florence Briggs Th... | female | 38.0 | 1 | 0 | PC 17599 | 71.2833 | C85 | C |
| 2 | 3 | 1 | 3 | Heikkinen, Miss. Laina | female | 26.0 | 0 | 0 | STON/O2. 3101282 | 7.9250 | NaN | S |
| 3 | 4 | 1 | 1 | Futrelle, Mrs. Jacques Heath (Lily May Peel) | female | 35.0 | 1 | 0 | 113803 | 53.1000 | C123 | S |
| 4 | 5 | 0 | 3 | Allen, Mr. William Henry | male | 35.0 | 0 | 0 | 373450 | 8.0500 | NaN | S |
# 不要なデータの削除
titanic_df.drop(['PassengerId','Name','Ticket','Cabin'], axis=1, inplace=True)
#削除後のデータを表示
titanic_df.head()
| Survived | Pclass | Sex | Age | SibSp | Parch | Fare | Embarked | |
| 0 | 0 | 3 | male | 22.0 | 1 | 0 | 7.2500 | S |
| 1 | 1 | 1 | female | 38.0 | 1 | 0 | 71.2833 | C |
| 2 | 1 | 3 | female | 26.0 | 0 | 0 | 7.9250 | S |
| 3 | 1 | 1 | female | 35.0 | 1 | 0 | 53.1000 | S |
| 4 | 0 | 3 | male | 35.0 | 0 | 0 | 8.0500 | S |
# nullを含んでいる行を表示
titanic_df[titanic_df.isnull().any(1)].head(5)
| Survived | Pclass | Sex | Age | SibSp | Parch | Fare | Embarked | |
| 5 | 0 | 3 | male | NaN | 0 | 0 | 8.4583 | Q |
| 17 | 1 | 2 | male | NaN | 0 | 0 | 13.0000 | S |
| 19 | 1 | 3 | female | NaN | 0 | 0 | 7.2250 | C |
| 26 | 0 | 3 | male | NaN | 0 | 0 | 7.2250 | C |
| 28 | 1 | 3 | female | NaN | 0 | 0 | 7.8792 | Q |
○NaN(null)行を含んでしまうと、正しく学習できないため、前処理の段階で必ず無くしておくこと
# Ageカラムのnullを中央値で補完（新規に「AgeFill」カラムを作成し、保存）
titanic_df['AgeFill'] = titanic_df['Age'].fillna(titanic_df['Age'].mean())
#再度確認
titanic_df[titanic_df.isnull().any(1)].head(5)
| Survived | Pclass | Sex | Age | SibSp | Parch | Fare | Embarked | AgeFill | |
| 5 | 0 | 3 | male | NaN | 0 | 0 | 8.4583 | Q | 29.699118 |
| 17 | 1 | 2 | male | NaN | 0 | 0 | 13.0000 | S | 29.699118 |
| 19 | 1 | 3 | female | NaN | 0 | 0 | 7.2250 | C | 29.699118 |
| 26 | 0 | 3 | male | NaN | 0 | 0 | 7.2250 | C | 29.699118 |
| 28 | 1 | 3 | female | NaN | 0 | 0 | 7.8792 | Q | 29.699118 |
○今回は年齢の中央値でNaNを埋めた。
その他にも最小値や最大値など様々なデータを入れれる
○まずは運賃のみで、生死フラグの判別を行えるように実装する
# 運賃だけのリストを作成
data1 = titanic_df.loc[:, ["Fare"]].values
# 生死フラグのみのリストを作成
label1 = titanic_df.loc[:, ["Survived"]].values
# モデルの作成
model = LogisticRegression()
# fit関数で学習
model.fit(data1,label1)
#実行結果
#/usr/local/lib/python3.7/dist-packages/sklearn/utils/validation.py:760: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().
# y = column_or_1d(y, warn=True)
#LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,
# intercept_scaling=1, l1_ratio=None, max_iter=100,
# multi_class='auto', n_jobs=None, penalty='l2',
# random_state=None, solver='lbfgs', tol=0.0001, verbose=0,
# warm_start=False)
○線形回帰モデルの時と同じように、説明変数、目的変数を分け、モデルの生成後、fit関数で学習させる
# 学習結果の確認(運賃が９０ドルの人は生き残ったか？)
model.predict([[90]])
#実行結果
#array([1])
# 学習結果の確認2(運賃が９０ドルの人は何％の確率で生き、何％の確率で死ぬか？)
model.predict_proba([[90]])
#実行結果
#array([[0.3949907, 0.6050093]])
○今回は試しで、運賃が９０ドルの人が生き残れるかの確認を行う
■model.predict では実際に生き残れるかを確認することができる
■model.predict_proba では何％の確率で生き、何％の確率で死ぬかを確認することができる
○次に2変数を使用して、生死を判別するように作成する
# 性別カラムを数値に変換（新規「Gender」カラムを作成し保存）
titanic_df['Gender'] = titanic_df['Sex'].map({'female' : 0, 'male' : 1}).astype(int)
#中身の確認
titanic_df.head()
| Survived | Pclass | Sex | Age | SibSp | Parch | Fare | Embarked | AgeFill | Gender | |
| 0 | 0 | 3 | male | 22.0 | 1 | 0 | 7.2500 | S | 22.0 | 1 |
| 1 | 1 | 1 | female | 38.0 | 1 | 0 | 71.2833 | C | 38.0 | 0 |
| 2 | 1 | 3 | female | 26.0 | 0 | 0 | 7.9250 | S | 26.0 | 0 |
| 3 | 1 | 1 | female | 35.0 | 1 | 0 | 53.1000 | S | 35.0 | 0 |
| 4 | 0 | 3 | male | 35.0 | 0 | 0 | 8.0500 | S | 35.0 | 1 |
#Pclass(階級)とGender(性別)の特徴量を足し合わせる
titanic_df["Pclass_Gender"] = titanic_df['Pclass'] + titanic_df['Gender']
#中身の確認
titanic_df.head()
| Survived | Pclass | Sex | Age | SibSp | Parch | Fare | Embarked | AgeFill | Gender | Pclass_Gender | |
| 0 | 0 | 3 | male | 22.0 | 1 | 0 | 7.2500 | S | 22.0 | 1 | 4 |
| 1 | 1 | 1 | female | 38.0 | 1 | 0 | 71.2833 | C | 38.0 | 0 | 1 |
| 2 | 1 | 3 | female | 26.0 | 0 | 0 | 7.9250 | S | 26.0 | 0 | 3 |
| 3 | 1 | 1 | female | 35.0 | 1 | 0 | 53.1000 | S | 35.0 | 0 | 1 |
| 4 | 0 | 3 | male | 35.0 | 0 | 0 | 8.0500 | S | 35.0 | 1 | 4 |
○ここでは新たな特徴量を作成している。
階級と性別の数字を足し合わせているが、
どのような考えで新たな特徴量作成に至ったかというと、
身近でよくある、階級が高い人を優先的に救助し、
さらにはレディファーストな考えで、女性を先に救助するのように
考えることができるのではないかと想定できる。
#不要なデータの削除
titanic_df = titanic_df.drop(['Pclass', 'Sex', 'Age', 'Gender'], axis=1)
#中身の確認
titanic_df.head()
| Survived | SibSp | Parch | Fare | Embarked | AgeFill | Pclass_Gender | |
| 0 | 0 | 1 | 0 | 7.2500 | S | 22.0 | 4 |
| 1 | 1 | 1 | 0 | 71.2833 | C | 38.0 | 1 |
| 2 | 1 | 0 | 0 | 7.9250 | S | 26.0 | 3 |
| 3 | 1 | 1 | 0 | 53.1000 | S | 35.0 | 1 |
| 4 | 0 | 0 | 0 | 8.0500 | S | 35.0 | 4 |
# 作成した特徴量の確認のため図をプロット
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, 'AgeFill'],
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, 'AgeFill'],
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('AgeFill')
ax.set_ylabel('Pclass_Gender')
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
ax.legend(bbox_to_anchor=(1.4, 1.03) )
○x軸には、年齢 y軸には新しく作成した特徴量でプロットを行った結果次のことが分かる。
■Pclass_Genderが低い(階級が高く、女性の可能性が高い）と生き残る人が多くなる
■全体的に年齢が低い（若い）方が生き残る人が多くなる
○上記の整理したデータでロジスティック回帰モデルを作成してみる
# 年齢と作成した特徴量(Pclass_Gender)でリストを作成
data2 = titanic_df.loc[:, ['AgeFill' , 'Pclass_Gender']].values
# 生死フラグのみのリストを作成
label2 = titanic_df.loc[:, 'Survived'].values
# モデルの作成と学習
model2 = LogisticRegression()
model2.fit(data2,label2)
# 学習結果の確認(年齢が２０歳で、特徴量が２の人)
model2.predict_proba([[20,2]])
# 実行結果
array([[0.21920089, 0.78079911]])
○試しに年齢が２０歳で、特徴量が２の人で予測してみると
７８％で生き残り、２１％で死亡してしまうという結果になった。
○この学習結果を先程の図に追加してみると
# 学習結果の確認（図）
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, 'AgeFill'],
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, 'AgeFill'],
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('AgeFill')
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--')
○境界線を引くことができる
この境界線より下に行った場合は生存確率が高くなり、
上に行った場合は生存確率が低くなっていることが分かる