コード実装 忘備録: "もう少し深めようかな？" ロジスティック回帰 編

初めに

機械学習のリサーチペーパーを読み、アルゴリズムをコンピューター言語で実装する力がほしいと思いこのシリーズを始めました。高度の数学は必要ありませんが、理解できていれば、データ構造、動作理解、更にはトラブルシューティングに役に立つはずです。

今回はロジスティック回帰モデルの実装に関する要点を見ていきます。

全ての実装コードを書きだすより、分からなかった場所を書きだし、忘備録とします。また、線形回帰モデルとの違いを比較すると、理解が一層進みます。

ここで取り上げたいのは、以下になります。

• 線形判別関数 (Linear Discriminant)
• ヘヴィサイド関数 (Heaviside Function)
• シグモイド関数(Sigmoid Function)
• 最尤推定(Maximum Likelihood Estimation:MLE)

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

ロジスティック回帰モデルのアルゴリズムは回帰 (Regression) のように読めますが、実際には、分類 (Binary classification) に使われます。

具体的には、目的変数 (Response variable) が 0 または 1 や 'Spam' または 'NotSpam' の二値となる場合に用いられる手法です。

数学的基礎

1. 線形回帰モデル (Linear Regression Model)

線形回帰モデルでは、説明変数 (Features)である $x_i$ の各測定値と係数 $β_i$ の線形結合が、そのまま目的変数 (Response variable) $y$（連続値）に対応しています。具体的には、以下のように表現されます。

y = β_0 + β_1 x_1 + β_2 x_2 + \cdots + β_d x_d + \varepsilon

• $x_1, x_2$, $\ldots$, $x_d$) : 説明変数の測定値
• $β_0$ : 切片 (intercept)
• $β_1$, $β_2$, $\ldots$, $β_d$ : 各説明変数に対応する係数 (重み)
• $\varepsilon$ : 誤差項（Residual）

2. ロジスティック回帰モデル (Logistic Regression Model)

一方、ロジスティック回帰モデルでは、各観測に対して説明変数（Features）$x_1, x_2, \dots, x_d$ を基に、$y=1$ である確率$p$を求めます。ここでいう $y=1$ は、二値の一つを表します。つまり、

p = P(y = 1 \mid x)

と定義できます。$y=1|x$は, $x$ の条件下で、可能な値の一つを取るという意味で、0 か 1 や、Spam, NotSpam を表しています。$P(y = 1 \mid x)$は、x の条件下でその事象が発生する確率 $(0≤𝑝≤1)$ を表しています。

この確率 $p$ は 0 から 1 の値を取る必要があるため、直接

z = \beta_0 + \beta_1 x_1 + \cdots + \beta_d x_d

という線形結合を用いるのではなく、シグモイド関数を用いて変換します。

3. シグモイド関数(Sigmoid Function)

シグモイド関数は

\sigma(z) = \frac{1}{1+\exp(-z)}

と定義され、どんな実数 $z$ に対しても $\sigma(z)$ は必ず 0 から 1 の範囲に収まります。ここで$z$は上記の線形結合を表しています。

実装では、線形結合の$z$を計算し、この値をシグモイド関数に適用し、$y=1$ である確率 $p$ は

p = \sigma(z) = \frac{1}{1+\exp\left(-\left(\beta_0 + \beta_1 x_1 + \cdots + \beta_d x_d\right)\right)}

の数式で求めます。これがロジスティック回帰モデルの主要な形です。

余談ですが、線形回帰の場合、$y$ はそのまま連続値の予測値として目的変数 (response variable) に利用されますが、ロジスティック回帰の場合は、$z$ は中間的な値となり、シグモイド関数のインプットになり、最終的にシグモイド関数の $p$ が目的変数 (response variable) になります。

実装

ここから実装に入ります。通常のデータを読み込んでライブラリ－を利用するのではなく、どのような関数が実装されているのかを見ていきます。

Part1 - 拡張ベクトル(Augumented Vector)

線形回帰でもよく使われる手法で、各観測値に「切片 (Intercept) 」を含めるために、最初の座標として常に1を追加します。各観測値は、もともとの$d$個の特徴量に加えて、ダミー変数（固定値1）を追加して、拡張ベクトルとします。以下が数式となります。

x_i^T \equiv (x_{i,0}, x_{i,1}, \ldots, x_{i,d-1}, 1)

• $x_{i,0}, x_{i,1}, \ldots, x_{i,d-1}$ は元々の $d$ 個の特徴（説明変数）の値
• 末尾の $1$ は切片項を導入するための固定値

次に行列Xを以下のように定義します：

\begin{align*}
X \equiv 
\begin{pmatrix}
x^{0T} \\
x^{1T} \\
\vdots \\
x^{m-1T}
\end{pmatrix}
\end{align*}

単純に、行ベクトルが列ベクトルになっただけと理解できます。

なぜこんなことをするのでしょうか？

1. 数式を単純化出来ます
2. データ行列と列ベクトルの計算に落とし込めますので、一括計算でき、コードが高速化でき、単純化できます

例をあげましょう。

観測値1: $x^{0T} = (2,5,1)$
観測値2: $x^{1T} = (3,6,1)$
観測値3: $x^{2T} = (7,8,1)$

の場合、各観測値の拡張ベクトルを行としてまとめると、データ行列 $X$ は以下のようになります。

X =
\begin{pmatrix}
2 & 5 & 1 \\
3 & 6 & 1 \\
7 & 8 & 1
\end{pmatrix}

この行列のサイズは $(3 \times 3)$ で 3行で、観測値が元々2つと1の切片項(Intercept)で3列になります。

例えば、モデルの係数 (Coefficents) $\beta$ が

\beta =
\begin{pmatrix}
\beta_0\\
\beta_1\\
\beta_2
\end{pmatrix}

とした場合、予測値 $\hat{y}$（線形予測子）は

\hat{y} = X\beta

と一括で計算できます。これは計算が速く、コードもシンプルになります。ロジスティック回帰の場合の $z$ でも同じです。

$\hat{y} = X\beta$ を演算するコードをPythonで実装してみましょう。

import numpy as np

# Design matrix X: each row is an augmented vector [x₁, x₂, 1] (example)
X = np.array([[2, 5, 1],
              [3, 6, 1],
              [7, 8, 1]])

# Coefficient vector beta (reshaped as a column vector)
beta = np.array([0.4, 0.3, 0.2]).reshape(-1, 1)

# Calculate the linear predictor y_pred (predicted values)
y_pred = X @ beta

print(y_pred)

np.arrayを利用し、@でdot.productを計算しています。

Part 2 - 線形判別関数 (Linear Discriminant)

この関数は、マトリックス$X$に係数ベクトル$\theta$の内積を計算し、その符号(正か負）で２つのクラス(１ か ０）のクラスに分けるシンプルな方法です。

まず、例として２つの観測値（２次元）を考えてみましょう。

• 拡張ベクトル $x$:

各観測値は、本来の特徴（例：$x_0, x_1$）に加えて、切片項（固定で1）を付加します。これを拡張ベクトルと呼びます。

$$
x = \begin{pmatrix}
x_0 \
x_1 \
1
\end{pmatrix}
$$

• 係数ベクトル $\theta$:

係数(Coefficients)のベクトルです。

$$
\theta = \begin{pmatrix}
\theta_0 \
\theta_1 \
\theta_2
\end{pmatrix}
$$

• 線形判別関数の計算:
  内積 $\theta^T x$ により、スカラー値 $z$ を得ます。つまり、

$$
z = \theta^T x = \theta_0 x_0 + \theta_1 x_1 + \theta_2
$$

決定境界と分類:

• 決定境界: $z = 0$ のとき、すなわち

$$
\theta_0 x_0 + \theta_1 x_1 + \theta_2 = 0
$$

となる点の集合は、入力空間上で直線を形成し、クラスの境界線となります。

• 分類ルール:
  • $z > 0$ → クラス1（例）
  • $z \leq 0 $ → クラス0（例）

内積 $\theta^T x$ の値が正か負か（つまり、入力 $x$ が重みベクトル $\theta$ とどの方向にあるか）によって、データ点がどちら側に位置するかを判断することができます。

具体的な数値で考えてみると、たとえば、$\theta$が

\theta = \begin{pmatrix}
1\\
1\\
-3
\end{pmatrix}

とすると、決定境界の式は

$$
1 \cdot x_0 + 1 \cdot x_1 - 3 = 0 \quad \Longrightarrow \quad x_0 + x_1 = 3
$$

となります。これを $x_1$ について解くと

x_1 = 3 - x_0

となります。これは、2次元平面上における直線の方程式です。緑の線となるはずです。

この直線 $x_1 = 3 - x_0$ は、 内積 $ \theta^T x $ が正（例えば、点 $(2,2)$ の場合、計算すると $2 + 2 - 3 = 1$ で正 なら、一方のクラスに分類し、内積が負（例えば、点 $(1,1)$ の場合、計算すると $1 + 1 - 3 = -1$ で負なら、もう一方のクラスに分類します。

Part 3 - ヘヴィサイド関数 (Heaviside Function)

ヘヴィサイド関数は、各観測値について、$x$ と $\theta$ との内積 $z = \theta^T x$ を計算し、その $z$ の値が入力となります。 Heaviside 関数は、この $z$ の各値について、

• $z > 0$ の場合は $1.0$ を出力し、
• $z \le 0$ の場合は $0.0$ を出力

します。クラスラベルの値はFloat型です。

つまり、各観測値が $\theta^T x$ によって得られるスコアがゼロ以上かどうかをチェックして、対応するクラスのラベル（0,1）に変換しています。

線形判別関数とヘヴィサイド関数をPythonで実装してみましょう。

import numpy as np

def linear_discriminant(X: np.ndarray, theta: np.ndarray) -> np.ndarray:
    # Returns
    # -------
    # numpy.ndarray, shape (m, 1)
    #    The result of the dot product (linear discriminant) for each observation.
    return X @ theta

# Define the coefficient vector θ (as a column vector)
theta = np.array([1, 1, -3]).reshape(-1, 1)

# Construct the design matrix X where each row is an augmented observation [x0, x1, 1]
X = np.array([
    [2, 2, 1],   # Observation 1: Expected z = 1
    [1, 1, 1],   # Observation 2: Expected z = -1
    [3, 0, 1],   # Observation 3: Expected z = 0
    [4, 2, 1]    # Observation 4: Expected z = 3
])

# Linear discriminant
# Calculate the linear discriminant for each observation
z = linear_discriminant(X,theta)

print("Design Matrix (X):")
print(X)
print("\nCoefficient Vector (θ):")
print(theta)
print("\nComputed Linear Discriminant (z = θ^T x):")
print(z)

# Define the Heaviside function (elementwise)
def heaviside(y: np.ndarray) -> np.ndarray:
    # Return 1.0 if element > 0, else 0.0.
    return np.where(y > 0, 1.0, 0.0)

# Apply the Heaviside function to the computed z-values to get final labels
labels = heaviside(z)
print("\nPredicted Labels (using Heaviside function):")
print(labels)

アウトプットは以下のようになります。

Design Matrix (X):
[[2 2 1]
 [1 1 1]
 [3 0 1]
 [4 2 1]]

Coefficient Vector (θ):
[[ 1]
 [ 1]
 [-3]]

Computed Linear Discriminant (z = θ^T x):
[[ 1]
 [-1]
 [ 0]
 [ 3]]

Predicted Labels (using Heaviside function):
[[1.]
 [0.]
 [0.]
 [1.]]

Part 4 - ロジスティック関数　(Logistic Function)

ロジスティック関数はシグモイド関数(Sigmoid Function)とも呼ばれる関数を利用します。1 又は 0 のクラスラベルを 0 から 1 までの確率に直す必要があります。

シグモイド関数の定義は

\sigma(z) = \frac{1}{1+\exp(-z)}

でした。

ロジスティック関数は、まず各観測に対して説明変数の線形結合 $z$ を計算し、この $z$ にシグモイド関数を適用します。すると $y=1$ である確率 $p$ は

p = \sigma(z) = \frac{1}{1+\exp\left(-\left(\beta_0 + \beta_1 x_1 + \cdots + \beta_d x_d\right)\right)}

と表せます。ここでいう $y=1$ は、クラスラベルの値を表します。

ロジスティック関数をPythonで実装してみましょう。

要点は、NumPy の vectorize 関数を利用することで、カスタム関数（ここでは sigmoid 関数）を配列全体に対して要素ごとに適用する関数を定義しています。

def logistic(Y: np.ndarray) -> np.ndarray:
    # Define a helper function that computes the sigmoid for a single value.
    def sigmoid(zi: float) -> float:
        from math import exp
        sig = 1 / (1 + exp(-zi))
        return sig

    # The defintion of the python function that uses np.vectorize to create a vectorized
    # version of the sigmoid function
    logistic_values = np.vectorize(sigmoid)

    # Apply the python function that calls sigmoid function and return the results
    return logistic_values(Y)

mpl.rc("savefig", dpi=120) # Adjust for higher-resolution figures
sns.set_style("darkgrid")
y_pos = y_values[y_values > 0]
y_rem = y_values[y_values <= 0]
plt.plot(y_rem, heaviside (y_rem), 'b')
plt.plot(y_pos, heaviside (y_pos), 'b')
plt.plot(y_values, logistic (y_values), 'r--')

実行できるコードは以下に置いておきます。Github code

Part 5 - 最尤推定関数(Maximum Likelihood Estimation Function)

ロジスティック回帰の最大尤度推定（MLE）は、観測データが実際に観測される確率（尤度）を最大にするパラメータ（係数）を求める推定方法です。具体的な流れと考え方は以下の通りです。

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

上記で説明していますが、ロジスティック回帰モデルでは、各観測値（例：各観測 拡張ベクトル $ \hat{x}^i $）に対して、線形結合を計算します。

z_i = \theta^T \hat{x}^i = \beta_0 + \beta_1 x_{i,1} + \cdots + \beta_d x_{i,d}

その値 $ z_i $ にシグモイド関数（ロジスティック関数） $ \sigma(z) $ を適用して、

g_i = \sigma(z_i) = \frac{1}{1+\exp(-z_i)}

$ y_i = 1 $ となる確率 $ g_i $ を得ます。

• $ \hat{x}^i $ : 観測データのi番の拡張ベクトル。元の観測データが[1, 2, 3]であれば、$ \hat{x}^1$ は =[1,2,3,1]

• $ y_i $ は各観測に対するラベル

• $ g_i $ は各観測 $ \hat{x}^i $ のシグモイド関数 G を適用した結果として計算される確率

観測データの確率モデル（ベルヌーイ分布)

各観測 $ y_i $ は、ベルヌーイ試行として扱います。つまり、$ y_i $ は 1 または 0 の値をとり、

• $ y_i = 1 $ であれば、その確率は $ g_i $
• $ y_i = 0 $ であれば、その確率は $ 1 - g_i $

となります。

これを一般化すると、各観測 $ i $ に対して条件付き確率は

\Pr(y_i \mid \hat{x}^i, \theta) = g_i^{\, y_i} \, (1-g_i)^{\, 1-y_i}

と表せます。

全データに対する尤度関数と対数尤度

全 $ m $ 個の観測が独立（i.i.d.）であると仮定すると、全体の尤度 $ L(\theta) $ は各観測の条件付き確率の積となります：

L(\theta) = \prod_{i=0}^{m-1} \Pr(y_i \mid \hat{x}^i, \theta)
          = \prod_{i=0}^{m-1} \Bigl[ g_i^{\, y_i} \, (1-g_i)^{\, 1-y_i} \Bigr]

計算を簡単にするために、対数をとって対数尤度（log-likelihood） $ \ell(\theta) $ を求めます：

\ell(\theta) = \ln L(\theta) = \sum_{i=0}^{m-1} \Bigl[ y_i \ln g_i + (1-y_i) \ln(1-g_i) \Bigr]

この対数尤度は、各観測での貢献（$ y_i \ln g_i $ や $ (1-y_i) \ln(1-g_i) $ の計算）を総和として表せます。

MLEの目的：対数尤度の最大化

MLE（最大尤度推定）の目的は、この対数尤度関数 $ \ell(\theta) $ を最大化する $ \theta $ を求めることです。 すなわち、実際の観測データ $ y $ が最も発生しやすくなるようなパラメータ $ \theta $ を探します：

\hat{\theta} = \arg\max_{\theta} \, \ell(\theta)

これにより、モデルは「このパラメータなら観測されたデータを最もよく説明できる」という最尤推定の考え方に基づいて学習されます。

行列記法を用いた表現

実際の計算では、多数の観測に対して同じ処理（各観測の $ y_i \ln g_i $ や $ (1-y_i) \ln(1-g_i) $ の計算）を行います。

この総和の操作を行列記法でコンパクトに書くために、

• 観測ラベルベクトル $ y $
• すべての観測に対する $ g = G(X\theta) $（ここで $ X $ はデザインマトリックス）
• そして、すべての成分が1 のベクトル $ u $

を用いて、いくつかの記法がありますが、たとえば次の形の一例があります（複数の導出方法が存在します)：

\ell(\theta) = y^T \ln\Bigl( G(X\theta) \Bigr) + (u-y)^T \ln\Bigl[ u - G(X\theta) \Bigr]

ここで $\ln$ や $G$ は、行列またはベクトルに対して要素ごとに適用されるものと解釈します。

また、もう一つの表現として、

\ell(\theta) = y^T X\theta + u^T \ln\left[ G(-X\theta) \right]

があります。

これは、

\ell(\theta) = y^T X\theta + u^T \ln\left[ G(-X\theta) \right]

以下の要素が含まれます。

• $y$: 観測されたクラスラベルの列ベクトルです。各成分 $y_i$ は通常、0 または 1 の値。

• $X$: デザインマトリックスです。各行は、拡張された入力ベクトル $\hat{x}^i$ を表しており、元の特徴量に加えて切片項として1が付加されています。

• $\theta$: パラメータ（係数）ベクトルです。このベクトルの値を最尤推定（MLE）により求め、観測データを最もよく説明する決定境界を作り出します。

• $X\theta$: 行列 $X$ とベクトル $\theta$ の積です。結果は各観測に対する内積 $\theta^T \hat{x}^i$ となり、各観測の線形予測子（スコア）を表します。

• $G(-X\theta)$: ここで $G(z)$ はシグモイド関数（ロジスティック関数）を意味し、

  G(z) = \frac{1}{1+\exp(-z)}
  

  です。$X\theta$ で求めた各観測の線形予測子に対して、マイナスをかけたものにシグモイド関数 $G$ を要素ごとに適用することで、各観測で「$y_i=1$ となる確率」の補数、つまり「$y_i=0$ となる確率」に対応する値が得られます。

• $\ln\Bigl[G(-X\theta)\Bigr]$
  シグモイド関数 $G(-X\theta)$ の各要素に対して自然対数 $\ln$ を要素ごとに適用したものです。これにより、各観測における対数尤度の「$y_i=0$ に関する部分」が表現されます。

• $u$: すべての成分が1の列ベクトルです。長さは観測数 $m$ に合わせたもので、

  u = \begin{pmatrix} 1 \\ 1 \\ \vdots \\ 1 \end{pmatrix}
  

  となります。これは、対数尤度の各観測における項の総和を取るために用いられます。たとえば、

  u^T \ln\Bigl[G(-X\theta)\Bigr]
  

  は、すべての観測の $\ln\Bigl[G(-X\theta)\Bigr]$ の値の和を計算する操作と同じです。

このように記号を使って総和の計算を一括して表現することで、対数尤度の式がよりコンパクトになり、数学的な解析や実装が楽になります.

MLE関数をPythonで実装してみましょう。

要点は、行列で理解できる上記のような数式にたどり着き、要素を理解し、numpy.NDArrayで実装することです。

import numpy as np
from numpy.linalg import multi_dot
from numpy.typing import NDArray

def log_likelihood(theta: NDArray, y: NDArray, X: NDArray) -> float:

    """
    Computes the log-likelihood of the logistic regression model.

    Parameters:
    - theta: NDArray of shape (d+1, 1)
        Coefficient (parameter) vector including intercept term.
    - y: NDArray of shape (m, 1)
        Vector of observed binary labels (0 or 1).
    - X: NDArray of shape (m, d+1)
        Design matrix where each row is an augmented feature vector (with bias term).

    Returns:
    - float
        The scalar value of the log-likelihood for the current theta.
    """

   # Vector of ones (used to compute the full sum across all observations)
    u = np.ones(y.shape)

    # Compute the linear scores z = -X @ theta
    minus_X_theta = np.dot(-X, theta)

    # Apply the sigmoid function element-wise: G = sigmoid(-X @ theta)
    # Ensure that we have double minus operations
    G = 1 / (1 + np.exp(-minus_X_theta))  # shape (m, 1)

    # Compute log-likelihood: yᵀXθ + uᵀ log(G(-Xθ))
    l_theta_b = multi_dot([y.T, X, theta]) + np.dot(u.T, np.log(G))
    # l_theta_b is <class 'numpy.ndarray'>

    return l_theta_b.item()  # return a scalar value (float) from numpy.ndarray

# Example usage:
if __name__ == "__main__":
    # Create an example design matrix X (each row is an augmented observation)
    X_example = np.array([[1, 2, 1],
                          [3, 4, 1]])
    
    # Example parameter vector theta (as a column vector)
    theta_example = np.array([0.5, 0.5, 0.1]).reshape(-1, 1)
    
    # Example observed labels y (as a column vector)
    y_example = np.array([[1], [0]])
    
    # Compute the log-likelihood for the example data
    ll_value = log_likelihood(theta_example, y_example, X_example)
    print("Log-likelihood:", ll_value)

アウトプットは以下になります。

Log-likelihood: -3.8108578338965464

終わりに

いかがでしたでしょうか？自分も手探りで書いておりますので、間違っている箇所があるかもしれません。もし間違いや指摘等があれば、コメントしていただけると幸いです。

数式の実装が面白いと感じられたらうれしいです。行列やベクトルで考え、Numpyでどうやって高速実装するかが肝のように思えます。

Happy Hacking!!!