Python でロジスティック回帰をするには、少なくとも3通りの方法があります：
ここでは scikit-learn を使うことにします。
Rのロジスティック回帰もご覧ください。
100点満点の模擬テストを20人が受験しました。点数は次の通りです：
x = [4, 5, 17, 17, 19, 19, 29, 41, 45, 47, 49, 54, 56, 58, 70, 76, 88, 88, 93, 97]
この20人がその後実際の試験を受験して、合否が判明しました。合 = 1、否 = 0 で表すと、次のようになりました：
y = [0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1]
プロットしてみましょう：
import matplotlib.pyplot as plt plt.figure(figsize=(8, 3)) # 図の横・縦。なくてもよい plt.plot(x, y, "ok")
x から y を予測する式を求めたいのですが、通常の線形回帰ではうまくあてはまりませんし、そもそも何を予測したいのかはっきりしません。
予測するものを、模擬テストの点数が x のときの合格確率だとしましょう。確率は負になったり 1 を超えたりしませんので、線形回帰ではうまくありません。うまく 0 から 1 の範囲におさまってくれる関数はないでしょうか。
例えば
\[ y = \frac{1}{1 + \exp(-x)} = \frac{1}{1 + e^{-x}} \]はどうでしょうか。これをグラフにすると、次のようになります：
import numpy as np a = np.linspace(-5, 5) b = 1 / (1 + np.exp(-a)) plt.figure(figsize=(8, 3)) plt.plot(a, b)
これなら縦軸の値は必ず 0 から 1 の範囲に収まりますので、確率と解釈することができます。実際には、このグラフの中心（変曲点の位置）と曲がり具合とを変えられるようにした
\[ y = \frac{1}{1 + \exp(-(w_1 x + w_0))} \]という式で表される曲線を使うことにします。
この形の曲線をロジスティック曲線といいます。より一般に、0 から 1 に滑らかに変化する S 字形をつぶしたような曲線をシグモイド曲線といいます。シグモイド曲線は、ニューラルネットの活性化関数として使われてきました（現在は ReLU の類を使うのが普通です）。
$x$ が与えられたとき $y = 1$ になる確率 $p$ が
\[ p = \frac{1}{1 + \exp(-(w_1 x + w_0))} \]で近似できるような係数 $w_1$、$w_0$ を求める方法をロジスティック回帰といいます。ただ、線形の回帰分析のときのような最小2乗法ではなく、最尤法を使います。つまり、$y = 1$ になるデータ点での $p$ と $y = 0$ になるデータ点の $1-p$ を全部掛け算したものが最大になるように係数を決めます。
上の式を書き直すと
\[ \log\left(\frac{p}{1-p}\right) = w_1 x + w_0 \]となりますが、この左辺に現れる $p$ の関数をリンク関数といいます。この場合のリンク関数はロジット関数と呼ばれ、$\mathrm{logit}(p)$ と書くことがあります。ちなみに $p/(1-p)$ はオッズ（odds）と呼ばれる量ですので、ロジット関数はオッズの対数ということになります。
scikit-learn ではロジスティック回帰は次のようにして行います：
from sklearn.linear_model import LogisticRegression X = np.array(x).reshape(-1, 1) model = LogisticRegression(penalty=None).fit(X, y) print(model.coef_, model.intercept_)
ここで、np.array() は Python のリストを NumPy の配列に変換する関数です。これに .reshape(m, n) を付けると m 行 n 列の行列になります。m と n はどちらか一つを指定すれば十分です。指定しない側を -1 で表します。つまり .reshape(-1, 1) は1列だけの行列に変換します。ここでは行列を表す変数は大文字にしています。np.array(x).reshape(-1, 1) は np.column_stack([x]) と書くこともできます。
結果は次のようになります：
[[0.08318673]] [-3.56667708]
これが $w_1$、$w_0$ の値です。つまり、
\[ p = \frac{1}{1 + \exp(-(0.08319 x - 3.56668))} \]ということになります。元のデータとこの確率を重ねてグラフにしてみます：
plt.figure(figsize=(8, 3)) plt.plot(x, y, "ok") a = np.arange(0, 101) plt.plot(a, 1 / (1 + np.exp(-(0.08319 * a - 3.56668))))
scikit-learn の LogisticRegression() は、デフォルトでは、係数が 0 に近くなるような正則化（regularization）が行われます。ベイズ推定のことばでいえば、係数が 0 に近くなるような事前分布が設定されています。正則化のパラメータ C はデフォルトでは 1 に設定されており、この値が小さいほど強い正則化が行われます。昔の scikit-learn では、正則化をやめるには C に非常に大きな値（例えば C=np.inf）を設定するしかありませんでした。scikit-learn 0.21 で penalty='none' というオプションで正則化を無効にできるようになりました。scikit-learn 1.2 で penalty='none' は廃止予定（deprecated）になり、今後は penalty=None を使うことが推奨されています。
次の表を考えましょう。ここで $y$ は入試の合否（合格なら1、不合格なら0）、$x_1$ は内申点、$x_2$ は模試偏差値を表すとします。
| $y$ | $x_1$ | $x_2$ |
|---|---|---|
| 0 | 3.6 | 60.1 |
| 1 | 4.1 | 52.0 |
| 0 | 3.7 | 62.5 |
| 0 | 4.9 | 60.6 |
| 1 | 4.4 | 54.1 |
| 0 | 3.6 | 63.6 |
| 1 | 4.9 | 68.0 |
| 0 | 4.8 | 38.2 |
| 1 | 4.1 | 59.5 |
| 0 | 4.3 | 47.3 |
さきほどと同様に、$x_1$、$x_2$ が与えられたとき $y = 1$ になる確率 $p$ が
\[ p = \frac{1}{1 + \exp(-(w_1 x_1 + w_2 x_2 + w_0))} \]で近似できるようにパラメータ $w_1$、$w_2$、$w_0$ を定めます。
まずはデータを Python で表します：
x1 = [3.6, 4.1, 3.7, 4.9, 4.4, 3.6, 4.9, 4.8, 4.1, 4.3] x2 = [60.1, 52.0, 62.5, 60.6, 54.1, 63.6, 68.0, 38.2, 59.5, 47.3] y = [0, 1, 0, 0, 1, 0, 1, 0, 1, 0]
このデータについてロジスティック回帰を行うには
X = np.column_stack([x1, x2]) model = LogisticRegression(penalty=None).fit(X, y) print(model.coef_, model.intercept_)
と打ち込みます。結果は
[[1.27158497 0.06424052]] [-9.44589515]
つまり $\mathrm{logit}(p) = 1.27158x_1 + 0.06424x_2 - 9.44590$ で予測されることがわかります。
Last modified: