Logistic回帰 – TauStation

1 概要
2 Logistic回帰の定式化
- 2.1 確率分布関数の仮定
- 2.2 最尤推定
3 1変数の場合
4 2変数の場合
5 特徴量の係数

概要

適用対象となる問題

ロジスティック回帰(Logistic regression)は「回帰」という名称だが、その機能は2クラス分類である。複数の特徴量に対して、2つのクラス1/0の何れになるかを判定する。

たとえば旅行会社の会員顧客20人の年齢と、温泉(SPA)とレジャーランド(LSR)のどちらを選んだかというデータが以下のように得られているとき、新たな顧客にどちらを勧めればよいか、といったような問題。

上のデータの”is SPA”列は、温泉を選んだ場合に1、レジャーランドを選んだ場合に0としている。これを散布図にすると以下の通り。

年齢が高いほど温泉を選ぶ傾向があるのは分かるが、「年齢が与えられたときに温泉とレジャーランドのどちらを選ぶか」というモデルを導出するのが目的。

(1) $\begin{equation*} (x_1, \ldots, x_m) \rightarrow \left \{ \begin{array}{cc} 1 \\ 0 \end{array} \right. \end{equation*}$

モデルを線形モデルとし、特徴量の線形和の値によってクラス1/0の何れになるかが決まるとすると

(2) $\begin{equation*} b + w_1 x_1 + \cdots + w_m x_m \rightarrow \left \{ \begin{array}{cc} 1 \\ 0 \end{array} \right. \end{equation*}$

線形回帰による場合

ここでそのままデータに対して線形回帰を適用して、たとえば回帰式の値0.5に対する年齢を境にして温泉とレジャーランドを判定してもよさそうに見える。ただ、この時の回帰式がデータに対してどのような意味を持つのかが定かではない。

import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression

path = os.path.join(os.path.dirname(__file__), "data_spa.csv")
df = pd.read_csv(path)

X = np.array(df['age']).reshape(-1, 1)
y = df['is_spa']

linreg = LinearRegression()
linreg.fit(X, y)
b = linreg.intercept_
w = linreg.coef_[0]
x_border = (0.5 - b) / w
print(x_border)

xl = 0
xr = 100
yl = b + w * xl
yr = b + w * xr

fig, ax = plt.subplots()

ax.scatter(X, y, c='tab:blue')
ax.plot([xl, xr], [yl, yr], c='tab:orange')
ax.plot([x_border, x_border], [0, 1], c='k', linestyle='dashed')

ax.set_xlim(xl, xr)
ax.set_ylim(-0.1, 1.1)
ax.set_yticks(np.linspace(0, 1, 11))
ax.grid(True)

plt.show()

import os

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

from sklearn.linear_model import LinearRegression

path = os.path.join(os.path.dirname(__file__), "data_spa.csv")

df = pd.read_csv(path)

X = np.array(df['age']).reshape(-1, 1)

y = df['is_spa']

linreg = LinearRegression()

linreg.fit(X, y)

b = linreg.intercept_

w = linreg.coef_[0]

x_border = (0.5 - b) / w

print(x_border)

xl = 0

xr = 100

yl = b + w * xl

yr = b + w * xr

fig, ax = plt.subplots()

ax.scatter(X, y, c='tab:blue')

ax.plot([xl, xr], [yl, yr], c='tab:orange')

ax.plot([x_border, x_border], [0, 1], c='k', linestyle='dashed')

ax.set_xlim(xl, xr)

ax.set_ylim(-0.1, 1.1)

ax.set_yticks(np.linspace(0, 1, 11))

ax.grid(True)

plt.show()

確率分布の仮定

単純な線形回帰ではなく、特徴量の線形和に対してクラス1となる確率分布を考える。

(3) $\begin{align*} \Pr (y=1) &= f(z) = f( b + w_1 x_1 + \cdots + w_m x_m ) \\ f(0) &= 0.5 \end{align*}$

このときの確率分布の関数形として以下が必要となる。

xの(−∞, +∞)の範囲に対して、値が(0, 1)となる
xが小さい→確率が0に近くなる→クラス0と判別されやすくなる
xが大きい→確率が1に近くなる→クラス1と判別されやすくなる

これをグラフにすると、以下のような形になる。

このような関数の回帰によるクラス分類は、確率分布に用いる関数形からLogistic回帰と呼ばれている。

Logistic回帰の定式化

確率分布関数の仮定

説明変数が $x_1, \ldots, x_m$ であり、それに対してある事象が発生する／発生しないの2通りがあるとする。

いま、 $i=1, \ldots, n$ の観察対象について $x_{i1}, \ldots, x_{im}$ のデータが観測されたとする。このとき、以下のように定式化する。

(4) $\begin{align*} & z_i = b + w_1 x_{i1} + \cdots + w_m x_{im} = b + {\boldsymbol a}{\boldsymbol x} \\ & y_i = \left\{ \begin{array}{cl} 1 & \textrm{(occured)}\\ 0 & \textrm{(not occured)} \end{array} \right. \\ & \Pr (y=1) = p = \frac{1}{1 + e^{- ( b + {\boldsymbol w}{\boldsymbol x} ) }} \end{align*}$

$x_{ij} \; (j=1, \ldots, m)$ が得られたとき、その多項式を使った確率によって事象が発生することを意味している。この確率分布の式はLogistic関数と呼ばれる。

ロジスティック関数(logistic function)はシグモイド関数(sigmoid function)とも呼ばれ、その値は(0, 1)の間にあることから確率分布として採用している。

(5) $\begin{equation*} {\rm logistic} (t) = \sigma (t) = \frac{1}{1 + e^{-t}} = p \end{equation*}$

なお、ロジスティック関数の逆関数はロジット関数(logit function)と呼ばれる。

(6) $\begin{equation*} {\rm logit} (p) = \ln \frac{p}{1 - p} = t \end{equation*}$

ここでLogit関数の対数の中はオッズの形になっている。

最尤推定

$x_{ij} (i=1, \ldots , n)$ のm×n個の説明変数データと、 $y_i (1 \; \textrm{or} \; 0)$ のターゲットデータが得られたとする。

このとき、ターゲットデータが得られたパターンとなる確率は、それぞれが独立事象であるとすれば確率の積になるから、

(7) $\begin{equation*} L &= \prod_{i=1}^n \Pr (y_i = 1) ^{y_i} \Pr (y_i = 0) ^{1 - y_i} \end{equation*}$

この確率は、パラメーター $b , w_1, \ldots, w_m$ に対する尤度関数であり、その対数尤度関数は以下のようになる。

(8) $\begin{align*} \ln L &= \sum_{i=1}^n \left( y_i \ln \Pr(y_i = 1) + (1 - y_i) \ln \Pr(y_i = 0) \right) \\ &= \sum_{i=1}^n \left( y_i \ln \Pr(y_i = 1) + (1 - y_i) \ln \left(1 - \Pr(y_i = 1) \right) \right) \\ \end{align*}$

ここで、

(9) $\begin{align*} \ln \left( 1 - \Pr(y_i = 1) \right) &= \ln \left(1 - \frac{1}{1 - e^{-(b + {\boldsymbol w}{\boldsymbol x})}} \right) \\ &= \ln \frac{e^{-(b + {\boldsymbol w}{\boldsymbol x})}}{1 + e^{-(b + {\boldsymbol a}{\boldsymbol x})}} \\ &= \ln e^{-(b + {\boldsymbol w}{\boldsymbol x})} + \ln \Pr (y_i = 1 ) \\ &= -b - {\boldsymbol w}{\boldsymbol x} + \ln \Pr (y_i = 1 ) \end{align*}$

これを使って、対数尤度関数は以下のように変形される。

(10) $\begin{align*} \ln L &= \sum_{i=1}^n \left( y_i \ln \Pr(y_i = 1) + (1 - y_i) \left( -b - {\boldsymbol w}{\boldsymbol x} + \ln \Pr (y_i = 1 ) \right) \right) \\ &= \sum_{i=1}^n \left( \ln \Pr (y_i = 1 ) + (1 - y_i) (-b - {\boldsymbol w}{\boldsymbol x}) \right) \\ &= \sum_{i=1}^n \left( \ln \frac{1}{1 + e^{-(b + {\boldsymbol w}{\boldsymbol x})}} + (1 - y_i) (-b - {\boldsymbol w}{\boldsymbol x}) \right) \\ &= \sum_{i=1}^n \left( -\ln \left( 1 + e^{-(b + {\boldsymbol w}{\boldsymbol x})} \right) + (1 - y_i) (-b - {\boldsymbol w}{\boldsymbol x}) \right) \\ \end{align*}$

得られたデータに対してこの対数尤度を最大とするような ${\boldsymbol w}$ を求めるのには、通常数値計算が用いられる。

1変数の場合

定式化

線形表現は以下のようになる。

(11) $\begin{equation*} z = b + w x \end{equation*}$

ロジスティック関数による確率表現は

(12) $\begin{equation*} \Pr (y=1) = p = \frac{1}{1 + e^{-(b + w x)}} \end{equation*}$

nこのデータセット $(x_i, y_i)$ が与えられたとき、この確率分布に関する尤度関数は

(13) $\begin{equation*} L(b, w) = \prod_{i=1}^n \left( \frac{1}{1 + e^{-(b + wx_i)}} \right)^{y_i} + \left( 1 - \frac{1}{1 + e^{-(b + wx_i)}} \right)^{1 - y_i} \end{equation*}$

その対数尤度関数は

(14) $\begin{equation*} \ln L(b, w) &=& \sum_{i=1}^{n} \left( y_i \ln \frac{1}{1 + e^{-b - w x_i}} + (1 - y_i) \ln \left( 1 - \frac{1}{1 + e^{-b - w x_i}} \right) \right) \end{equation*}$

ここでデータセット(x_i, y_i)として、が与えられたとき、対数尤度を最大とするようなb, wを計算する。

LogisticRegressionモデル

温泉のデータセットに対して、scikit-learnのLogisticRegressionモデルを適用した結果が冒頭のグラフで、コードは以下の通り。

LogisticRegressionのコンストラクターの引数としてsolver='lbfgs'を指定しているが、2020年4月時点ではこれを指定しないと警告が出るため。収束アルゴリズムにはいくつかの種類があって、データサイズや複数クラス分類・1対他分類の別などによって使い分けるようだ。

import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LogisticRegression

path = os.path.join(os.path.dirname(__file__), "data_spa.csv")
df = pd.read_csv(path)

X = np.array(df['age']).reshape(-1, 1)
y = df['is_spa']

logreg = LogisticRegression(C=1e5, solver='lbfgs')
logreg.fit(X, y)
b = logreg.intercept_[0]
w = logreg.coef_[0][0]
x_border = 0.0 - b / w
print("intercept  : {}".format(b))
print("coefficient: {}".format(w))
print("x_border   : {}".format(x_border))

xl = 0
xr = 100

x_graph = np.linspace(xl, xr)
y_graph = 1 / (1 + np.exp(- b - w * x_graph))

fig, ax = plt.subplots()

ax.scatter(X, y, c='tab:blue')
ax.plot(x_graph, y_graph, c='tab:orange')
ax.plot([x_border, x_border], [0, 1], c='k', linestyle='dashed')

ax.set_xlim(xl, xr)
ax.set_ylim(-0.1, 1.1)
ax.set_yticks(np.linspace(0, 1, 11))
ax.grid(True)

plt.show()

# intercept  : -13.44995797571776
# coefficient: 0.23116842117614878
# x_border   : 58.18250566961732

import os

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

from sklearn.linear_model import LogisticRegression

path = os.path.join(os.path.dirname(__file__), "data_spa.csv")

df = pd.read_csv(path)

X = np.array(df['age']).reshape(-1, 1)

y = df['is_spa']

logreg = LogisticRegression(C=1e5, solver='lbfgs')

logreg.fit(X, y)

b = logreg.intercept_[0]

w = logreg.coef_[0][0]

x_border = 0.0 - b / w

print("intercept : {}".format(b))

print("coefficient: {}".format(w))

print("x_border : {}".format(x_border))

xl = 0

xr = 100

x_graph = np.linspace(xl, xr)

y_graph = 1 / (1 + np.exp(- b - w * x_graph))

fig, ax = plt.subplots()

ax.scatter(X, y, c='tab:blue')

ax.plot(x_graph, y_graph, c='tab:orange')

ax.plot([x_border, x_border], [0, 1], c='k', linestyle='dashed')

ax.set_xlim(xl, xr)

ax.set_ylim(-0.1, 1.1)

ax.set_yticks(np.linspace(0, 1, 11))

ax.grid(True)

plt.show()

# intercept : -13.44995797571776

# coefficient: 0.23116842117614878

# x_border : 58.18250566961732

また、同じデータセットについてExcelで解いた例はこちら。

正則化

LogisticeRegressionモデルのコンストラクターの引数Cは正則化の強さを正の実数で指定する。大きな値を設定するほど正則化が弱く、小さくするほど正則化が効いてくる。

温泉のデータセットに対してCを変化させてみたところ、1000あたりから上はほとんど関数形が変わらず、ほぼ正則化がない状態。デフォルトのC=1(赤い線)は結構強い正則化のように思える。

興味深いのはこれらのグラフがほぼ1点で交わっており、その点での確率が0.6近いことである。温泉かレジャーランドかを選ぶ年齢と確率が正則化の程度に寄らず一定だということになるが、これが何を意味しているか、今のところよくわからない。

一方で、確率0.5に相当する年齢は正則化を強めていくにしたがって低くなっていて、これを判定基準とするとC=0.001ではほぼ全年齢で温泉が選択されることになってしまう。

さらにCが更に小さい値になると、確率曲線はどんどんなだらかになっていき、C=1e-6では全年齢にわたって0.5、すなわち温泉、レジャーランドのいずれを選ぶかは年齢に関わらず1/2の賭けのようになっている。

2変数の場合

2変数の例として、”O’REILLYの”Pythonではじめる機械学習”にあるforgeデータの例をトレースしたものをまとめた。

Logistic回帰～forgeデータ～Pythonではじめる機械学習より

特徴量の係数

式(4)によるとz_iの値が大きいほどy = 1となる確率が高くなり、小さいほどy = 0となる確率が高くなる。また、特徴量x_i ≥ 0とすると、係数w_iが正のときにはz_iを増加させる効果、負の時にはz_iを減少させる効果がある。y = 0, 1にターゲット0, 1が割り当てられ、それらがクラス0、クラス1を表すとすると、係数が正のときにはその特徴量の増加がクラス1を選択する方向に働き、係数が負の時にはクラス0をを選択する方向に働く。

以下は特徴量数2、クラス数2の簡単な例で、Feature-0は係数が正なのでクラス0を選択する確率を高め、Feature-1の係数は係数が負なのでクラス1を選択する確率を高める。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_blobs
from sklearn.linear_model import LogisticRegression

X, y = make_blobs(n_samples=100, n_features=2, centers=2, random_state=9)

logreg = LogisticRegression().fit(X, y)
coef = logreg.coef_.reshape(logreg.coef_.size)
f0 = np.linspace(-15, 3)
f1 = (- logreg.intercept_ - coef[0] * f0) / coef[1]
print("intercept:{}".format(logreg.intercept_))
print("coef:{}".format(coef))

fig, ax = plt.subplots()
ax.scatter(X[y==0][:, 0], X[y==0][:, 1], label="Class-0")
ax.scatter(X[y==1][:, 0], X[y==1][:, 1], label="Class-1")
ax.plot(f0, f1, c='g')
ax.set_xlim(-14, 4)
ax.set_ylim(-12, 6)
ax.set_xlabel("Feature-0")
ax.set_ylabel("Feature-1")
ax.legend()
plt.show()

# intercept:[2.11704743]
# coef:[ 0.89424625 -0.66371147]

import numpy as np

import matplotlib.pyplot as plt

from sklearn.datasets import make_blobs

from sklearn.linear_model import LogisticRegression

X, y = make_blobs(n_samples=100, n_features=2, centers=2, random_state=9)

logreg = LogisticRegression().fit(X, y)

coef = logreg.coef_.reshape(logreg.coef_.size)

f0 = np.linspace(-15, 3)

f1 = (- logreg.intercept_ - coef[0] * f0) / coef[1]

print("intercept:{}".format(logreg.intercept_))

print("coef:{}".format(coef))

fig, ax = plt.subplots()

ax.scatter(X[y==0][:, 0], X[y==0][:, 1], label="Class-0")

ax.scatter(X[y==1][:, 0], X[y==1][:, 1], label="Class-1")

ax.plot(f0, f1, c='g')

ax.set_xlim(-14, 4)

ax.set_ylim(-12, 6)

ax.set_xlabel("Feature-0")

ax.set_ylabel("Feature-1")

ax.legend()

plt.show()

# intercept:[2.11704743]

# coef:[ 0.89424625 -0.66371147]

概要

適用対象となる問題

線形回帰による場合

確率分布の仮定

Logistic回帰の定式化

確率分布関数の仮定

最尤推定

1変数の場合

定式化

LogisticRegressionモデル

正則化

2変数の場合

特徴量の係数

コメントを残す コメントをキャンセル

コメントを残すコメントをキャンセル