次に、例題を作って、手順を理解してみる。
5教科の成績から入学試験の合否を予測する
5教科の成績から入学試験の合否を予測するロジスティック回帰モデルを考えてみる。
ここでは、5教科の成績を入力変数、入学試験の合否を出力変数にして、ロジスティック回帰モデルのパラメータを、重み（係数）$w$とバイアス（切片）$b$からなる$\theta$とする。 そして、入力変数を反復学習して、ロジットを計算し、シグモイド関数を適用して確率を求め、コスト関数からパラメータを推定する。
データセットの準備
20人分のデータセットを作成して、最初の15人分をロジスティック回帰モデルの構築に使う学習セット、最後の5人分を評価用いる新しいデータセットにする。
このとき、5教科の成績を0から100点からなる入力変数、入学試験の合否を合格なら1、不合格なら0にした出力変数にする。なお、実際には、データのスケーリングや正規化などの前処理も必要ですが、ここでは省略する。
| データセット | 英 | 数 | 国 | 理 | 社 | 合否 | |
|---|---|---|---|---|---|---|---|
| 1 | 学習 | 75 | 82 | 90 | 68 | 78 | 1 |
| 2 | 学習 | 92 | 88 | 77 | 80 | 85 | 1 |
| 3 | 学習 | 63 | 70 | 75 | 58 | 65 | 0 |
| 4 | 学習 | 81 | 78 | 86 | 92 | 88 | 1 |
| 5 | 学習 | 70 | 65 | 72 | 68 | 74 | 0 |
| 6 | 学習 | 95 | 92 | 88 | 97 | 90 | 1 |
| 7 | 学習 | 72 | 75 | 68 | 64 | 70 | 0 |
| 8 | 学習 | 88 | 92 | 84 | 78 | 80 | 1 |
| 9 | 学習 | 60 | 58 | 63 | 52 | 55 | 0 |
| 10 | 学習 | 78 | 80 | 76 | 82 | 85 | 1 |
| 11 | 学習 | 68 | 72 | 64 | 58 | 62 | 0 |
| 12 | 学習 | 83 | 87 | 92 | 89 | 94 | 1 |
| 13 | 学習 | 76 | 70 | 72 | 78 | 80 | 1 |
| 14 | 学習 | 55 | 62 | 58 | 50 | 60 | 0 |
| 15 | 学習 | 90 | 85 | 92 | 88 | 94 | 1 |
| 16 | 評価 | 72 | 74 | 68 | 76 | 80 | 1 |
| 17 | 評価 | 62 | 58 | 65 | 60 | 64 | 0 |
| 18 | 評価 | 80 | 82 | 85 | 78 | 84 | 1 |
| 19 | 評価 | 67 | 72 | 66 | 58 | 70 | 0 |
| 20 | 評価 | 94 | 90 | 88 | 92 | 96 | 1 |
X = c( 75, 82, 90, 68, 78, 1, 92, 88, 77, 80, 85, 1, 63, 70, 75, 58, 65, 0, 81, 78, 86, 92, 88, 1, 70, 65, 72, 68, 74, 0, 95, 92, 88, 97, 90, 1, 72, 75, 68, 64, 70, 0, 88, 92, 84, 78, 80, 1, 60, 58, 63, 52, 55, 0, 78, 80, 76, 82, 85, 1, 68, 72, 64, 58, 62, 0, 83, 87, 92, 89, 94, 1, 76, 70, 72, 78, 80, 1, 55, 62, 58, 50, 60, 0, 90, 85, 92, 88, 94, 1, 72, 74, 68, 76, 80, 1, 62, 58, 65, 60, 64, 0, 80, 82, 85, 78, 84, 1, 94, 90, 88, 92, 96, 1 ) X <- data.frame(matrix(X, ncol = 6, byrow = TRUE)) colnames(X) = c("英", "数", "国", "理", "社", "合否") y <- factor(X$"合否") X <- X[, 1:5] new_X <- X[-c(1:15), ] new_y <- y[-c(1:15)] X <- X[1:15, ] y <- y[1:15]
ロジスティック回帰モデルの考え方
回帰を考えるとき、予測対象が合格（1）、不合格（0）1と0の二値だと扱いにくい。予測対象を合格確率と考えて、テストの点数が高いほど1に近づき、低いほど0に近づくと考えます。
線形予測子を作る
まず、入力変数$x$を5教科の点数、$P(x)$を合格確率と考えます。
このとき、5教科の点数を個別に見るよりも、一つにまとめて、その大小で合格確率を考える方が楽です。そこで、5教科の点数をまとめる式$z$を考えます。
式$z$は係数が$w$、切片が$b$になる多項式で、出力変数$y$の予測に用いる線形予測子になります。
$$ z = w^{T}x + b $$
ここで、$w$と$b$は、ロジスティック回帰モデルのパラメータ$\theta$の重みとバイアスになります。
ロジット関数を作る
次に、出力変数$y$は合格確率なので範囲に上下限があります。これは使いにくいので、オッズ比の対数をとるロジット関数を考えます。
$$ logit(y)=log \frac{y}{1-y} $$
$logit(y)$の予測にすることで、予測範囲は$- \infty \sim \infty$になります。
リンク関数を設定する
次に、多項式$z$が出力変数$y$のロジット関数$logit(y)$になると考えます。このとき、ロジット関数は$x$と$y$の関係を示すリンク関数になります。
$$ logit(y) = z $$
これを変換すると、ロジスティック関数（シグモイド関数）になります。
$$ y = \frac{1}{1+e^{(w^{T}x + b)}} $$
パラメータ推定、学習、誤差計算
実計算では、パラメータ$\theta$を繰り返し修正しながら、作成したロジスティック回帰モデルの予測値と実測値との乖離を最小化して行きます。
このパラメータ$\theta$の修正は、最尤推定法でパラメータ$\theta$を推定する、予測値と実測値のクロスエントロピー誤差を求める、勾配降下法でこの誤差を最小化する流れで進みます。
# ロジット関数の定義 logits <- function(X, w, b) X %*% w + b # シグモイド関数の定義 sigmoid <- function(z) 1 / (1 + exp(-z)) # クロスエントロピー誤差（コスト関数）の定義 cost <- function(y, y_pred) { sum(y * log(y_pred) + (1 - y) * log(1 - y_pred)) }
この定義では、sigmoid(logit(X, w, b))が予測値の計算になります。
また、クロスエントロピー誤差のコスト関数はyが実際の合否、y_predが合否の予測値で、y_predがyに近いほど小さくなります。まず、y * log(y_pred)で合格（1）の予測確率の対数を求め、次に、(1 - y) * log(1 - y_pred)で不合格（0）の予測確率の対数を求めます。次に、この2つを足し合わせ、最後に合計して全体誤差を求めています。
# パラメータの初期化 w <- matrix(0, nrow = ncol(X), ncol = 1) b = matrix(0, nrow = 1, ncol = 1) # ハイパーパラメータの設定 learning_rate <- 0.01 num_iterations <- 100 # 最尤推定法によるパラメータの最適化 cost_iter = c() for (i in 1:num_iterations) { # 予測値の計算 y_pred <- sigmoid(logit(X, w, b)) # 勾配の計算 dw <- (1 / length(x)) * sum((y_pred - y) * x) db <- (1 / length(x)) * sum(y_pred - y) # パラメータの更新 w <- w + learning_rate * dw b <- b + learning_rate * db cost_iter = c(cost_iter, cost(y, y_pred)) }
ここでは、勾配降下法の反復計算でパラメータ$\theta$を最適化しています。
まず、予測値y_predを求めて、実際の合否yとの差分から勾配を計算します。次に、この勾配に学習率learning_rateを掛け合わせて、パラメータを更新する程度を決めます。最後に、これを重み$w$とバイアス$b$に加えて、パラメータを更新します。指定の回数num_iterationsの反復を行って、最後に得た重み$w$とバイアス$b$が最適化したパラメータ$\theta$と考えます。
このとき、反復計算の各回のクロスエントロピー誤差をcost_iterに貯めています。これは、計算後に誤差の値の減少を確認するためです。この値が減少していなければ、学習率や反復回数の見直しが必要です。
また、誤差の減少が重要であれば、固定値でなく、誤差が指定の値以下になるまで反復計算する方法もあります。ですが、過度な反復計算は過学習（overfitting）のリスクが上げ、未知のデータに対する汎化性能を下げてしまいます。ここには注意が必要です。
一般的には、反復計算の回数を適切な設定の他、正則化（regularization）や交差検証（cross-validation）などの手法を使うことになるでしょう。
予測
学習が終了した$\theta$を使用して新しいデータを予測します。ロジスティック回帰モデルは確率を出力するため、0.5以上なら1と判定し、0.5未満なら0と判定します。
# 正解率、適合率、再現率、F1スコアを計算 pred_score = function(predictions, true_labels){ # 正解率（Accuracy）の計算 accuracy <- sum(predictions == true_labels) / length(predictions) # 適合率（Precision）と再現率（Recall）の計算 TP <- sum(predictions == 1 & true_labels == 1) # True Positive FP <- sum(predictions == 1 & true_labels == 0) # False Positive FN <- sum(predictions == 0 & true_labels == 1) # False Negative precision <- TP / (TP + FP) recall <- TP / (TP + FN) # F1スコア（F1 Score）の計算 f1_score <- 2 * (precision * recall) / (precision + recall) score = list(accuracy, precision, recall, f1_score) return(score) } # 予測結果と真の値のサンプルデータ（0: Negative, 1: Positive） pred_y = sigmoid(logit(new_X, w, b)) pred_y <- ifelse(pred_y >= .5, 1, 0) pred_score(prediction=pred_y, true_labels=new_y)
