Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

ロジスティック回帰

モデル

目的変数y∈{0,1}y\in\{0,1\}の二値分類において、y=1y=1である確率をpp、y=0y=0である確率を1−p1-pとする。

ロジスティック回帰（logistic regression） はy=1y=1の確率を表現するモデル

P(y=1∣β)=p=11+exp⁡(−β⊤x)P(y = 1 \mid \boldsymbol{\beta}) = p = \frac{1}{1 + \exp(- \boldsymbol{\beta}^\top \boldsymbol{x} )}
(1)

である。ここで 11+exp⁡(−z)\frac{1}{1 + \exp(-z)} は (−∞,∞)(-\infty, \infty) の入力を (0,1)(0, 1)の範囲にする関数で、 （ロジスティック）シグモイド関数 という。

また、ロジスティック回帰はオッズ比の対数を線形モデルで説明するモデル

\renewcommand{\b} when command \b does not yet exist; use \newcommand

\renewcommand{\b}[1]{\boldsymbol{#1}}
\log\left( \frac{p}{1-p} \right) = \boldsymbol{\beta}^\top \boldsymbol{x}

ともいえる。

モデル表現の導出

（参考） 2つのモデル表現について
log⁡(p1−p)=β⊤x\log\left( \frac{p}{1-p} \right) = \boldsymbol{\beta}^\top \boldsymbol{x}
(3)

の両辺の指数をとると

p=exp⁡(β⊤x)1+exp⁡(β⊤x)p = \frac{ \exp( \boldsymbol{\beta}^\top \boldsymbol{x} ) }{ 1 + \exp( \boldsymbol{\beta}^\top \boldsymbol{x} ) }
(4)

となり、

exp⁡(β⊤x)1+exp⁡(β⊤x)=11+exp⁡(−β⊤x)\frac{ \exp( \boldsymbol{\beta}^\top \boldsymbol{x} ) }{ 1 + \exp( \boldsymbol{\beta}^\top \boldsymbol{x} ) } = \frac{ 1 }{ 1 + \exp( -\boldsymbol{\beta}^\top \boldsymbol{x} ) }
(5)

であるため、2つのモデル表現が同値である。

前半部の導出
Undefined control sequence: \b at position 37: …}{1-p} &= \exp(\̲b̲{\beta}^\top \b…

\begin{align}
\frac{p}{1-p} &= \exp(\b{\beta}^\top \b{x})\\
\implies
p &= \exp(\b{\beta}^\top \b{x}) (1-p)\\
  &= \exp(\b{\beta}^\top \b{x}) - p \exp(\b{\beta}^\top \b{x})
\end{align}

両辺をppで割ると

Undefined control sequence: \b at position 26: …gn}
1  &= \exp(\̲b̲{\beta}^\top \b…

\begin{align}
1  &= \exp(\b{\beta}^\top \b{x}) \frac{1}{p} - \exp(\b{\beta}^\top \b{x})\\
1 + \exp(\b{\beta}^\top \b{x}) &= \exp(\b{\beta}^\top \b{x}) \frac{1}{p}\\
\implies p &= \frac{ \exp(\b{\beta}^\top \b{x}) }{ 1 + \exp(\b{\beta}^\top \b{x}) }
\end{align}
後半部の導出
exp⁡(z)1+exp⁡(z)=11+exp⁡(−z)\frac{\exp(z)}{1 + \exp(z)} = \frac{1}{1 + \exp(-z)}
(8)

については、分子・分母にexp⁡(z)\exp(z)を掛けると

exp⁡(z)/exp⁡(z)(1+exp⁡(z))/exp⁡(z)=1(1+exp⁡(z))/exp⁡(z)\frac{\exp(z)/\exp(z)}{ (1 + \exp(z)) /\exp(z)} = \frac{1}{ (1 + \exp(z)) /\exp(z)}
(9)

であり、分母部分 (1+exp⁡(z))/exp⁡(z)(1 + \exp(z)) /\exp(z) については

1+exp⁡(z)exp⁡(z)=1exp⁡(z)+exp⁡(z)exp⁡(z)=exp⁡(−z)+1\begin{aligned} \frac{ 1 + \exp(z)}{ \exp(z) } &= \frac{ 1 }{ \exp(z) } + \frac{ \exp(z)}{ \exp(z) }\\ &= \exp(-z) + 1\\ \end{aligned}
(10)

となる。

※1exp⁡(z)=exp⁡(−z)\frac{ 1 }{ \exp(z) } = \exp(-z)については、 1exp⁡(z)\frac{ 1 }{ \exp(z) } は exp⁡(z)\exp(z) と乗じると 1 になるため、指数法則 en⋅em=enme^{n} \cdot e^{m} = e^{nm} より、 exp⁡(−z)\exp(-z) となる

誤差関数

ロジスティック回帰は、統計学的な言い方だと最尤推定法でパラメータを推定する。

機械学習的な言い方をすると交差エントロピー誤差を最小化するようにパラメータを推定する。

ベルヌーイ分布

ロジスティック回帰はP(y=1)=p,P(y=0)=1−pP(y=1) = p, P(y = 0) = 1-pのベルヌーイ分布に従う。

P(y)={pif y=11−pif y=0P(y) = \begin{cases} p & \text{if } y = 1\\ 1-p & \text{if } y = 0 \end{cases}
(11)

この確率質量関数は一括で書くと

P(Y=y)=py(1−p)1−yP(Y=y) = p^y (1 - p)^{1-y}
(12)

と書くことができる。

尤度関数

尤度関数L(θ)L(\theta)とは一般に確率(密度/質量)関数f(x∣θ)f(x| \theta)の積

L(θ)=∏i=1nf(xi∣θ)L(\theta) = \prod^n_{i=1} f(x_i| \theta)
(13)

である（独立に得られたサンプルを仮定するので単純な積が同時確率を意味する）。

そのため、ベルヌーイ分布の尤度関数は

L(p)=∏i=1npyi(1−p)1−yiL(p) = \prod^n_{i=1} p^{y_i} (1 - p)^{1- y_i}
(14)

となる。

ロジスティック回帰で使う場合、p∈[0,1]p \in [0, 1]はロジスティック回帰の予測値pi=σ(β⊤xi)p_i=\sigma(\beta^\top x_i)、y∈{0,1}y \in \{0, 1\}は実測値である。

L(β)=∏i=1npiyi(1−pi)1−yiL(\beta) = \prod^n_{i=1} p_i^{y_i} (1 - p_i)^{1- y_i}
(15)

確率の積だとすごく小さい値になって計算が大変なので、通常は対数を取った対数尤度を使う。

ℓ(β)=ln⁡L(β)=∑i=1n[yiln⁡pi+(1−yi)ln⁡(1−pi)]\ell(\beta)=\ln L(\beta)=\sum_{i=1}^n \left[y_i \ln p_i + (1-y_i) \ln (1-p_i)\right]
(16)

交差エントロピー誤差

ベルヌーイ分布の対数尤度関数の符号を負に反転させたものを交差エントロピー誤差（cross entropy loss）という。log lossやlogistic lossとも呼ばれる

L(β)=−ℓ(β)=−∑i=1n{yiln⁡pi+(1−yi)ln⁡(1−pi)}\mathcal{L}(\beta) = -\ell(\beta) = -\sum^n_{i=1} \{ y_i \ln p_i + (1 - y_i) \ln (1 - p_i) \}
(17)

単に言葉の問題だが、統計学系の分野では「対数尤度の最大化（最尤推定法）」という言い方をして、機械学習では「交差エントロピー誤差（=負の対数尤度）の最小化」とか言う。やってることは同じ。

スクラッチ実装例

# 対数尤度
import numpy as np

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

def log_likelihood(X, y, beta):
    p = sigmoid(X @ beta)
    c = 0.0001
    p = np.clip(p, a_min=c, a_max=1 - c) # log(p), log(1-p)でゼロをいれないようにする
    return np.sum(y * np.log(p) + (1 - y) * np.log(1 - p))

2次元の特徴量からなる、次のようなデータがあるとする

Source
from sklearn.datasets import make_blobs
import matplotlib.pyplot as plt

# データ生成
X, y = make_blobs(n_samples=300, n_features=2, centers=2, cluster_std=1.7, random_state=42)

# プロット
fig, ax = plt.subplots(figsize=(4, 3))
ax.scatter(X[y == 0, 0], X[y == 0, 1], s=20, label="y = 0")
ax.scatter(X[y == 1, 0], X[y == 1, 1], s=20, label="y = 1")
ax.legend()
ax.set(xlim=(X[:, 0].min()-1, X[:, 0].max()+1), ylim=(X[:, 1].min()-1, X[:, 1].max()+1))
ax.set(xlabel="x1", ylabel="x2")
plt.show()
<Figure size 400x300 with 1 Axes>

さまざまなβ1,β2\beta_1, \beta_2の値のもとでの対数尤度を計算すると次のような等高線になる

Source
# 対数尤度のplot
import matplotlib.pyplot as plt

# --- βのグリッド ---
beta_range = np.linspace(-50, 50, 100)
B1, B2 = np.meshgrid(beta_range, beta_range)

# --- 尤度を格子上で計算 ---
LL = np.zeros_like(B1)
for i in range(B1.shape[0]):
    for j in range(B1.shape[1]):
        beta = np.array([B1[i, j], B2[i, j]])
        LL[i, j] = log_likelihood(X, y, beta)

# --- contour plot ---
fig, ax = plt.subplots(figsize=(4, 4))
cs = ax.contour(B1, B2, LL, levels=10)
ax.clabel(cs)
ax.set_xlabel(r"$\beta_1$")
ax.set_ylabel(r"$\beta_2$")
ax.set_title("Log-likelihood contour")
plt.show()
/tmp/ipykernel_40358/1805774186.py:5: RuntimeWarning: overflow encountered in exp
  return 1 / (1 + np.exp(-x))
<Figure size 400x400 with 1 Axes>
from scipy.optimize import minimize

beta0 = np.zeros(X.shape[1])
res = minimize(
    fun=lambda beta: -log_likelihood(X, y, beta), # minimizeの関数なので負にする
    x0=beta0,
    method="BFGS",
    options={"gtol": 1e-6, "maxiter": 1000},
)

beta_hat = res.x
print("success:", res.success)
print("message:", res.message)
print("beta_hat:", beta_hat)
print("max loglik:", -res.fun)
success: True
message: Optimization terminated successfully.
beta_hat: [ 41.07515312 -13.2522211 ]
max loglik: -0.0300015001000042

推定値での識別境界は次のようになる

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

b1, b2 = beta_hat

# --- 散布図 ---
fig, ax = plt.subplots(figsize=(4, 3))
ax.scatter(X[y == 0, 0], X[y == 0, 1], s=20, label="y = 0")
ax.scatter(X[y == 1, 0], X[y == 1, 1], s=20, label="y = 1")
ax.legend()
ax.set(xlim=(X[:, 0].min()-1, X[:, 0].max()+1), ylim=(X[:, 1].min()-1, X[:, 1].max()+1))

# --- 識別境界 ---
x1 = np.linspace(X[:, 0].min(), X[:, 0].max(), 100)
eps = 1e-12
if abs(b2) > eps:
    # 切片なしモデルだと識別境界は b1*x1 + b2*x2 = 0 → x2 = - (b1/b2) * x1
    x2 = -(b1 / b2) * x1
    ax.plot(x1, x2, linewidth=2)
else:
    # β2 = 0 → 垂直線 β1 x1 = 0
    ax.axvline(0.0, linewidth=2)

ax.set(xlabel="x1", ylabel="x2", title="Decision boundary")
plt.show()
<Figure size 400x300 with 1 Axes>

勾配

交差エントロピーの勾配は

∇L(β)=∂L(β)∂β=∑i=1n(pi−yi)xi\nabla \mathcal{L}(\beta) = \frac{ \partial \mathcal{L}(\beta) } { \partial \beta } = \sum^n_{i=1} (p_i - y_i) x_i
(18)
導出

総和を取るまえの1レコード単位のものを使う。

\renewcommand{\s} when command \s does not yet exist; use \newcommand

\renewcommand{\s}{ \sigma(\beta^\top x) }
\ell(\beta) = y \ln \s + (1 - y) \ln (1 - \s)

\sに関する微分

Undefined control sequence: \s at position 27: …n}
\frac{d \ln \̲s̲}{d \beta}
&= \…

\begin{align}
\frac{d \ln \s}{d \beta}
&= \frac{d \ln \s}{d \s} \frac{d \s}{d (\beta^\top x)} \frac{d (\beta^\top x)}{d \beta}\\
&= \frac{1}{\s} \cdot \s (1 - \s) \cdot x\\
&= 1 - \s x
\\
\frac{d \ln (1 -\s)}{d \beta}
&= \frac{d \ln (1 -\s)}{d \s} \frac{d \s}{d (\beta^\top x)} \frac{d (\beta^\top x)}{d \beta}\\
&= -\frac{1}{1-\s} \cdot \s (1 - \s) \cdot x\\
&= - \s x
\end{align}

より、

Undefined control sequence: \s at position 54: …eta}
&= y (1 - \̲s̲)x - (1 - y) \s…

\begin{align}
\frac{d \ell(\beta)}{d\beta}
&= y (1 - \s)x - (1 - y) \s x\\
&= (y - y\s - \s + y\s) x\\
&= (y - \s ) x\\
\end{align}

（参考）使った微分

dlog⁡xdx=1xdσ(x)dx=σ(x)(1−σ(x))dσ(β⊤x)dx=dσ(β⊤x)d(β⊤x)×d(β⊤x)dx=σ(β⊤x)(1−σ(β⊤x))×x\begin{align} \frac{d \log x}{dx} &= \frac{1}{x}\\ \frac{d \sigma(x)}{dx} &= \sigma(x) (1- \sigma(x))\\ \frac{d \sigma(\beta^\top x)}{dx} &= \frac{d \sigma(\beta^\top x)}{d(\beta^\top x)} \times \frac{d (\beta^\top x)}{dx} = \sigma(\beta^\top x) (1- \sigma(\beta^\top x)) \times x \\ \end{align}
(22)
from sklearn.datasets import load_breast_cancer
X, y = load_breast_cancer(return_X_y=True)

y∈{−1,1}y \in \{-1, 1\}にする場合

上記の例ではy∈{0,1}y \in \{0, 1\}としていた。y∈{−1,1}y\in \{-1, 1\}とする場合は少し書き方が変わる

y=1y=1の確率とy=−1y=-1の確率がそれぞれ

p(y=1∣x)=exp⁡(β⊤x)1+exp⁡(β⊤x)=11+exp⁡(−β⊤x)p(y=−1∣x)=1−p(y=1∣x)=1+exp⁡(β⊤x)1+exp⁡(β⊤x)−exp⁡(β⊤x)1+exp⁡(β⊤x)=11+exp⁡(β⊤x)\begin{align} p(y=1|x) &= \frac{\exp(\beta^\top x)}{1 + \exp(\beta^\top x)} = \frac{ 1 }{1 + \exp(- \beta^\top x)}\\ p(y=-1|x) &= 1 - p(y=1|x)\\ &= \frac{ 1 + \exp(\beta^\top x)}{1 + \exp(\beta^\top x)} - \frac{\exp(\beta^\top x)}{1 + \exp(\beta^\top x)} = \frac{ 1 }{1 + \exp(\beta^\top x)}\\ \end{align}
(23)

で表されるとする。y∈{−1,1}y\in \{-1, 1\}のとき、これらを1つにまとめて、yyの確率を

p(y∣x)=11+exp⁡(−yβ⊤x)p(y|x) = \frac{1}{1 + \exp(- y \beta^\top x )}
(24)

と書くことができる。

尤度は

∏i=1n11+exp⁡(−yiβ⊤xi)\prod^n_{i=1} \frac{1}{1 + \exp(- y_i \beta^\top x_i )}
(25)

負の対数尤度は

∑i=1nlog⁡(1+exp⁡(−yiβ⊤xi))\sum^n_{i=1} \log \big( 1 + \exp(- y_i \beta^\top x_i ) \big)
(26)

と書くことができる。機械学習の分野だとこちらの表現のほうが目にするかも。

線形分離可能性

機械学習として（目的変数の予測が目的）のロジスティック回帰では、線形分離可能な問題であることが嬉しい

統計学としては最尤推定量が存在しない（解が一意に定まらない）という扱いになる様子

機械学習
線形判別モデル
機械学習
Support Vector Machine