PRML 4.3.2. ロジスティック回帰

はじめに

パターン認識と機械学習（PRML）4章においては、ロジスティック回帰を用いることによって2クラス分類が可能であることが説明されていた。輪読の際に4.3.2-4.3.3の内容が少し難しかったため、実際にRでの具体例を交えながら理解を深めることとした。

4.3.2-4.3.3. ロジスティック回帰

本章においてはロジスティック回帰モデルのパラメータ$\mathbf{w}$を最尤法を用いて決定しようとしている。最尤法を用いてパラメータを推定するには、尤度関数（4.89）を算出し、算出された尤度関数をもとに交差エントロピー誤差関数（4.90）を導出する。この（4.90）式が最小となるパラメータ$\mathbf{w}$を求めることができれば最尤法の完了となる。$\mathbf{w}$を求める最も簡単な方法は（4.90）式を$\mathbf{w}$で偏微分し、傾きが0となる$\mathbf{w}$を求める方法である。これは（3.1.1節）での最尤推定と全く同じ流れとなる。

\begin {align*}
p(\mathbf{t} \mid \mathbf{w})=\prod_{n=1}^N y_n^{t_n}\left\{1-y_n\right\}^{1-t_n}
\tag{4.89}
\end {align*}

\begin {align*}
E(\mathbf{w})=-\ln p(\mathbf{t} \mid \mathbf{w})=-\sum_{n=1}^N\left\{t_n \ln y_n+\left(1-t_n\right) \ln \left(1-y_n\right)\right\}
\tag{4.90}
\end {align*}

\begin {align*}
\nabla E(\mathbf{w})=\sum_{n=1}^N\left(y_n-t_n\right) \boldsymbol{\phi}_n
\tag{4.91}
\end {align*}

しかし、（4.91）式における$\ y_{n}$は$\ y_{n} = \sigma(\mathbf{w}^{T}\mathbf{\phi}_{n})$であり、$\mathbf{w}$に依存している。そのため、（4.91）式の左辺を0とおいても容易に$\mathbf{w}$が求まらないことが推測され、実際解析的に求めることはできない。ただ、（演習4.15）によって、（4.90）の交差エントロピー誤差関数は凸関数であり、唯一の最小解を持つことが証明される。解析的に求まらずとも、ニュートンラフソン法を用いて$\mathbf{w}$を繰り返し更新することによって唯一の最小解を求めることができる。
上記の一連の流れを実際にR上でデータのシミュレーションを行い、実行してみる。

まず、2次元データのシミュレーションにおいてはRパッケージの「MASS」を活用して作成する。

library(MASS)
# set the parameters
n <- 50 # the number of data
sigma <- diag(2) * 0.3 # covariance matrix

# simulate the data
set.seed(20230109)
blue1 <- mvrnorm(n = n, 
                 mu = c(-1, -1), 
                 Sigma = sigma)

blue2 <- mvrnorm(n = n, 
                 mu = c(1, 1), 
                 Sigma = sigma)

red1 <- mvrnorm(n = 2 * n, 
                mu = c(0, 0), 
                Sigma = sigma)

data <- rbind(blue1, blue2, red1)
class <- rep(c("Blue", "red"), each = 2 * n)
df <- data.frame(data, Class = class)
plot(df$X1, df$X2, col = df$Class, pch = 19, cex = 0.7, 
     xlab = "x1", ylab = "x2")

上図のデータに対して、ガウス基底関数を適用し、データを新たな特徴空間に写す。

ガウス基底関数

\begin {align*}
\begin {array}{l}
\phi_{1n} = \exp\Big\lbrace{-\frac{(x_{1n} - \mu_{11}) ^ 2}{2s ^ 2}\Big\rbrace} + \exp\Big\lbrace{-\frac{(x_{2n} - \mu_{12}) ^ 2}{2s ^ 2}\Big\rbrace} \\
\phi_{2n} = \exp\Big\lbrace{-\frac{(x_{1n} - \mu_{21}) ^ 2}{2s ^ 2}\Big\rbrace} + \exp\Big\lbrace{-\frac{(x_{2n} - \mu_{22}) ^ 2}{2s ^ 2}\Big\rbrace}
\end {array}
\end {align*}

$\mathbf{\phi}^T_1 = (\phi_{11}, ..., \phi_{1N})$, $\ \mathbf{\phi}^T_2 = (\phi_{21}, ..., \phi_{2N})$,
$\mathbf{x}^T_1 = (x_{11}, ..., x_{1N})$, $\ \mathbf{x}^T_2 = (x_{21}, ..., x_{2N})$,
$\mu_{11} = \mu_{12} = -1$, $\ \mu_{21} = \mu_{22} = 0$, $s = 0.6$　（$N$は総データ数）
上述の特徴空間に新たに投影してみる。

# Gaussian function
Gauss <- function(x, mu, s) {
  y <- exp(-(x - mu)^2 / (2 * s^2))
  return(y)
}

s <- 0.6
gauss1 <- Gauss(df$X1, mu = -1, s = s) / 2 + Gauss(df$X2, mu = -1, s = s) / 2
gauss2 <- Gauss(df$X1, mu = 0, s = s) / 2 + Gauss(df$X2, mu = 0, s = s) / 2
plot(gauss1, gauss2, col = df$Class, pch = 19, cex = 0.7, 
     main = "Gaussian function", xlab = "phi1", ylab = "phi2")

本データを用いて実際に逐次アルゴリズムとニュートンラフソン法をそれぞれ適用してクラス分類を行ってみる。

逐次アルゴリズム

4.3.2章で少し触れられていたが、データが与えられる度にモデルの更新を行う逐次学習も可能である。その際には、(4.91)式を用いて

\begin {align*}
\mathbf{w}^{(\tau+1)}=\mathbf{w}^{(\tau)}-\eta \nabla E_n
\tag{3.22}
\end {align*}

上式（3.22）を作成し、$\mathbf{w}$を逐次更新することができる。ただし、本手法は初期値と$\eta$を適切な値に調整する必要がある。とりあえず$\mathbf{w} = 0$, $\eta = 0.01$で実行してみる。

w <- c(0, 0, 0)
eta <- 0.01

# Logistic sigmoid function
LogSigma <- function(a) {
  logSigma <- 1 / (1 + exp(-a))
  return(logSigma)
}

phai0 <- rep(1, 4 * n)
phai1 <- gauss1
phai2 <- gauss2
phai <- cbind(phai0, phai1, phai2)
# 0 means blue, 1 means red
t <- rep(c(0, 1), each = 2 * n)
# The stochastic gradient descent algorithm
for (n in sample(4 * n, replace = F)) {
  # n <- 1
  a_n <- t(w) %*% phai[n, ]
  y_n <- LogSigma(a_n)
  deltaE_n <- c((y_n - t[n])) * phai[n, ]
  w <- w - eta * deltaE_n
}

# calculating the predicted value
a <- t(w) %*% t(phai)
y <- LogSigma(a)

# visualize the border (p = 0.5)
p1 <- seq(0, 1, 0.01)
p2 <- -(p1 * w[2] + w[1]) / w[3]
plot(gauss1, gauss2, col = df$Class, pch = 19, cex = 0.7, 
     main = "Gaussian function", xlab = "phi1", ylab = "phi2")
points(p1, p2, type = "l")

黒色の直線が決定境界（$p = 0.5$）となるが、あまりうまくクラス分類ができていないことが分かる。当然初期値や$ \ \eta \ $の値を妥当なものにすればある程度はうまくいくが、手間がかかる。

ニュートンラフソン法

こちらも繰り返しパラメータの更新を行うが、先ほどの逐次更新と異なり最初からすべてのデータを活用する。ニュートンラフソン法の更新式は以下の式で表される。

\begin {align*}
\mathbf{w}^{(\text {new })}=\mathbf{w}^{\text {(old) }}-\mathbf{H}^{-1} \nabla E(\mathbf{w})
\tag{4.92}
\end {align*}

（4.92）式における$\nabla E(\mathbf{w})\ $および$ \ \mathbf{H} \ $は以下の式で表される。

\begin {align*}
\nabla E(\mathbf{w})=\sum_{n=1}^N\left(y_n-t_n\right) \boldsymbol{\phi}_n=\boldsymbol{\Phi}^{\mathrm{T}}(\mathbf{y}-\mathbf{t})
\tag{4.96}
\end {align*}

\begin {align*}
\mathbf{H}=\nabla \nabla E(\mathbf{w})=\sum_{n=1}^N y_n\left(1-y_n\right) \boldsymbol{\phi}_n \boldsymbol{\phi}_n^{\mathrm{T}}=\boldsymbol{\Phi}^{\mathrm{T}} \mathbf{R} \Phi
\tag{4.97}
\end {align*}

\begin {align*}
R_{n n}=y_n\left(1-y_n\right)
\tag{4.98}
\end {align*}

（4.92）（4.96）（4.97）を用いることで更新式を求めることができる（計算省略）。実際に実装してみる。

#### 4.3.3. Iterative reweighted least squares ####
w <- c(0, 0, 0)
diff <- 10
while (diff > 0.0001) {
  print(w)
  a <- t(w) %*% t(phai)
  y <- LogSigma(a)
  
  # weighted diagonal matrix
  R <- diag(c(y * (1 - y)))
  z <- phai %*% w - solve(R) %*% t(y - t)
  
  # calculating new "w"
  new_w <- solve(t(phai) %*% R %*% phai) %*% t(phai) %*% R %*% z
  diff <- sum(abs(new_w - w))
  w <- new_w
}

# visualize the border (p = 0.5)
p1 <- seq(0, 1, 0.01)
p2 <- -(p1 * w[2] + w[1]) / w[3]
plot(gauss1, gauss2, col = df$Class, pch = 19, cex = 0.7, 
     main = "Gaussian function", xlab = "phi1", ylab = "phi2")
points(p1, p2, type = "l")

非常にうまく決定境界を求められていることが分かる。また、繰り返し回数も5回程度であり収束は非常に速い。データが段階的に得られる場合には、最初にニュートンラフソン法を適用してある程度妥当な$\mathbf{w}$の値を求めておき、データが得られる度に逐次更新を行う手法も考えられる。