眠れない夜にロジスティック回帰についてまとめてみた

はじめに

最近ロジスティック回帰について、少し悩むことが多いでのでちょっと自分なりにまとめてみることを試みました。

本編

ロジスティック回帰とは

ロジスティック回帰は確率値を出力する回帰器と理解しています。2値分類でよく使われるので分類器としてのイメージが強いですが、あくまで回帰器をそう使っている、というだけです。
ちなみに、式は以下です。

$y = a_1x_1 + a_2x_2 + \cdots + a_nx_n + \beta$
$p = \frac{1}{1+exp(-y)}$

式だけ見てもよくわからないので、いくつか可視化してみようと思います。

モデルの用意

今回は単純にirisデータセットにある3種類の花データから2つを取り出して2値分類モデル化します。
なお、わかりやすいようにきれいに分かれるところを使っています。

import numpy as np
import matplotlib.pyplot as plt
[f:id:marufeuillex:20191023231808p:plain]
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
import seaborn as sns
from sklearn import datasets

dataset = datasets.load_iris()
data = pd.DataFrame(dataset.data, columns=dataset.feature_names)
data["Species"] = dataset.target

train = data[(data.Species == 0) | (data.Species == 1)]

sns.scatterplot(data=train, x="petal length (cm)", y="petal width (cm)", hue="Species")

f:id:marufeuillex:20191023230738p:plain — 可視化

とりあえず可視化しました。なんかこれならif-thenで分けられそうですが、そういう話をしたいわけではないので先に進みます。

特徴量の抽出

それぞれX, Yに必要な特徴量のみ、絞り込んで投入します。非常に恣意的ですがご容赦ください。
あと、今回は特に制度の評価がしたいわけでもないので、データは全て学習用ということにしてしまいます。

X = train.drop(columns=["Species", "sepal length (cm)", "sepal width (cm)"])
Y = train["Species"]

モデルの作成

いよいよモデルの作成です。が、非常に単純です。

logi_model = LogisticRegression()
logi_model = logi_model.fit(X, Y)

これで学習まで完了します。

推定されたパラメータの確認

最初に書いたとおり、ロジスティック回帰の中では $y = a_1x_1 + a_2x_2 + \cdots + a_nx_n + \beta$ という式が使われるので、この式を確認していきます。

係数は以下で確認可能です。

logi_model.coef_ # 係数
logi_model.intercept_ # 切片

可視化1: 分類の線を引いてみる

まずは最初に書いた散布図に分類の線を引いてみましょう。
今回は $x_0, x_1$ の中に直線を引く必要があるため、 $x_0$ をx軸、 $x_1$ がy軸となるように式変形しますが、そうなるとyが邪魔です。
確率が50%を超えたらクラス1の方に所属すると考えると、先程の式のp = 0.5となりますが、それは直接 $a_1x_1 + a_2x_2 + \cdots + a_nx_n + \beta$ とリンクせず、 $p = \frac{1}{1+exp(-y)}$ を介してリンクすることに注意してください。
これを変形して得られる$y = \log\frac{p}{1-p}$という変換を行う必要があります。

その上で式を変換をした上で実装します。

まずは、式を関数として定義します。

def div_line(x):
    threshold_p = 0.5
    threshold = np.log(threshold_p+1e-10) - np.log(1-threshold_p+1e-10)
    a = logi_model.coef_[0][0]
    b = logi_model.coef_[0][1]
    c = logi_model.intercept_[0]
    return (-a * x - c + threshold)/b

実際のプロットは以下のように行います。

# プロット
setosa = X[Y == 0]
versicolor = X[Y ==1]
plt.scatter(setosa["petal length (cm)"],  setosa["petal width (cm)"], color="red", label="setosa")
plt.scatter(versicolor["petal length (cm)"],  versicolor["petal width (cm)"], color="blue", label="versicolor")

# 回帰直線の計算
min_x = np.min(X["petal length (cm)"])
max_x = np.max(X["petal length (cm)"])
min_y = np.min(X["petal width (cm)"])
max_y = np.max(X["petal width (cm)"])
x_range = np.arange(min_x, max_x, 0.1)
y = div_line(x_range)
display_range = (min_y < y) & (y < max_y)

# 回帰直線の表示
plt.plot(x_range[display_range],y[display_range], color="gray", label="divided(th=0.5)")
plt.legend()

f:id:marufeuillex:20191023231808p:plain

概ねこんな絵になったかと思います。pを変えることで直線の位置が変わるので是非試してください。

可視化2: 確率値の計算とデータの分布を見る

では次に、 $p = \frac{1}{1+exp(-y)}$ の直線を書き、実際のデータを当てはめて眺めてみようと思います。

y(x)はyを計算する関数で、inv_logitはyからpを計算する関数です。

def y(x):
    a = logi_model.coef_[0][0]
    b = logi_model.coef_[0][1]
    c = logi_model.intercept_[0]
    return a * x[:,0] + b * x[:, 1] + c

def inv_logit(x):
    return (1 /( 1 + np.exp(-x)))

そして、上記関数を利用して以下のように絵を書きます。

preds = logi_model.predict(X)
preds_proba = logi_model.predict_proba(X)[:, 1]
values = y(X.values)
setosa = Y == 0
versicolor = Y ==1
plt.scatter(values[setosa], preds[setosa], color="red", label="setosa")
plt.scatter(values[versicolor], preds[versicolor], color="blue", label="versicolor")

plt.scatter(values[setosa], preds_proba[setosa], color="orange", label="setosa(proba)")
plt.scatter(values[versicolor], preds_proba[versicolor], color="cyan", label="versicolor(proba)")
x_range = np.arange(np.min(values)-1, np.max(values)+1, 0.1)
plt.plot(x_range, inv_logit(x_range), color="gray", label="logi")

plt.legend()

f:id:marufeuillex:20191023232250p:plain

すべてのデータの計算結果yがどのようになっていて、それがどのように確率値pに紐づくかが見えると思います。
要はこの絵は横軸が線形式の計算結果で、それを確率値に変換する動きをしているということです。

可視化3: 何と何が線形なのか？

$y = a_1x_1 + a_2x_2 + \cdots + a_nx_n + \beta$ と $p = \frac{1}{1+exp(-y)}$ の日本の式をくっつけると、 $\log\frac{p}{1-p}= a_1x_1 + a_2x_2 + \cdots + \beta$ という式となります。
これは左辺が対数オッズで、右辺が線形の式です。