「自然科学の統計学」を読む　～第8章質的データの統計分析 ①ロジット・モデルとプロビット・モデル～

はじめに

東京大学教養学部統計学教室編「自然科学の統計学」は，1992年発行のやや古典的な文献であるが，自然科学に関わる統計学的テーマが簡潔にまとめられている。数理統計学の復習も兼ねて，本書を読むこととした。
ただ，基本的なことは他書で学んできたのと，本書自体がかなり細かく説明されているので，本書内の内容や数式を細かく追うというより，実務や統計検定の受験において有用そうなことを選んでまとめてみたい。

本記事は，「第8章質的データの統計分析」における，ロジット・モデルとプロビット・モデルに関する読書メモである。

第8章質的データの統計分析

本章では質的変数，すなわち「買う・買わない」，「合格する・合格しない」，「(医学・生物学での)生きる・死ぬ」といった変数を扱う。また，これらを扱うための手法であるプロビットモデル・ロジットモデルなどについて説明している。

8.1 二値データ

先ほど示した「買う・買わない」，「合格する・合格しない」，「(医学・生物学での)生きる・死ぬ」のように，対象がある状態やカテゴリーの項目に属しているかどうかのみを知ることができるデータを質的データと呼ぶ。

質的データを扱う方法として代表的なものは，分割表を用いる方法や，回帰分析の考え方を応用した確率モデルがある。この確率モデルとして代表的なものが，

プロビット・モデル
ロジット・モデル(ロジスティック・モデル)

である。

8.2 ロジット・モデルとプロビット・モデル

プロビット・モデル

プロビット・モデルは，標準正規分布の累積分布関数

$\begin{align} \Phi(z) = \int_{- \infty}^z \frac{1}{\sqrt{2} \pi} e^{-x^2/2} dx \\ \\ \end{align}$

を使い， $Y_i=1$ となる確率を

$\begin{align} F^*(X_i) = \Phi(\beta_0 + \beta_1 X_i) \\ \\ \end{align}$

とするものである。

ロジット・モデル

ロジット・モデルは，ロジスティック分布の累積分布関数

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

を使い， $Y_i=1$ となる確率を

$\begin{align} F^*(X_i) = \Lambda(\beta_0 + \beta_1 X_i) \\ \\ \end{align}$

とするものである。

ロジットモデルにおいて，

$\begin{align} P_{0i}=P(Y_i=0 \mid X_i), \quad P_{1i}=P(Y_i=1 \mid X_i) \\ \\ \end{align}$

とすると，

$\begin{align} P_{0i} &= 1 - P_{1i} = \frac{1}{1 + \exp (\beta_0 + \beta_1 X_i)} \\ \\ \Rightarrow \frac{ P_{1i} }{ P_{0i} } &= \frac{ P_{1i} }{ 1- P_{1i} } = \exp (\beta_0 + \beta_1 X_i) \\ \\ \therefore \log \left( \frac{ P_{1i} }{ P_{0i} } \right) &= \beta_0 + \beta_1 X_i \\ \\ \end{align}$

となり，これを対数オッズと呼ぶ。ロジットモデルは，対数オッズが線形関数になるモデルである。

標準正規分布およびロジスティック分布のそれぞれの累積分布関数を並べて描くと，左図のようになる。
また標準正規分布および分散1のロジスティック分布のそれぞれの確率密度関数を並べて描くと，右図のようになる。

描画用コードはクリックで展開

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, logistic

# Define the functions
def logistic_cdf(z):
    return 1 / (1 + np.exp(-z))

def logistic_pdf(z):
    d = np.sqrt(3)/np.pi
    return (1/d) * np.exp(-z/d) / (1 + np.exp(-z/d))**2

# Generate x values
z = np.linspace(-5, 5, 1000)

# Calculate y values for each function
logistic_cdf_y = logistic_cdf(z)
norm_cdf_y = norm.cdf(z)
logistic_pdf_y = logistic_pdf(z)
norm_pdf_y = norm.pdf(z)

# Create the plot
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Plot CDFs
axes[0].plot(z, logistic_cdf_y, label='Logistic CDF')
axes[0].plot(z, norm_cdf_y, label='Normal CDF', linestyle='dashed')
axes[0].set_title('Cumulative Distribution Functions')
axes[0].set_xlabel('z')
axes[0].set_ylabel('Probability')
axes[0].legend()
axes[0].grid(True)

# Plot PDFs
axes[1].plot(z, logistic_pdf_y, label='Logistic PDF')
axes[1].plot(z, norm_pdf_y, label='Normal PDF', linestyle='dashed')
axes[1].set_title('Probability Density Functions')
axes[1].set_xlabel('z')
axes[1].set_ylabel('Density')
axes[1].legend()
axes[1].grid(True)

plt.tight_layout()
plt.show()

プロビット・モデルとロジット・モデルの比較

確率密度関数

ロジスティック分布の分散は $\sigma^2 = \pi^2/3$ になることが知られており，分散が1になるように正規化すると，ロジスティック分布の確率密度関数は $d = \sqrt{3}/\pi$ として，

$\begin{align} f_L(z) = \frac{1}{d} \cdot \frac{\exp(z/d)}{(1 + \exp(z/d))^2} \\ \\ \end{align}$

となる。

パラメータの推定値

図に示す通り，変数 $z$ の標準正規分布の確率密度関数と，分散1のロジスティック分布の確率密度関数はかなり近い形になる。分散1のロジスティック分布は，変数 $z/d$ のロジスティック分布とみなせるので，ロジット・モデルによるパラメータの推定値は，プロビット・モデルのパラメータの推定値のほぼ $1/d = \pi / \sqrt{3}$ となる。