『機械学習のエッセンス(http://isbn.sbcr.jp/93965/)』のPythonサンプルをJuliaで書き換えてみる。（第05章06ロジスティック回帰）

# はじめに

[『機械学習のエッセンス(http://isbn.sbcr.jp/93965/)』のPythonサンプルをJuliaで書き換えてみる。（第05章05ラッソ回帰）](https://qiita.com/yshutaro/items/5818e4fd2ab9e976c3b1)の続きです。

# ロジスティック回帰

※だいぶ省略しているので、詳細は書籍をご参照ください。

与えられた特徴量のサンプル$x \in \mathbb{R}^d$に対して、ラベル$y$が1になる確率を$P(Y = 1|X = x)$で表し、$y$が0になる確率を$(Y = 0|X = x)$で表す。

```math
P(Y = 1|X = x) = \sigma(w_0 + \sum_{j=1}^d x_jw_j) = \sigma(w^T\tilde{x}^T)
```

$\sigma$はシグモイド関数と呼ばれ、下記で定義される。

```math
\sigma(\xi) = \frac{1}{1 + e^{-\xi}}
```

ラベルの値が$y$になる確率は次のようになる。

```math
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$が生じる確率を考える。

```math
P(y|X) = \prod_{k=1}^n[\sigma(\tilde{x}_kw)^{y_k}(1-\sigma(\tilde{x}_kw))^{1-y_k}]
```

これの確率を最大化することを考える。$log$を取ってマイナスを付けたものを$E(w)$として考える。

```math
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))]
```

この最適化をニュートン法によって求める。つまり

```math
\nabla E(w) = 0
```

という方程式の解をニュートン法で求める。$E(w)$のヘッセ行列$H = \nabla^2E(w)$を計算する。
$E(w)$の1階微分は

```math
\nabla E(w) = \sum_{k=1}^n (\sigma(\tilde{x}_k w) - y_k)\tilde{x}_k^T
```

となる。ここで

```math
\hat{y} = (\sigma(w^T \tilde{x}_1),\sigma(w^T \tilde{x}_2),\dots,\sigma(w^T \tilde{x}_n))^T
```

とおくと、行列$\tilde X$を使って次のように書ける。

```math
\nabla E(w) = \tilde{X}^T(\hat{y} - y)
```

ヘッセ行列$H$の$(i, j)$成分$H_{ij}$は次のようになる。

```math
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$を次の式で定義する。

```math
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$は次で表せる。

```math
H = \tilde{X}^T R \tilde{X}
```

$\nabla E(w) = 0$となる$w$をニュートン法で求める。


```math
w^{new} = w^{old} - H^{-1} \nabla E(w^{old})
```

これを変形して

```math
w^{new} = (X^T R X)^{-1}(X^T R)[X w^{old} - R^{-1}(\hat{y} - y)]
```

# ロジスティック回帰の実装

```logisticreg.jl
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サンプルと比較していきました。

```julia
XR = Xtil' .* r'
```

のところは、Xtilのサイズが`(469, 31)`なのでXtil'だと転置なので`(31, 469)`、rのサイズが`(469,)`となります。
PythonのサンプルでXRの型を表示すると`(31, 469)`となっています。単純に`*`でかけ算するとXRの型が`(31,)`になってしまい一致しません。そのため、ｒも転置して`.*`によりブロードキャストでかけ算することで`(31, 469)`になります。このあたりは数式とサンプルを読んだだけではわからず難しいところでした。

# 実際のデータに適用

本と同じく、下記のデータを取得し、プログラムと同じディレクトリに置きます。

https://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/wdbc.data

```logisticreg_wdbc.jl
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
julia> include("logisticreg_wdbc.jl")
Accuracy: 97/100 = 0.97
```

本と同じく100件中97件が当たっている結果となりました。


## 補足

### Numpyの`clip`の部分

Numpyの`clip`は、Juliaでは`clamp`でした。

https://docs.julialang.org/en/v1.1/base/math/#Base.Math.clamp

```julia
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
julia> clamp(x , 2, 7)
ERROR: MethodError: no method matching isless(::Int64, ::Array{Int64,1})
（略）
```

### 配列の比較

本のサンプル`logisticreg_wdbc.py`では最後の配列の一致を下記のように完結に書いています。

```python
n_hits = (y_test == y_predict).sum()
```

これはNumPｙの配列比較を利用しています。

```python
>>> 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
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
julia> a .== b
6-element BitArray{1}:
  true
 false
  true
 false
 false
  true

julia> sum(a .== b)
3
```