ロジスティック回帰

ロジスティック回帰についてのメモです．内容はPythonで学ぶあたらしい統計学の教科書 ，機械学習のエッセンスを参考にしています．

ロジスティック回帰は，確率分布に二項分布を用いて，リンク関数にロジッド関数を用いた一般線形化モデルで，主に二値分類に用いられるアルゴリズムです．説明変数は複数あっても構わず，連続型・カテゴリがたが混在していても支障はありません．

ここで，ラベルの値は0 or 1とします．与えられた特徴量のサンプル${\bf x} \in {\mathcal{R}}^{d}$に対して，ラベル$y$が1になる確率を$P(Y = 0 | X = x)＄で表します．
ロジスティック回帰は次式を前提としたモデルです．

$P(Y = 1 | X = x)  = \sigma (w_{0} + \sum_{j = 1}^{d} x_{j} w_{j}) = \sigma ({\bf{w}}^{T} {\bf{\tilde{x}}}^{T})$

$\sigma$の中は${\bf{x}}$の線形関数となっており，このアイデアは，サンプル$x$がラベル1に属する確からしさは線形関数の値が大きいほど大きくなるという仮定です．ここで，線形関数であることから，その値は$- \infty$から$\infty$の間をとることになります．その線形関数の値にシグモイド関数を作用させることで確率として扱うことができるようになります．
ここで，$\sigma$はシグモイド関数で

$\sigma(\eta) = \cfrac{1}{1 + e^{-\eta}}$

で定義されます．

特徴量行列$X$とラベルのベクトル${\bf{y}}$が与えられたときに，その$X$から${\bf{y}}$が生じる確率を考えます．
行列$X$のi番目のサンプルがラベル$y_{i}$に分類される確率を全て掛け算したものとして考えるこができて，以下のようになります．

$P(y | x) =  \prod_{k=1}^{n} [ \sigma ({\bf{\tilde{x_{k}}}} {\bf{w}})^{y_{k}} (1 - \sigma({\bf{\tilde{x_{k}}}}) {\bf{w}}))^{1-y_{k}} ]$

この確率を最大化することを考えます（このままでは考えにくいので，対数をとってマイナスをつけると下式のようになります）．

$E({\bf{w}}) = - \log P ( {\bf{y}} | X) = - \sum_{k = 1}^{n} [ y_{k} \log \sigma ( {\bf{w}}^{T} {\tilde{\bf{x_{k}}}}) + (1 - y_{k}) \log ( 1 - \sigma ({\bf{w}}^{T} {\tilde{x_{k}}})) ]$

この$E({\bf{w}})$の最適値をニュートン法で求めます．つまり$\nabla E({\bf{w}}) = 0$の解をニュートン法で求めることになります．

これをPythonで実装すると以下のようになります．
ファイル名：logisticreg.py

import numpy as np

from scipy import linalg

THRESHMIN = 1e-10

def sigmoid(x):

    return 1 / (1 + np.exp(-x))

class LogisticRegression:

    def __init__(self, tol=0.001, max_iter=3, random_seed=0):

        self.tol = tol

        self.max_iter = max_iter

        self.random_state = np.random.RandomState(random_seed)

        self.w_ = None

    def fit(self, X, y):

        self.w_ = self.random_state.randn(X.shape[1] + 1)

        Xtil = np.c_[np.ones(X.shape[0]), X]

        diff = np.inf

        w_prev = self.w_

        iter = 0

        while diff > self.tol and iter < self.max_iter:

            yhat = sigmoid(np.dot(Xtil, self.w_))

            r = np.clip(yhat * (1 - yhat),

                        THRESHMIN, np.inf)

            XR = Xtil.T * r

            XRX = np.dot(Xtil.T * r, Xtil)

            w_prev = self.w_

            b = np.dot(XR, np.dot(Xtil, self.w_) -

                       1 / r * (yhat - y))

            self.w_ = linalg.solve(XRX, b)

            diff = abs(w_prev - self.w_).mean()

            iter += 1

    def predict(self, X):

        Xtil = np.c_[np.ones(X.shape[0]), X]

        yhat = sigmoid(np.dot(Xtil, self.w_))

        return np.where(yhat > .5, 1, 0)

ファイル名：logisticreg_wdbc.py

import logisticreg

import csv

import numpy as np

n_test = 100

X = []

y = []

with open("wdbc.data") as fp:

    for row in csv.reader(fp):

        if row[1] == "B":

            y.append(0)

        else:

            y.append(1)

        X.append(row[2:])

y = np.array(y, dtype=np.float64)

X = np.array(X, dtype=np.float64)

y_train = y[:-n_test]

X_train = X[:-n_test]

y_test = y[-n_test:]

X_test = X[-n_test:]

model = logisticreg.LogisticRegression(tol=0.01)

model.fit(X_train, y_train)

y_predict = model.predict(X_test)

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

print("Accuracy: {}/{} = {}".format(n_hits, n_test, n_hits / n_test))

実行すると以下のようになります（実行する際には，logisticreg_wdbc.pyとwdbc.dataが同じディレクトリにある必要があります）．

$ python3 logisticreg_wdbc.py

Accuracy: 97/100 = 0.97

Juliaで実装すると以下のようになります．
ファイル名：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

ファイル名：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)")

実行すると以下のようになります（実行する際には，logisticreg_wdbc.pyとwdbc.dataが同じディレクトリにある必要があります）．

$ Julia logisticreg_wdbc.jl
Accuracy: 97/100 = 0.97