機械学習と情報技術

ロジスティック回帰の理論と最尤推定

多変量解析

患者の年齢・血圧・コレステロール値から、心臓病を発症するかしないか（0か1か）を予測したいとします。重回帰分析を使って予測値を計算すると、0未満や1を超える値が出てしまうかもしれません。確率として解釈するには、出力が0から1の間に収まる必要があります。

あるいは、メールがスパムかどうかを判定する問題を考えてみましょう。メールの長さ、特定のキーワードの出現回数、送信者の情報といった特徴量から「スパムである確率」を推定したいとき、通常の線形回帰では確率としての制約を満たせません。「予測確率が$-0.3$」や「確率が$1.5$」といった値は意味を持ちません。

ロジスティック回帰（Logistic Regression）は、二値の応答変数（はい/いいえ、合格/不合格、発症/非発症など）を予測するための回帰モデルです。線形回帰の出力をシグモイド関数（ロジスティック関数）に通すことで、出力を確率として解釈できる0から1の範囲に変換します。名前に「回帰」と付いていますが、実質的には分類の手法であり、機械学習における最も基本的かつ重要なアルゴリズムの一つです。

ロジスティック回帰を理解すると、以下のような応用が開けます。

  • 医療: 疾患の発症リスク予測。オッズ比によるリスク要因の定量化
  • 信用スコアリング: ローンのデフォルト確率の推定
  • マーケティング: 購入確率の予測と顧客セグメンテーション
  • 自然言語処理: テキスト分類の基本モデル
  • 機械学習: ニューラルネットワークの出力層（ソフトマックス関数）の基礎
  • 疫学: リスクファクターの効果をオッズ比として定量化する

本記事では、ロジスティック回帰モデルを対数オッズ（ロジット）の線形モデルとして定式化し、パラメータの最尤推定を導出します。対数尤度の凹性を証明し、勾配降下法とIRLSの両方を導出・実装して理論を確認します。

本記事の内容

  • ロジスティック関数（シグモイド関数）の性質と導出の動機
  • ロジットモデルと対数オッズの解釈
  • オッズ比の統計的意味
  • 最尤推定量の導出（対数尤度、スコア関数、ヘッセ行列）
  • 対数尤度の凹性の証明
  • 勾配降下法と反復重み付き最小二乗法（IRLS）
  • 多クラスへの拡張（ソフトマックス回帰）
  • Pythonでの実装と判別境界の可視化

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

シグモイド関数とロジットモデル

なぜ線形回帰では不十分か

二値の目的変数 $y \in \{0, 1\}$ に対して、直接 $y = \bm{x}^\top\bm{\beta} + \epsilon$ というモデルを当てはめると（線形確率モデルと呼ばれます）、深刻な問題が生じます。

  1. 値域の制約違反: 予測値 $\hat{y}$ が0未満や1を超えてしまい、確率としての解釈が破綻する
  2. 等分散性の違反: $y$ が0か1しか取らないため、$\text{Var}(\epsilon | \bm{x}) = p(\bm{x})(1 – p(\bm{x}))$ となり、$p(\bm{x})$ に依存する。等分散性の仮定が成り立たない
  3. 正規性の違反: 誤差項は正規分布に従わない（ベルヌーイ分布に従う）

必要なのは、線形予測子 $z = \bm{x}^\top\bm{\beta}$ を確率の範囲 $[0, 1]$ に写す関数です。このような関数はいくつか候補がありますが（プロビット関数、補対数対数関数など）、最も広く使われるのがシグモイド関数です。

シグモイド関数（ロジスティック関数）

シグモイド関数（sigmoid function）は次のように定義されます。

$$ \begin{equation} \sigma(z) = \frac{1}{1 + e^{-z}} = \frac{e^z}{1 + e^z} \end{equation} $$

この関数はS字型（シグモイド = S字の意味）のカーブを描き、以下の重要な性質を持ちます。

  • 定義域: $z \in (-\infty, \infty)$
  • 値域: $\sigma(z) \in (0, 1)$ — 常に0と1の間に収まるため、確率として解釈できます
  • 対称性: $\sigma(-z) = 1 – \sigma(z)$ — 原点に関して点対称
  • $z = 0$ での値: $\sigma(0) = 0.5$ — 判別境界に対応
  • 微分: $\sigma'(z) = \sigma(z)(1 – \sigma(z))$ — 微分が自分自身で表せる美しい性質
  • 単調増加: $\sigma'(z) > 0$ — $z$ が増えれば確率も増える

直感的には、$z$ が大きな正の値のとき確率は1に近づき、大きな負の値のとき確率は0に近づきます。$z = 0$ が確率0.5の境界（判別境界）に対応します。

微分の性質 $\sigma'(z) = \sigma(z)(1 – \sigma(z))$ は特に重要です。この性質のおかげで、最尤推定の勾配計算が非常にシンプルになります。導出は次の通りです。

$$ \sigma'(z) = \frac{e^{-z}}{(1 + e^{-z})^2} = \frac{1}{1 + e^{-z}} \cdot \frac{e^{-z}}{1 + e^{-z}} = \sigma(z) \cdot (1 – \sigma(z)) $$

ロジスティック回帰モデル

ロジスティック回帰は、$y = 1$ となる確率を次のようにモデル化します。

$$ \begin{equation} P(y = 1 | \bm{x}) = \sigma(\bm{x}^\top\bm{\beta}) = \frac{1}{1 + \exp(-\bm{x}^\top\bm{\beta})} \end{equation} $$

$p(\bm{x}) = P(y = 1 | \bm{x})$ と書くと、対数オッズ（log-odds, logit）は次のように線形になります。

$$ \begin{equation} \ln\frac{p(\bm{x})}{1 – p(\bm{x})} = \bm{x}^\top\bm{\beta} = \beta_0 + \beta_1 x_1 + \dots + \beta_p x_p \end{equation} $$

この変換をロジット変換と呼びます。ロジスティック回帰は、確率そのものではなく対数オッズが説明変数の線形関数であると仮定するモデルです。

なぜ対数オッズなのでしょうか。確率 $p$ を直接線形モデルにすると値域の制約に違反します。オッズ $p/(1-p)$ にすると正の値に制限されますが、まだ非負制約があります。対数を取ることで値域が $(-\infty, \infty)$ に拡がり、線形予測子と整合します。

オッズ比の解釈

対数オッズの線形性から、偏回帰係数 $\beta_j$ には明確な解釈があります。$x_j$ が1単位増加すると、対数オッズが $\beta_j$ だけ増加します。オッズ（$p/(1-p)$）で見ると、$\exp(\beta_j)$ 倍になります。

$$ \begin{equation} \text{Odds ratio} = e^{\beta_j} \end{equation} $$

たとえば $\beta_j = 0.5$ なら、$x_j$ が1単位増えるとオッズは $e^{0.5} \approx 1.65$ 倍になります。$\beta_j = 0$ ならオッズに影響せず（$e^0 = 1$）、$\beta_j < 0$ ならオッズが減少します（$e^{\beta_j} < 1$）。

このオッズ比の解釈は、医学統計学においてリスク要因の効果を定量化するために広く使われています。たとえば「喫煙者は非喫煙者に比べて肺がんのオッズが3.5倍」という表現は、$\exp(\beta_{\text{smoking}}) = 3.5$、つまり $\beta_{\text{smoking}} = \ln 3.5 \approx 1.25$ であることを意味します。

重要な注意点として、オッズ比はリスク比（相対リスク）とは異なります。稀な疾患ではオッズ比はリスク比によく近似しますが、発症率が高い場合にはオッズ比はリスク比よりも大きくなります。

パラメータ $\bm{\beta}$ を推定する方法を見ていきましょう。

最尤推定量の導出

尤度関数

$n$ 個の独立なデータ $(\bm{x}_i, y_i)$（$y_i \in \{0, 1\}$）に対する尤度関数を構成します。$p_i = P(y_i = 1 | \bm{x}_i) = \sigma(\bm{x}_i^\top\bm{\beta})$ とすると、$y_i$ の確率は次のように書けます。

$$ P(y_i | \bm{x}_i) = p_i^{y_i}(1 – p_i)^{1 – y_i} $$

これは $y_i = 1$ のとき $p_i$、$y_i = 0$ のとき $1 – p_i$ を返す簡潔な表現です。

独立性の仮定より、尤度関数は各データの確率の積です。

$$ \begin{equation} L(\bm{\beta}) = \prod_{i=1}^{n}p_i^{y_i}(1 – p_i)^{1 – y_i} \end{equation} $$

対数尤度関数

対数を取ると積が和に変わり、扱いやすくなります。

$$ \begin{equation} \ell(\bm{\beta}) = \sum_{i=1}^{n}\left[y_i\ln p_i + (1 – y_i)\ln(1 – p_i)\right] \end{equation} $$

$p_i = \sigma(\bm{x}_i^\top\bm{\beta})$ を代入し、シグモイド関数の性質を使って整理します。$\ln\sigma(z) = z – \ln(1 + e^z)$ と $\ln(1 – \sigma(z)) = -\ln(1 + e^z)$ を使うと次のようになります。

$$ \ell(\bm{\beta}) = \sum_{i=1}^{n}\left[y_i\bm{x}_i^\top\bm{\beta} – \ln(1 + \exp(\bm{x}_i^\top\bm{\beta}))\right] $$

この対数尤度関数は交差エントロピー損失（cross-entropy loss）の負に等しく、機械学習ではこの負の対数尤度を最小化する問題として定式化されます。

勾配（スコア関数）の導出

対数尤度を $\bm{\beta}$ で微分して勾配を求めます。$z_i = \bm{x}_i^\top\bm{\beta}$ とおくと、各項の微分は次のようになります。

$$ \frac{\partial}{\partial\bm{\beta}}\left[y_i z_i – \ln(1 + e^{z_i})\right] = y_i\bm{x}_i – \frac{e^{z_i}}{1 + e^{z_i}}\bm{x}_i = (y_i – \sigma(z_i))\bm{x}_i = (y_i – p_i)\bm{x}_i $$

全データの和を取ると、勾配ベクトル（スコア関数）が得られます。

$$ \begin{equation} \frac{\partial \ell}{\partial \bm{\beta}} = \sum_{i=1}^{n}(y_i – p_i)\bm{x}_i = \bm{X}^\top(\bm{y} – \bm{p}) \end{equation} $$

ここで $\bm{p} = (p_1, p_2, \dots, p_n)^\top$、$p_i = \sigma(\bm{x}_i^\top\bm{\beta})$ です。

この勾配の形は驚くほどシンプルです。残差 $(y_i – p_i)$ を重みとして説明変数 $\bm{x}_i$ を加重平均しています。重回帰分析の正規方程式 $\bm{X}^\top(\bm{y} – \bm{X}\bm{\beta}) = \bm{0}$ と比較すると、$\bm{X}\bm{\beta}$ が $\bm{p}$（非線形な予測確率）に置き換わっただけです。

勾配をゼロとおいた方程式 $\bm{X}^\top(\bm{y} – \bm{p}) = \bm{0}$ は $\bm{\beta}$ について非線形であるため、解析的な閉形式解は存在しません。これは重回帰分析との大きな違いです。反復的な数値最適化が必要です。

ヘッセ行列の導出

ニュートン法やIRLSに必要なヘッセ行列（対数尤度の2次微分）を計算します。

$p_i$ の $\bm{\beta}$ に関する微分は次の通りです。

$$ \frac{\partial p_i}{\partial\bm{\beta}} = \sigma'(z_i)\bm{x}_i = p_i(1-p_i)\bm{x}_i $$

したがって勾配の各成分の微分は次のようになります。

$$ \frac{\partial^2\ell}{\partial\bm{\beta}\partial\bm{\beta}^\top} = -\sum_{i=1}^{n}p_i(1-p_i)\bm{x}_i\bm{x}_i^\top = -\bm{X}^\top\bm{W}\bm{X} $$

ここで $\bm{W} = \text{diag}(p_1(1-p_1), p_2(1-p_2), \dots, p_n(1-p_n))$ は対角重み行列です。

対数尤度の凹性の証明

ヘッセ行列が常に半負定値であることを示します。任意の非ゼロベクトル $\bm{v}$ に対して次が成り立ちます。

$$ \bm{v}^\top\left(-\bm{X}^\top\bm{W}\bm{X}\right)\bm{v} = -\sum_{i=1}^{n}p_i(1-p_i)(\bm{x}_i^\top\bm{v})^2 $$

$0 < p_i < 1$ であるため $p_i(1-p_i) > 0$ が成り立ち、$(\bm{x}_i^\top\bm{v})^2 \geq 0$ です。したがって右辺は常に非正であり、ヘッセ行列は半負定値です。

これは対数尤度が凹関数であることを意味します。凹関数では局所最大が大域最大でもあるため、任意の初期値から出発しても勾配法やニュートン法は確実に大域最適解に収束します。これはロジスティック回帰の非常に好ましい性質です。

ニュートン=ラフソン法（IRLS）

ニュートン法の更新式は次のようになります。

$$ \bm{\beta}^{(t+1)} = \bm{\beta}^{(t)} – \left(\frac{\partial^2\ell}{\partial\bm{\beta}\partial\bm{\beta}^\top}\right)^{-1}\frac{\partial\ell}{\partial\bm{\beta}} $$

ヘッセ行列と勾配を代入すると次のようになります。

$$ \begin{equation} \bm{\beta}^{(t+1)} = \bm{\beta}^{(t)} + (\bm{X}^\top\bm{W}^{(t)}\bm{X})^{-1}\bm{X}^\top(\bm{y} – \bm{p}^{(t)}) \end{equation} $$

この式を整理すると、反復重み付き最小二乗法（IRLS）の形になります。

$$ \bm{\beta}^{(t+1)} = (\bm{X}^\top\bm{W}^{(t)}\bm{X})^{-1}\bm{X}^\top\bm{W}^{(t)}\bm{z}^{(t)} $$

ここで $\bm{z}^{(t)} = \bm{X}\bm{\beta}^{(t)} + (\bm{W}^{(t)})^{-1}(\bm{y} – \bm{p}^{(t)})$ は作業応答変量（working response）です。各反復で、重み $\bm{W}$ で重み付けした $\bm{z}$ に対する重回帰分析を行う、というのがIRLSの直感的な解釈です。

IRLSはニュートン法の一種であるため二次収束（各反復で有効桁数がおよそ倍増）し、通常は5〜10回の反復で十分な精度に収束します。

勾配降下法との比較

勾配降下法は1次の情報（勾配）のみを使う方法で、更新式は次の通りです。

$$ \bm{\beta}^{(t+1)} = \bm{\beta}^{(t)} + \alpha\bm{X}^\top(\bm{y} – \bm{p}^{(t)}) $$

ここで $\alpha > 0$ は学習率です。勾配降下法はIRLSに比べて各反復の計算コストは低いですが、収束は線形（1次）であり、より多くの反復が必要です。大規模データでは確率的勾配降下法（SGD）が使われ、全データの一部（ミニバッチ）で勾配を近似することで計算を高速化します。

多クラスへの拡張 — ソフトマックス回帰

ロジスティック回帰は2クラスの分類問題を扱いますが、$K > 2$ クラスに拡張した多項ロジスティック回帰（ソフトマックス回帰）も広く使われています。

クラス $k$ に属する確率は次のように定義されます。

$$ \begin{equation} P(y = k | \bm{x}) = \frac{\exp(\bm{x}^\top\bm{\beta}_k)}{\sum_{j=1}^{K}\exp(\bm{x}^\top\bm{\beta}_j)} \end{equation} $$

この関数はソフトマックス関数と呼ばれ、ニューラルネットワークの出力層で標準的に使用されています。$K = 2$ の場合、ソフトマックス回帰は通常のロジスティック回帰に帰着します。

理論をPythonで実装して確認しましょう。

Pythonによる実装と可視化

勾配降下法によるロジスティック回帰

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

# --- 合成データの生成 ---
n = 300
# クラス0
X0 = np.random.multivariate_normal([1, 1], [[1, 0.3], [0.3, 1]], n // 2)
# クラス1
X1 = np.random.multivariate_normal([3, 3], [[1, 0.3], [0.3, 1]], n // 2)
X_raw = np.vstack([X0, X1])
y = np.array([0] * (n // 2) + [1] * (n // 2))

# 切片を追加
X = np.column_stack([np.ones(n), X_raw])

# --- シグモイド関数 ---
def sigmoid(z):
    return 1.0 / (1.0 + np.exp(-np.clip(z, -500, 500)))

# --- 対数尤度 ---
def log_likelihood(X, y, beta):
    z = X @ beta
    return np.sum(y * z - np.log(1 + np.exp(np.clip(z, -500, 500))))

# --- 勾配降下法 ---
beta = np.zeros(X.shape[1])
lr = 0.01
n_iter = 500
ll_history = []

for iteration in range(n_iter):
    p = sigmoid(X @ beta)
    gradient = X.T @ (y - p)
    beta += lr * gradient
    ll_history.append(log_likelihood(X, y, beta))

# --- IRLS（ニュートン法） ---
beta_irls = np.zeros(X.shape[1])
for iteration in range(20):
    p = sigmoid(X @ beta_irls)
    W = np.diag(p * (1 - p))
    H = X.T @ W @ X
    grad = X.T @ (y - p)
    beta_irls += np.linalg.solve(H, grad)

print("=== ロジスティック回帰の結果 ===")
print(f"勾配降下法: beta = {beta}")
print(f"IRLS:       beta = {beta_irls}")
print(f"差: {np.max(np.abs(beta - beta_irls)):.6f}")

# --- 予測精度 ---
p_pred = sigmoid(X @ beta)
y_pred = (p_pred >= 0.5).astype(int)
accuracy = np.mean(y_pred == y)
print(f"正解率: {accuracy:.4f}")

勾配降下法とIRLSの両方の結果が近い値に収束していることが確認できます。IRLSは二次の情報（ヘッセ行列）を使うため、はるかに少ない反復回数で高精度に収束します。勾配降下法は500回の反復が必要でしたが、IRLSは通常10回以内で十分です。

判別境界と確率の可視化

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)
n = 300
X0 = np.random.multivariate_normal([1, 1], [[1, 0.3], [0.3, 1]], n // 2)
X1 = np.random.multivariate_normal([3, 3], [[1, 0.3], [0.3, 1]], n // 2)
X_raw = np.vstack([X0, X1])
y = np.array([0]*(n//2) + [1]*(n//2))
X = np.column_stack([np.ones(n), X_raw])

def sigmoid(z):
    return 1.0 / (1.0 + np.exp(-np.clip(z, -500, 500)))

# IRLS
beta = np.zeros(X.shape[1])
for _ in range(20):
    p = sigmoid(X @ beta)
    W = np.diag(p * (1 - p))
    beta += np.linalg.solve(X.T @ W @ X, X.T @ (y - p))

fig, axes = plt.subplots(1, 3, figsize=(16, 5))

# (a) 判別境界
ax = axes[0]
xx, yy = np.meshgrid(np.linspace(-2, 6, 200), np.linspace(-2, 6, 200))
grid = np.column_stack([np.ones(xx.size), xx.ravel(), yy.ravel()])
Z = sigmoid(grid @ beta).reshape(xx.shape)
ax.contourf(xx, yy, Z, levels=20, cmap="RdBu_r", alpha=0.4)
ax.contour(xx, yy, Z, levels=[0.5], colors="k", linewidths=2)
ax.scatter(X0[:, 0], X0[:, 1], c="blue", alpha=0.4, s=15, label="Class 0")
ax.scatter(X1[:, 0], X1[:, 1], c="red", alpha=0.4, s=15, label="Class 1")
ax.set_xlabel("$x_1$", fontsize=11)
ax.set_ylabel("$x_2$", fontsize=11)
ax.set_title("(a) Decision boundary (P=0.5)", fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# (b) シグモイド関数
ax = axes[1]
z = np.linspace(-6, 6, 200)
ax.plot(z, sigmoid(z), "b-", linewidth=2)
ax.axhline(0.5, color="gray", linestyle="--", linewidth=0.8)
ax.axvline(0, color="gray", linestyle="--", linewidth=0.8)
ax.fill_between(z, sigmoid(z), 0, where=(z > 0), alpha=0.1, color="red")
ax.fill_between(z, sigmoid(z), 0, where=(z < 0), alpha=0.1, color="blue")
ax.set_xlabel("$z = \\mathbf{x}^T\\boldsymbol{\\beta}$", fontsize=11)
ax.set_ylabel("$\\sigma(z) = P(y=1|\\mathbf{x})$", fontsize=11)
ax.set_title("(b) Sigmoid function", fontsize=12)
ax.grid(True, alpha=0.3)

# (c) 対数尤度の収束
beta_gd = np.zeros(X.shape[1])
lr = 0.01; ll_hist = []
for _ in range(500):
    p_val = sigmoid(X @ beta_gd)
    beta_gd += lr * X.T @ (y - p_val)
    ll = np.sum(y * np.log(p_val + 1e-15) + (1-y) * np.log(1 - p_val + 1e-15))
    ll_hist.append(ll)
ax = axes[2]
ax.plot(ll_hist, "b-", linewidth=1.5)
ax.set_xlabel("Iteration", fontsize=11)
ax.set_ylabel("Log-likelihood", fontsize=11)
ax.set_title("(c) Convergence of gradient descent", fontsize=12)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig("logistic_regression.png", dpi=150, bbox_inches="tight")
plt.show()

この可視化から、以下の知見が得られます。

  1. 判別境界（図a）: ロジスティック回帰の判別境界は直線（$\bm{x}^\top\bm{\beta} = 0$ の超平面）です。背景色は予測確率を表し、赤い領域ほどクラス1の確率が高く、青い領域ほどクラス0の確率が高くなっています。境界から離れるほど確率が0か1に近づくグラデーションが見られ、これはシグモイド関数の飽和特性を反映しています

  2. シグモイド関数（図b）: S字型のカーブが、線形予測子 $z$ を確率に変換する仕組みを示しています。$z = 0$ で確率0.5（判別境界）、$z > 0$ で0.5以上（クラス1に分類）、$z < 0$ で0.5未満（クラス0に分類）です。勾配は $z = 0$ の近くで最大であり、判別境界付近でのデータの変化がモデルに最も大きな影響を与えます

  3. 対数尤度の収束（図c）: 勾配降下法による対数尤度が単調に増加し、収束していることが確認できます。初期の急激な上昇の後、次第に平坦になっていく形は、凹関数の最適化の典型的なパターンです

IRLSと勾配降下法の収束速度の比較

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)
n = 300
X0 = np.random.multivariate_normal([1, 1], [[1, 0.3], [0.3, 1]], n // 2)
X1 = np.random.multivariate_normal([3, 3], [[1, 0.3], [0.3, 1]], n // 2)
X_raw = np.vstack([X0, X1])
y = np.array([0]*(n//2) + [1]*(n//2))
X = np.column_stack([np.ones(n), X_raw])

def sigmoid(z):
    return 1.0 / (1.0 + np.exp(-np.clip(z, -500, 500)))

# 真の最適解をIRLSで事前に求める
beta_opt = np.zeros(X.shape[1])
for _ in range(50):
    p = sigmoid(X @ beta_opt)
    W = np.diag(p * (1 - p))
    beta_opt += np.linalg.solve(X.T @ W @ X, X.T @ (y - p))

# 勾配降下法
beta_gd = np.zeros(X.shape[1])
gd_errors = []
for _ in range(100):
    p = sigmoid(X @ beta_gd)
    beta_gd += 0.01 * X.T @ (y - p)
    gd_errors.append(np.linalg.norm(beta_gd - beta_opt))

# IRLS
beta_irls = np.zeros(X.shape[1])
irls_errors = []
for _ in range(10):
    p = sigmoid(X @ beta_irls)
    W = np.diag(p * (1 - p) + 1e-10)
    beta_irls += np.linalg.solve(X.T @ W @ X, X.T @ (y - p))
    irls_errors.append(np.linalg.norm(beta_irls - beta_opt))

fig, ax = plt.subplots(figsize=(8, 5))
ax.semilogy(range(1, len(gd_errors)+1), gd_errors, "b-", linewidth=1.5,
            label="Gradient descent ($\\alpha=0.01$)")
ax.semilogy(range(1, len(irls_errors)+1), irls_errors, "r-o", linewidth=2,
            label="IRLS (Newton)")
ax.set_xlabel("Iteration", fontsize=11)
ax.set_ylabel("$\\|\\beta^{(t)} - \\beta^*\\|$", fontsize=11)
ax.set_title("Convergence comparison: GD vs IRLS", fontsize=12)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("logistic_convergence.png", dpi=150, bbox_inches="tight")
plt.show()

収束速度の比較から、2つの最適化手法の根本的な違いが明確に読み取れます。

  1. IRLSの二次収束: 対数スケールで見ると、IRLSの誤差は各反復でおよそ半分の傾きで直線的に減少しています。これは二次収束（有効桁数が毎回倍増）の特徴です。わずか5-6回の反復で機械精度に達しています

  2. 勾配降下法の線形収束: 勾配降下法は対数スケールで一定の傾きで減少しており、線形収束（各反復で誤差が一定の比率で減少）を示しています。100回の反復でもIRLSの数回分の精度に到達していません

  3. 計算量のトレードオフ: IRLSは各反復で $p \times p$ の連立方程式を解くコスト（$O(p^3)$）がかかりますが、反復回数が少ないため全体の計算量は小さくなります。$p$ が大きい場合は勾配降下法（各反復 $O(np)$）のほうが有利になることがあり、問題の規模に応じた選択が重要です

まとめ

本記事では、ロジスティック回帰の理論を解説しました。

  • モデル: $P(y=1|\bm{x}) = \sigma(\bm{x}^\top\bm{\beta})$ で、対数オッズが説明変数の線形関数です。シグモイド関数が線形予測子を確率に変換します
  • 対数オッズ: $\ln(p/(1-p)) = \bm{x}^\top\bm{\beta}$ という変換により、値域の制約と線形モデルの整合性を両立しています
  • 最尤推定: 対数尤度 $\ell = \sum[y_i z_i – \ln(1+e^{z_i})]$ を最大化します。勾配は $\bm{X}^\top(\bm{y} – \bm{p})$ という簡潔な形で、閉形式解がないため反復法が必要です
  • 凹性: ヘッセ行列 $-\bm{X}^\top\bm{W}\bm{X}$ が半負定値であるため、対数尤度は凹関数であり、局所最適は大域最適です
  • IRLS: ニュートン法は反復重み付き最小二乗法として実装でき、二次収束により高速に収束します
  • オッズ比: $e^{\beta_j}$ は変数 $x_j$ が1単位増えたときのオッズの倍率を表し、リスク要因の効果の定量化に使われます
  • 多クラス拡張: ソフトマックス関数により $K > 2$ クラスへの自然な拡張が可能です

次のステップとして、以下の記事も参考にしてください。

関連記事

前の記事

多変量正規分布の性質と条件付き分布の完全導出

次の記事

重回帰分析の理論と最小二乗法の導出

関連タグ: シグモイド関数 ロジスティック回帰 分類 多変量解析 最尤推定 統計学