「データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC」（通称「緑本」）を読み直しました。一般化線形モデルと最尤法によるパラメータ推定の説明がとても分かりやすかったです。
この記事では二項ロジスティック回帰を例にとって最尤法をPythonで実装してみて、statsmodelsと同じパラメータの推定値が求められることを示します。
セットアップ
import numpy as np
import scipy.optimize as opt
import statsmodels.api as sm
import statsmodels.formula.api as smf
import pandas as pd
import plotnine as p9
from mizani.breaks import breaks_width
記事を書いた環境は以下のとおりです。
- python: 3.13.5
- mizani: 0.14.1
- numpy: 2.3.2
- scipy: 1.15.3
- statsmodels: 0.14.4
- pandas: 2.3.1
- plotnine: 0.14.6
二項ロジスティック回帰
今回使うサンプルデータの説明です。
いま、生徒20人のクラスでテストを行ったとします。次のDataFrameは生徒ごとのテストの合否と勉強状況を記録したものです。idは生徒の出席番号、xはテスト前日の勉強時間（単位は時間）、yはテストの合否（合格なら1、不合格なら0）とします。
df = pd.DataFrame({
"id": [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20],
"x": [0.2, 0.6, 1.6, 1.6, 1.8, 1.8, 2.1, 2.9, 3.0, 3.7, 4.3, 5.2, 6.0, 6.0, 7.1, 7.3, 8.3, 8.7, 9.5, 9.7],
"y": [0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1],
})
(
p9.ggplot(df, p9.aes(x="x", y="y"))
+ p9.theme_classic()
+ p9.geom_point()
+ p9.scale_x_continuous(breaks=breaks_width(2))
+ p9.scale_y_continuous(breaks=breaks_width(1))
)
テスト前日の勉強時間が長いほど合格する傾向にあるように見えるので、テスト前日の勉強時間から合否を予測するモデルが作れそうです。
いま、生徒$i (i = 1, \dots, 20)$について、以下のモデルを仮定します。
\begin{align}
z_{i} &= \beta_{0} + \beta_{1} x_{i} \\
p_{i} &= \frac {1} {1 + \exp(-z_{i})} \\
y_{i} &\sim Bernoulli(p_{i})
\end{align}
これが二項ロジスティック回帰です。今回は簡単のため、また図で分かりやすく示す都合上説明変数はxの一つとしましたが、二つ以上あっても構いません。
$y$は、合格か不合格かという「コイン投げ」のような結果として得られる値なので、ベルヌーイ分布に従います。そのため、$y$を$x$で線形回帰するのではなく、$y$の分布であるベルヌーイ分布の母数である確率$p$を一旦挟むのがポイントです。
$p$を$x$の線形和で表すと、$p$は確率の条件である$0 \leq p \leq 1$が保証されません。ロジスティック関数$f(x) = 1 / (1 + \exp(-x))$は$0 \leq f(x) \leq 1$ですから、$z$を$x$の線形和で表し、$z$をロジスティック関数に入れてあげることで$0 \leq p \leq 1$を保証するわけです。
二項ロジスティック回帰は、Pythonならstatsmodelsを用いて簡単に実装できます。
# add_constantしないと、beta_0がないモデルを推定する
X = sm.add_constant(df["x"])
y = df["y"]
mod = sm.Logit(y, X)
fit = mod.fit()
# print(fit.summary())でもよいが、パラメータの表だけ表示したいので
print(fit.summary().tables[1])
Optimization terminated successfully.
Current function value: 0.394625
Iterations 7
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const -2.8610 1.299 -2.202 0.028 -5.407 -0.315
x 0.7522 0.315 2.389 0.017 0.135 1.369
==============================================================================
Rのformula記法で書きたい場合はこちらでも構いません。
mod = smf.logit("y ~ x + 1", data=df)
fit = mod.fit()
print(fit.summary().tables[1])
printされたconstとxはそれぞれ式の$\beta_{0}, \beta_{1}$にあたります。
それでは実際のところ、$\beta_{0}, \beta_{1}$はどうやって推定されるのでしょうか？
最尤法
パラメータを与えたときの手持ちのデータへのフィット度合いを尤度といいます。今回の例では、パラメータは$\beta_{0}, \beta_{1}$です。
パラメータ（1個でも複数でも構いません）を$\theta$とすると、尤度は以下の式で表されます。
$$
L(\theta) = \prod_{i} P(Y = y_{i} | \theta)
$$
ベルヌーイ分布は、$y=1$となる確率が$p$, $y=0$となる確率が$1-p$の分布なので、ベルヌーイ分布の確率質量関数は以下のように表されます。
$$
P(Y = y) = p^{y} (1 - p)^{(1 - y)} (y = 0, 1)
$$
すなわち、今回のモデルの尤度は次で計算されます。
$$
L(\beta_{0}, \beta_{1}) = \prod_{i} p_{i}^{y_{i}} (1 - p_{i})^{(1 - y_{i})}
$$
尤度はパラメータの関数であり、今回の例では$p$は$\beta_{0}, \beta_{1}$の関数なので、尤度関数も$\beta_{0}, \beta_{1}$の関数です。
尤度関数が最大になる$\beta_{0}, \beta_{1}$が最尤法で求められる最尤推定量です。一見難しいことをしているように見えますが、例えば偏ったコインを100回投げて表が30回出たならこのコインの表が出る確率は30%と見積もるのと同じ考え方です。
尤度関数は積の形であり計算しづらいため、実装のときは対数を取った以下の対数尤度で計算します。
\begin{align}
\log L &= \log (\prod_{i} p_{i}^{y_{i}} (1 - p_{i})^{(1 - y_{i})}) \\
&= \sum_{i} (y_{i} \log p_{i} + (1 - y_{i}) \log (1 - p_{i}))
\end{align}
なお、線形回帰では最小二乗法が有名ですが、線形回帰で誤差項が独立で同一な正規分布に従う場合、最尤法と最小二乗法は同じ結果が導かれます（参考:
回帰分析におけるパラメータ推定：最尤法と最小二乗法 #統計学 -
Qiita）。
最尤法の実装
求めたいのは対数尤度$\log L$が最大となるときの$\beta_{0}, \beta_{1}$ですが、解析解を導出できるとは限らないので数値計算で求めます。scipy.optimize.minimizeは目的関数を最小化するものなので、$-\log L$を最小化する問題として扱います。
最適化アルゴリズムには、勾配の計算が不要なNelder-Mead法を用います。あとでstatsmodelsと比較できるよう、収束判定に用いるパラメータも指定しておきます。
def logistic(x):
return 1 / (1 + np.exp(-x))
def negative_log_likelihood(w, x, y):
"""
対数尤度関数の-1倍を計算する
Args:
w: 長さ2（z = \beta_0 + \beta_1 x の \beta_0と\beta_1）
x: 長さn
y: 長さn
Returns:
float
"""
p = logistic(w[0] + w[1] * x)
return -1 * np.sum(y * np.log(p) + (1-y) * np.log(1-p))
res_opt = opt.minimize(
negative_log_likelihood,
x0=[0, 0],
args=(df["x"], df["y"]),
method="Nelder-Mead",
options={
"maxiter": 1000,
"xatol": 1e-6,
"fatol": 1e-6,
}
)
print(f"beta_0 = {res_opt.x[0]:.4f}, beta_1 = {res_opt.x[1]:.4f}, logL = {-res_opt.fun:.4f}")
beta_0 = -2.8610, beta_1 = 0.7522, logL = -7.8925
ロジスティック回帰の説明変数の個数がいくつであってもnegative_log_likelihoodが動くように行列形式で書くこともできますが、行列ではなくスカラーで書く方が個人的には説明が平易に思えるので、wは説明変数が1個と切片があるモデルにのみ対応するようにしました。
実行してみると$\beta_{0} = -2.8610, \beta_{1} = 0.7522$と推定されました。そのときの対数尤度は$-7.8925$でした。
次にstatsmodelsで実装してみましょう。上で示したコードのとおり、fitの引数は省略してデフォルト値のままでも大抵うまくいきますが、今回は自作実装での対数尤度関数の最適化と条件を揃えるために引数を指定します。
statsmodelsのドキュメントを読むと、statsmodels.Logit.fitはデフォルトではNewton-Raphson法を用いますが、method="minimize"を指定するとscipyのminimizeをバックエンドで使用でき、そのあとに**kwargsでminimizeの引数を与えられるとあります。自作実装と揃えるため、method="minimize"としてoptionsも同じ値を入れてあげます。また、最適化探索のための初期値start_paramsも同じ値を入れます。
X = sm.add_constant(df["x"])
y = df["y"]
mod = sm.Logit(y, X)
fit = mod.fit(
start_params=[0, 0],
method="minimize",
min_method="Nelder-Mead",
maxiter=1000,
xatol=1e-6,
fatol=1e-6
)
print(fit.summary().tables[1])
print(f"logL = {fit.llf:.4f}")
Optimization terminated successfully.
Current function value: 0.394625
Iterations: 89
Function evaluations: 169
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const -2.8610 1.299 -2.202 0.028 -5.407 -0.315
x 0.7522 0.315 2.389 0.017 0.135 1.369
==============================================================================
logL = -7.8925
constとx（それぞれ$\beta_{0}, \beta_{1}$）ともに自作実装と同じ値が得られました。対数尤度も一致していますね。
参考文献
Join Us
株式会社JPX総研では、データとデジタルの力で証券取引所の新しいビジネスを作る仲間を募集しております。Wantedlyをご覧ください。