はじめに
『機械学習のエッセンス(http://isbn.sbcr.jp/93965/)』のPythonサンプルをJuliaで書き換えてみる。（第05章05ラッソ回帰）の続きです。
ロジスティック回帰
※だいぶ省略しているので、詳細は書籍をご参照ください。
与えられた特徴量のサンプル$x \in \mathbb{R}^d$に対して、ラベル$y$が1になる確率を$P(Y = 1|X = x)$で表し、$y$が0になる確率を$(Y = 0|X = x)$で表す。
P(Y = 1|X = x) = \sigma(w_0 + \sum_{j=1}^d x_jw_j) = \sigma(w^T\tilde{x}^T)
$\sigma$はシグモイド関数と呼ばれ、下記で定義される。
\sigma(\xi) = \frac{1}{1 + e^{-\xi}}
ラベルの値が$y$になる確率は次のようになる。
P(Y = y | X = x) = P(Y = 1 | X = x)^y P(Y = 0 | X = x)^{1-y} \\
= \sigma(\tilde{x}^Tw)^y(1-\sigma(\tilde{x}^Tw))^{1-y}
特徴量行列$X$とラベルベクトル$y$が与えられたとき、その$X$から$y$が生じる確率を考える。
P(y|X) = \prod_{k=1}^n[\sigma(\tilde{x}_kw)^{y_k}(1-\sigma(\tilde{x}_kw))^{1-y_k}]
これの確率を最大化することを考える。$log$を取ってマイナスを付けたものを$E(w)$として考える。
E(w) = -logP(y|X) = -\sum_{k=1}^n[y_k log\sigma(w^T\tilde{x}_k) + (1 - y_k)log(1-\sigma(w^T\tilde{x}_k))]
この最適化をニュートン法によって求める。つまり
\nabla E(w) = 0
という方程式の解をニュートン法で求める。$E(w)$のヘッセ行列$H = \nabla^2E(w)$を計算する。
$E(w)$の1階微分は
\nabla E(w) = \sum_{k=1}^n (\sigma(\tilde{x}_k w) - y_k)\tilde{x}_k^T
となる。ここで
\hat{y} = (\sigma(w^T \tilde{x}_1),\sigma(w^T \tilde{x}_2),\dots,\sigma(w^T \tilde{x}_n))^T
とおくと、行列$\tilde X$を使って次のように書ける。
\nabla E(w) = \tilde{X}^T(\hat{y} - y)
ヘッセ行列$H$の$(i, j)$成分$H_{ij}$は次のようになる。
H_{ij} = \sum_{k=1}^n \hat{y}_k(1 - \hat{y}_k)x_{ki}x_{kj}
ここで任意の$i = 1,2,\dots,n$に対して、$x_{i0} = 1$と定義する。
単純化のため対角行列$R$を次の式で定義する。
R = \left(
\begin{array}{ccccc}
\hat{y}_1(1 - \hat{y}_1) \\
& & \hat{y}_2(1 - \hat{y}_2) & & \\
& & & \ddots & \\
& & & & \hat{y}_n(1 - \hat{y}_n)
\end{array}
\right)
$H$は次で表せる。
H = \tilde{X}^T R \tilde{X}
$\nabla E(w) = 0$となる$w$をニュートン法で求める。
w^{new} = w^{old} - H^{-1} \nabla E(w^{old})
これを変形して
w^{new} = (X^T R X)^{-1}(X^T R)[X w^{old} - R^{-1}(\hat{y} - y)]
ロジスティック回帰の実装
module logisticreg
using LinearAlgebra
using Random
using Statistics
THRESHMIN = 1e-10
sigmoid(x) = 1 ./ (1 .+ ℯ.^(-x))
mutable struct LogisticRegression
tol
max_iter
random_seed
w_
function LogisticRegression(tol, max_iter=3, random_seed=0)
new(tol, max_iter, random_seed, Nothing)
end
end
function fit(s::LogisticRegression, X, y)
Random.seed!(s.random_seed)
s.w_ = randn(size(X)[2] + 1)
Xtil = hcat(ones(size(X)[1]), X)
diff = Inf
w_prev = s.w_
iter = 0
while diff > s.tol && iter < s.max_iter
yhat = sigmoid(Xtil * s.w_)
r = clamp.(yhat .* (1 .- yhat), THRESHMIN, Inf)
XR = Xtil' .* r'
XRX = (Xtil' .* r') * Xtil
w_prev = s.w_
b = XR * (Xtil * s.w_ .- (1 ./ r) .* (yhat - y))
s.w_ = XRX \ b
diff = mean(abs.(w_prev - s.w_))
iter = iter + 1
end
end
function predict(s::LogisticRegression, X)
Xtil = hcat(ones(size(X)[1]), X)
yhat = sigmoid(Xtil * s.w_)
return [ifelse(x .> .5, 1, 0) for x in yhat]
end
end
補足
下記の実際のデータに適用で実行を行ったところ、途中の行列演算で型が合わずエラーが出たので、途中でそれぞれの型をPythonサンプルと比較していきました。
XR = Xtil' .* r'
のところは、Xtilのサイズが(469, 31)なのでXtil'だと転置なので(31, 469)、rのサイズが(469,)となります。
PythonのサンプルでXRの型を表示すると(31, 469)となっています。単純に*でかけ算するとXRの型が(31,)になってしまい一致しません。そのため、ｒも転置して.*によりブロードキャストでかけ算することで(31, 469)になります。このあたりは数式とサンプルを読んだだけではわからず難しいところでした。
実際のデータに適用
本と同じく、下記のデータを取得し、プログラムと同じディレクトリに置きます。
include("logisticreg.jl")
using .logisticreg
using CSV
n_test = 100
dataframe = CSV.read("wdbc.data", header=false)
row,col=size(dataframe)
X = Float64[dataframe[r,c] for r in 1:row, c in 3:col]
y = Float64[ifelse(dataframe[r,2] == "B", 0, 1) for r in 1:row]
y_train = y[1:end-n_test]
X_train = X[1:row-n_test, :]
y_test = y[end-n_test+1:end]
X_test = X[end-n_test+1:end, :]
model = logisticreg.LogisticRegression(0.01)
logisticreg.fit(model, X_train, y_train)
y_predict = logisticreg.predict(model, X_test)
n_hits = sum(y_test .== y_predict)
println("Accuracy: $n_hits/$n_test = $(n_hits / n_test)")
実行結果
julia> include("logisticreg_wdbc.jl")
Accuracy: 97/100 = 0.97
本と同じく100件中97件が当たっている結果となりました。
補足
Numpyのclipの部分
Numpyのclipは、Juliaではclampでした。
julia> x = collect(0:9)
10-element Array{Int64,1}:
0
1
2
3
4
5
6
7
8
9
julia> clamp.(x , 2, 7)
10-element Array{Int64,1}:
2
2
2
3
4
5
6
7
7
7
.（ドット）なしのclampだけだとエラーとなります。
julia> clamp(x , 2, 7)
ERROR: MethodError: no method matching isless(::Int64, ::Array{Int64,1})
（略）
配列の比較
本のサンプルlogisticreg_wdbc.pyでは最後の配列の一致を下記のように完結に書いています。
n_hits = (y_test == y_predict).sum()
これはNumPｙの配列比較を利用しています。
>>> import numpy as np
>>> a = np.array([1,2,3,4,5,6])
>>> b = np.array([1,5,3,3,2,6])
>>> a == b
array([ True, False, True, False, False, True])
>>> (a == b).sum()
3
同じことをJuliaで書くと下記のようになってうまくいきません。（NumPyの配列ではなくPythonの配列も下記と同じでした。） 上記のNumPyの様に書けないか調べたのですが、わからず、for文で書きました。 @yatra9 さんにご教示いただきました!下に追記してあります。logisticreg_wdbc.jl も修正済みです。
julia> a = [1,2,3,4,5,6]
6-element Array{Int64,1}:
1
2
3
4
5
6
julia> b = [1,5,3,3,2,6]
6-element Array{Int64,1}:
1
5
3
3
2
6
julia> a == b
false
julia> sum(a == b)
0
（追記）
a .== b として.（ドット）付きのブロードキャストで比較できました。
julia> a .== b
6-element BitArray{1}:
true
false
true
false
false
true
julia> sum(a .== b)
3