今回は機械学習の最も基礎的な手法にであるロジスティック回帰を使ってあやめのデータである iris の分類器を作ってみます。
まず、機械学習やロジスティック回帰の内容について説明してから、最後にR言語を用いて実装を行なっています。実装は機械学習ライブラリをあまり使わず行列計算を使って実装しているので数式と結びやすいのではないかと思います。
Contents
0. なぜロジスティック回帰を行うか
さまざまなライブラリ（PyTorch, Tensor Flow, Scikit-learn）などを用いれば簡単に畳み込みニューラルネットワーク（CNN）や時系列データ学習に用いるRecurrent Neural Network （RNN）、ニューラルネット以外ではランダムフォレストやサポートベクターマシーン（SVM）などを使って機械学習を行うことができます。しかし、基本的な仕組みがわからなかったり、理解しようとすると初学では難しかったりします。そこで今回は深層ニューラルネットワークにつなげる初歩としてロジスティック回帰について扱ってみようと思いました。ロジスティック回帰は単層のパーセプトロンと等価でありこれを多層化したものが多層パーセプトロン（MLP）です。
1. 機械学習についておさらい
軽く機械学習についておさらいしてみます。機械学習とはデータからパターンなどを自動で学習し性能を改善するプログラムです。
機械学習は主なものに分けると教師あり学習、教師なし学習、強化学習に分けることができます。最も理解しやすくはじめに学習するのは教師あり学習であることが多いです。教師あり学習の中では分類問題、回帰問題と分かれています。今回は分類問題でも最も簡単な2クラス分類を扱います。分類問題では離散化されたものに分類するものでしたが回帰問題は連続値を予測します。分類問題は区切る仕切りを学習するのに対して回帰問題は乗っている関数を見つけるというような感じで、実験データのプロットから近似曲線を書くのと似ています。
それぞれの用途に応じて機械学習の手法はさまざまあり上で例として出した、ニューラルネットやランダムフォレスト、サポートベクターマシーン、k-means、ロジスティク回帰、線形回帰などがあります。近年ではニューラルネットを深層化することで精度が向上しておりdeep learnig がよく話題に出ますが深層化も機械学習の手法の一つです。
2. ロジスティック回帰を用いた学習
2.1 ロジスティック回帰の構造
今回用いる iris のデータは4つの特徴量があります。4つの情報を0か1をきめる一つの0~1の1つの値に出力します。0~1の値に対してどこまでが0でどこからが1なのかを設定しその値に対して予測結果を出力します。その際それぞれの入力を結果に反映する大きさとして重み wとバイアス b を加えてありこれが学習するパラメータとなっています。
数式的にはこのように計算されます。
1つのデータに対し4つのデータがあり重みwをかけバイアスbを足し合わせたものを活性化関数に通します。今回用いる活性化関数はシグモイド関数といい2値分類をするときに用います。
データ数に対して行列計算で表すとこのようになります。
2.2 活性化関数(sigmoid関数)
シグモイド関数は以下のように定義される関数で入力の値に対して0~1の値に変換することができます。
グラフで表すと
このようなグラフであり、単調増加し入力に対する値が一対一対応で0~1の範囲の値に変えることができます。シグモイド関数の特徴は分類する際の境界を設定することができるところです。通常は0.5を境界にして2つの値を分類します。code に起こした場合は次のようになります。
ここでオーバーフローを回避するために入力がせいの場合は定義通り計算し負の場合は
と変形し右側の式を用いて計算します。
2.3 目的関数（交差エントロピー誤差関数）
機械学習のモデルと実際のデータとの関係がわからなければ学習することができません。そこで機械学習モデルと実際のデータとの違いを表すのが目的関数です。目的関数は損失関数と呼ばれ最小にしていく場合や、生成モデルを扱う場合では尤度関数と呼ばれ最大化させることもありますが本質的には同じことを意味していて機械学習モデルがデータとどれくらい違うのか近いのかを表しこれに従ってパラメータを更新していきます。目的関数も非常に多くのものがありモデルによって様々ですが分類問題では交差エントロピー誤差関数を用います。交差エントロピー誤差関数は二つの確率分布の差が小さいほど小さな値を返し確率分布が離れているほど大きな値を返す関数となっています。2値分類の場合以下のように書けます。
2.4 モデルの学習（勾配降下法）
精度を上げていくために、目的関数を最小の値にする方向に勾配降下法を用いてパラメータw、b の値を更新していきます。
目的関数の勾配を求め値が小さくなる方向に更新していくことを勾配降下法といいます。交差エントロピー誤差関数は凸関数であり勾配が 0 となるところで最小値を取ることがわかっています。つまり理論的に勾配が 0 となるパラメータ w と b が求められれば最も最適なモデルといえます。交差エントロピー誤差関数は勾配が 0 で最小値ですが、他のモデルや高次元なモデルの場合勾配が 0 になるところが最小値にならず最適な解が得られない場合があります。よって勾配降下法にもこのような問題を解決するために様々な手法が存在します。
今回の場合は最も単純に勾配降下法であり、勾配降下に関係するパラメータは
であり学習率を表していて一度に変化させる大きさを指定しています。
具体的に計算していくと交差エントロピー誤差関数の勾配
は非常に簡単な式で表すことができ以下のように計算することができます。
以下で計算してみます。まずそれぞれの勾配はこのように書けます。
それぞれについて求めると
これを行列計算しそれぞれのデータで求めデータの平均をとることで更新するパラメータの勾配とします。
3. 実装
では、実際にR言語でロジスティク回帰を使いあやめのデータセット iris の分類を行なっていきます。
3.1 前処理
まず取得データの確認を行います。
4つの特徴量がありデータラベルは3種類あり、それぞれ25個のデータとなっています。今回は2クラス分類のため上の50個のデータを用いてsetosaとversicolorの分類を行います。
codeは以下のようになります。また下のエディタで前処理後のデータを確認してみてください。
#ライブラリのインポート
library(dplyr)
# irisデータを使用
df <- iris
# 2クラス分類のためsetosaとversicolorを用いる
data <- df[df$Species == c("setosa","versicolor"),]
# クラスのカテゴリ変数を数値に変換
data <- mutate(data,
Species = case_when(
Species == "setosa" ~ "0"
, Species == "versicolor" ~ "1"
),
Species = as.numeric(Species)
)
# 数値データの標準化
data <- cbind(scale(data[1:4]), data[5])
# データフレームを行列に変換
data <- matrix(as.matrix(data), nrow(data), ncol(data))
# データの順番をシャッフル
set.seed(42)
shuffle <- sample(1:50)
data <- data[shuffle,]
# trainデータとtestデータの分類
split_data <- 1:40
x_train <- data[split_data, 1:4]
y_train <- data[split_data, 5]
x_test <- data[-split_data, 1:4]
y_test <- data[-split_data, 5]順番に説明すると、まずデータを使う部分(setosaとversicolor)だけ取り出します。次にこのラベルを数値(0,1)に変換します。学習しやすいようにデータを標準化し、データフレームのままでは計算できないので行列に変換します。最後にデータをシャッフルしtrainデータとtestデータに分けます。
3.2 関数の定義
3.2.1 sigmoid, log 関数の設定
まず最初に必要な関数を定義します。sigmoid関数は先ほど説明した通りです。また、交差エントロピー誤差関数を定義するときに logを用いますがこれも0以下の値だと定義できず、大きすぎると発散するため最大最小を設定しその値に収まるようにしています。コードは以下の通りです。
#sigmoidの定義
sigmoid <- function(x){
return (exp(pmin(x, 0)) / (1+exp(-abs(x))))
}
# logの発散を防ぐ
log_p <- function(x){
return (log(pmax(pmin(x, 1e-10), 1e+10)))
}3.2.2 重み, バイアスの設定
重み w とバイアス b は今回学習するバラメータです。このパラメータを勾配降下法を用いて更新していきます。はじめにこの重み w とバイアス b の設定をします。設定学習して更新されていきますので設定する値はなんでもいいです（設定によっては学習する回数が増えたりする場合はありますが…）。今回はw を正規分布、b を 0 に設定しました。
3.2.3 train関数, 評価関数(val)の設定
今回の設定で最も重要なところです。これまで求めてきた式を用いてcodeを書いていきます。式をまとめると
このようになります。
これをもとにcodeをかいていきます。
####### train関数 ######
train_logistic <- function(x, y, W, b, eps=0.1){
# 予想
y_p <- sweep(x %*% W, 2, b, FUN='+')
y_p <- sigmoid(y_p)
# 目的関数
cost <- mean(- y * log_p(y_p) - (1 - y) * log_p(1 - y_p))
delta <- y_p - y
#パラメータ更新
dW <- (t(x) %*% delta) / 50
db <- sum(delta) / 50
W <<- W - eps * dW
b <<- b * db
return (cost)
}
b は 1つの値なので sweep() を使って全ての要素に足し合わせています。t()は転置です。
w と b は更新する変数なので最初に設定した値にグローバル変数定義（<<-）を用いて更新します。
同様にval関数の code は以下のようになります。予測するだけなので cost は必要ありません。
######### val関数 #########
val_logistic <- function(x){
y_p <- sweep(x %*% W, 2, b, FUN='+')
y_p <- sigmoid(y_p)
return(y_p)
}3.3 学習と推論
長い設定が終わったので実際に学習を回していきます。今回は50エポック学習を回します。あまり多く回すと過学習を起こす可能性があります。過学習を回避する手法はミニバッチ学習やL1,L2正規化やk分割交差検証などたくさんありますが今回は特に行いませんでした。
では学習、推論、そして評価をしてみます。半分以上は上で書いた code なので付け加えたのは以下の code です。
#学習(50エポック)
epoch <- 50
for (i in 1:epoch){
cost <- train_logistic(x_train, y_train, W , b)
}
#推論
y_pred <- val_logistic(x_test)
for (i in 1:10){
if (y_pred[i] <= 0.5){
y_pred[i] <- 0
}else {
y_pred[i] <- 1
}
}予想した10個の結果をみてみます。
全て正解になったと思います。
4. まとめ
最も簡単なモデルであるロジスティック回帰を用いて2クラス分類をしてみました。出力を1つではなく増やすことで多クラス分類をすることができます。また、層を増やしていくことでMLPのとして扱うことができます。
数式を追ったりしたためとても長くなってしまいましたが、原理は少し理解できたと思います。お疲れ様でした。
