今回はk近傍法 (k-nearest neighbor algorithm; kNN)を使った分類を紹介します。今回も前回と同様、なるべく外部パッケージを使わずkNNを実装してみたいと思います1。
仮想データの作成
今回も仮想データの作成から始めましょう。やはり楽しいのは実際のデータをいじるのですが、既知のデータ生成過程から作成したデータセットは「正解」が分かっているので、シミュレーションなどに向いています。
今回も結果変数が(0, 1)の二値であるとし、以下のデータ生成過程にしたがうとします。
\[\begin{eqnarray} Y & = & \begin{cases} 1, & \text{if $y^* \geq 30$}, \\ 0, & \textrm{otherwise}, \end{cases}\\ y^* & \sim & 1 - 0.8 X1 + 2 X1^2 - 3 X2 + 0.2 X2^2 + \varepsilon, \nonumber \\ \varepsilon & \sim & \textrm{Normal}(0, 10) \nonumber \\ X1, X2 & \sim & \textrm{Uniform}(-5, 5) \nonumber \end{eqnarray}\]それでは早速データセットを作成し、dfと名付けましょう。
set.seed(19861008)
x1 <- runif(100, -5, 5)
x2 <- runif(100, -5, 5)
e <- rnorm(100, 0, 10)
ya <- 1 - (0.8 * x1) + (2 * x1^2) - 3 * x2 + (0.2 * x2^2) + e
y <- ifelse(ya >= 30, 1, 0)
df <- data.frame(x1 = x1, x2 = x2, y = y)
rm(x1, x2, e, ya, y)それでは、説明変数 (\(X1\)と\(X2\))と結果変数 (\(Y\))の関係を確認しましょう。
前回の例のように一本の線で2つを分けるのはなかなか難しいようです。
kNN法について
前回のロジスティック回帰は予め関数型が決められており、数個のパラメーター2だけを推定するだけで済みましたよね。このようなモデルをパラメトリック・モデル (parametric model)と呼びます。今回のkNNはノンパラメトリック・モデル (nonparametric model)の代表的なモデルです。推定するパラメーターの数も指定せず、パラメトリック・モデルに比べ計算量が多いです。しかし、回帰分析のような線形を前提としていないため、非線形の関係などにも柔軟に対応できるというメリットがあります。
kNN法ある地点がどのカテゴリーに属するかを決めるアルゴリズムです。どのカテゴリーに属するかは、ある点を中心にした最も近いk個の点を用いて多数決で決めます。ちなみに、距離はユークリッド距離を用います。たとえば、\(k = 3\)にした場合、(-3.5, -0.75)の点はどのカテゴリーに属するでしょうか。
ある点、(-3.5, -0.75)の近傍にある5つの点は赤 (1)が2つ、青 (0)が1つです。つまり、(-3.5, -0.75)は\(Y = 1\)の領域であることを意味します。このような作業をプロット全体に対して行います。プロット内の感覚を細かく刻めば刻むほど、より正確に計算できるようになります。
もし、\(k = 5\)ならどうでしょうか。
今回は赤 (1)が2つ、青 (0)が3つとなり、(-3.5, -0.75)は\(Y = 0\)の領域であるという結果が得られました。このようにkNNの性能を最大化するためには、適切なkを指定する必要があります。そのためには複数のkに対して分類機を生成し、その性能を評価・比較する必要があります。これについては後半で説明します。
kNN分類器の実装
まずは、二点間の距離を求める関数を作成します。ユークリッド距離なので関数名はCalc_EDと名付けましょう。必要な引数は二点の座標、(x, y)と(p, q)です。ユークリッド距離は以下のように計算します。
Calc_ED <- function(x, y, p, q){
return(sqrt((x - p)^2 + (y - q)^2))
}つづいて、任意の点 (p, q)の近傍にあるk個の点を抽出し、多数決を行う関数 (CountK)を作成します。この関数は以下のような手順で動きます。
- 距離を格納するベクトル、
D.vecを用意する。 - データ内のすべての点、(x1, x2)と(p, q)の計算し、
D.vecに格納する。 - データを距離が短い順で並び替える。
- データをk行まで切り取る。
- 残ったデータ内の\(Y\)の平均が0.5以上なら1を、未満なら0を返す。
CountK <- function(df, p, q, k){
D.vec <- rep(NA, nrow(df)) # Step 1
for (i in 1:nrow(df)) {
D.vec[i] <- Calc_ED(df$x1[i], df$x2[i], p, q) # Step 2
}
df <- cbind(df, ED = D.vec)
df <- df[order(df$ED), ] # Step 3
df <- df[1:k, ] # Step 4
return(list(df = df,
result = ifelse(mean(df$y) >= 0.5, 1, 0))) # Step 5
}実際にやってみましょう。先ほどの例、(-3.5, -0.75)でそれぞれ\(k = 3\)、\(k = 5\)で試してみます。
Test.K3 <- CountK(df = df, p = -3.5, q = -0.75, k = 3)
print(Test.K3)## $df
## x1 x2 y ED
## 79 -3.382601 -1.395158 1 0.6557525
## 84 -2.838075 -1.053057 0 0.7280032
## 89 -3.911556 -1.612889 1 0.9560103
##
## $result
## [1] 1(-3.5, -0.75)に最も近い3個の点は(-3.383, -1.395)、(-2.838, -1.053)、(-3.912, -1.613)であり、それぞれ(-3.5, -0.75)からの距離は0.656、0.728、0.956です。この中で\(Y = 1\)のケースは2個であり、多数決の結果、(-3.5, -0.75)は\(Y = 1\)に属することが分かりますね。
Test.K5 <- CountK(df = df, p = -3.5, q = -0.75, k = 5)
print(Test.K5)## $df
## x1 x2 y ED
## 79 -3.382601 -1.39515795 1 0.6557525
## 84 -2.838075 -1.05305692 0 0.7280032
## 89 -3.911556 -1.61288883 1 0.9560103
## 16 -3.566713 0.37596630 0 1.1279409
## 45 -4.515107 0.06657151 0 1.3027784
##
## $result
## [1] 0続いて、\(k = 5\)の場合、(-3.5, -0.75)に最も近い3個の点は(-3.383, -1.395)、(-2.838, -1.053)、(-3.912, -1.613)、(-3.567, 0.376)、(-4.515, 0.067)であり、それぞれ(-3.5, -0.75)からの距離は0.656、0.728、0.956、1.128、1.303です。この中で\(Y = 1\)のケースは2個であり、多数決の結果、(-3.5, -0.75)は\(Y = 0\)に属することが分かりますね。
この手順はプロット内のすべての点において実施します。\(X1\)と\(X2\)の範囲はそれぞれ(-5, 5)ですね。-5から5まで0.001刻みでやってもいいですが、あまりにも非効率的です。今回は、0.1刻みでやってみましょう。まずは、結果を格納するマトリックス、Contour.Mat3を作成します。
# -5から5まで0.1に刻みます
x1.vec <- seq(-5, 5, 0.1)
x2.vec <- seq(-5, 5, 0.1)
# 予測結果を格納する行列
Contour.Mat3 <- matrix(rep(NA, length(x1.vec) * length(x2.vec)),
nrow = length(x1.vec))続いて、すべてのx1.vec、x2.vecの値に対してどの領域に属するかを計算し、Contour.Mat3に格納します。まずは\(k = 3\)の設定でします。
for (i in 1:length(x1.vec)) { # すべてのx1.vecに対して
for (j in 1:length(x2.vec)) { # すべてのx2.vecに対して
# k = 3のkNNを実行
Contour.Mat3[i, j] <- CountK(df = df,
p = x1.vec[i],
q = x2.vec[j],
k = 3)$result
}
}いよいよ楽しい可視化の時間です。まずは、背景を予測した\(Y\)の値に応じて、赤と青に塗ります。その上に、実際のデータの散布図をオーバラップし、境界線を引くために等高線図を乗せます。これは前回のロジスティック回帰分類器の可視化と同じ手順です。ただし、今回は確率ではなく、0か1かなので背景がグラデーションではなく、二色のみとなります。等高線も0と1の境界線になりますね。
image(x = x1.vec, y = x2.vec, z = Contour.Mat3,
col = c(rgb(0, 0, 1, 0.5), rgb(1, 0, 0, 0.5)),
xlab = "X1", ylab = "X2")
points(x = df$x1, y = df$x2, pch = 21,
col = "black", bg = c("blue", "red")[df$y + 1])
contour(x = x1.vec, y = x2.vec, z = Contour.Mat3,
levels = 1, labels = "", add = TRUE)
legend("topright", pch = 19, col = c("blue", "red"),
legend = c("Y = 0", "Y = 1"))
いくつか誤分類された点がありすが、全体的に悪くない性能です。それでは\(k = 5\)にしてみましょう。\(k\)を奇数にするのは多数決の際に同率の可能性を排除するためです。
Contour.Mat5 <- matrix(rep(NA, length(x1.vec) * length(x2.vec)),
nrow = length(x1.vec))
for (i in 1:length(x1.vec)) {
for (j in 1:length(x2.vec)) {
Contour.Mat5[i, j] <- CountK(df = df,
p = x1.vec[i],
q = x2.vec[j],
k = 5)$result
}
}image(x = x1.vec, y = x2.vec, z = Contour.Mat5,
col = c(rgb(0, 0, 1, 0.5), rgb(1, 0, 0, 0.5)),
xlab = "X1", ylab = "X2")
points(x = df$x1, y = df$x2, pch = 21,
col = "black", bg = c("blue", "red")[df$y + 1])
contour(x = x1.vec, y = x2.vec, z = Contour.Mat5,
levels = 1, labels = "", add = TRUE)
legend("topright", pch = 19, col = c("blue", "red"),
legend = c("Y = 0", "Y = 1"))
\(k = 3\)の時に比べ、誤分類されているケースが多いようです。これを見ると\(k\)が小さいほど、分類器の性能は上がるようです。それでは\(k = 1\)にしてみましょう。
Contour.Mat1 <- matrix(rep(NA, length(x1.vec) * length(x2.vec)),
nrow = length(x1.vec))
for (i in 1:length(x1.vec)) {
for (j in 1:length(x2.vec)) {
Contour.Mat1[i, j] <- CountK(df = df,
p = x1.vec[i],
q = x2.vec[j],
k = 1)$result
}
}image(x = x1.vec, y = x2.vec, z = Contour.Mat1,
col = c(rgb(0, 0, 1, 0.5), rgb(1, 0, 0, 0.5)),
xlab = "X1", ylab = "X2")
points(x = df$x1, y = df$x2, pch = 21,
col = "black", bg = c("blue", "red")[df$y + 1])
contour(x = x1.vec, y = x2.vec, z = Contour.Mat1,
levels = 1, labels = "", add = TRUE)
legend("topright", pch = 19, col = c("blue", "red"),
legend = c("Y = 0", "Y = 1"))
なんと、予測精度が100%に達しました！背景の色と点の色がすべて一致していますね。
kNNは\(k = 1\)の時に最強でした。めでたしめでたしー
モデル評価によるkの選択
と言いたいわけですが、実は違います。kNNの場合、\(k\)が小さいほどバイアスが小さくなりますが、分散が大きくなります。バイアスは実測値と予測値の乖離であり、\(k = 1\)の場合、実測値と予測値は一致しているため、バイアスがない状態です。分散はモデルの安定性です。ケースに変動があった場合、\(k = 5\)の分類器は\(k = 1\)の分類器より変動が小さいです。\(k\)が小さいほど、縮小不可能な誤差 (irreducible error)に敏感となり、モデルが不安定になります。詳細はJames, Witten, Hastie, and Tibshirani (2013)3の第2章を参照して下さい。
どのモデルが良いモデルか、今回は以下の指標で確認してみたいと思います。
- エラー率
- \(\kappa\)統計量
- 4-fold Cross Validation
テストするモデルは\(k = \{1, 3, 5, 7, 9\}\)の5つのモデルです。
エラー率
エラー率は1から的中率を引いたものですね。的中率の計算方法はこれまでとほぼ同様です。ただし、先はx1.vecとx2.vecの値に応じて予測値を出すのではなく、データの値に応じて予測値も出すという点ですね。早速やってみましょう。
まずは、予測値を格納する空のベクトルを作成し、そこにdf内のすべての点に対してkNNを実行します。
K1.Pred <- rep(NA, nrow(df))
K3.Pred <- rep(NA, nrow(df))
K5.Pred <- rep(NA, nrow(df))
K7.Pred <- rep(NA, nrow(df))
K9.Pred <- rep(NA, nrow(df))
for (i in 1:nrow(df)) {
K1.Pred[i] <- CountK(df, df$x1[i], df$x2[i], k = 1)$result
K3.Pred[i] <- CountK(df, df$x1[i], df$x2[i], k = 3)$result
K5.Pred[i] <- CountK(df, df$x1[i], df$x2[i], k = 5)$result
K7.Pred[i] <- CountK(df, df$x1[i], df$x2[i], k = 7)$result
K9.Pred[i] <- CountK(df, df$x1[i], df$x2[i], k = 9)$result
}それでは各分類機のエラー率を見てみましょう。今後、\(\kappa\)統計量の計算にも使いたいので今のうちに混同行列も作っておきます。
ConMat.K1 <- table(df$y, K1.Pred)
ConMat.K3 <- table(df$y, K3.Pred)
ConMat.K5 <- table(df$y, K5.Pred)
ConMat.K7 <- table(df$y, K7.Pred)
ConMat.K9 <- table(df$y, K9.Pred)一個一個行列を見るのは面倒くさいので、エラー率の折れ線グラフで出します。
Accuracy.Rate <- c((sum(diag(ConMat.K1)) / sum(ConMat.K1)),
(sum(diag(ConMat.K3)) / sum(ConMat.K3)),
(sum(diag(ConMat.K5)) / sum(ConMat.K5)),
(sum(diag(ConMat.K7)) / sum(ConMat.K7)),
(sum(diag(ConMat.K9)) / sum(ConMat.K9)))
plot(x = c(1, 3, 5, 7, 9), y = 1 - Accuracy.Rate, type = "b",
xlab = "k", ylab = "Error Rate", xaxt = "n")
axis(1, at = c(1, 3, 5, 7, 9), labels = c(1, 3, 5, 7, 9))
\(k\)が大きくなるほどエラー率が上がりますね。つまり、ここまでの結果を見るとベストモデルは\(k = 1\)のkNNになります。
\(\kappa\)統計量
次は\(\kappa\)統計量です。\(\kappa\)統計量の関数は前回作成したCalc_Kappaをそのまま使いましょう。
Calc_Kappa <- function(tab){
p0 <- sum(diag(tab)) / sum(tab) # 的中率のことです
pe1 <- (sum(tab[, 2]) / sum(tab)) * (sum(tab[2, ]) / sum(tab))
pe2 <- (sum(tab[, 1]) / sum(tab)) * (sum(tab[1, ]) / sum(tab))
pe <- pe1 + pe2
k <- (p0 - pe) / (1 - pe)
return(round(k, 3))
}今回も折れ線グラフです。
Kappa.Vec <- c(Calc_Kappa(ConMat.K1),
Calc_Kappa(ConMat.K3),
Calc_Kappa(ConMat.K5),
Calc_Kappa(ConMat.K7),
Calc_Kappa(ConMat.K9))
plot(x = c(1, 3, 5, 7, 9), y = Kappa.Vec, type = "b",
xlab = "k", ylab = expression(paste(kappa, " statistics")),
xaxt = "n")
axis(1, at = c(1, 3, 5, 7, 9), labels = c(1, 3, 5, 7, 9))
ここでも、\(k\)が大きくなるほど\(\kappa\)統計量は小さくなります。やはりこの結果を見るとベストモデルは\(k = 1\)のkNNになります。
4-fold Cross Validation
エラー率と\(\kappa\)統計量を見る限り、\(k = 1\)のkNNが最強のモデルであることは疑いの余地がないように見えます。しかし、今回は一味違うと思います。CVの出番です。k-fold CVのkは一般的に5か10が使われますが、今回は可視化の便宜上、\(k = 4\)にしました。
手順はロジスティック分類器の回と同じです。まずは、\(k = 1\)の時の4-fold CVからです。今回はCVの結果の可視化まで含めて示したいと思います。単純に4-fold CVのエラー率を計算したいだけなら、ここは飛ばして後半へ進んで下さい。
# まずデータの順番をランダム化します
k.index <- sample(1:nrow(df), nrow(df), replace = FALSE)
# 次は、25個ずつのデータに分割します。
k1.index <- k.index[1:25]
k2.index <- k.index[26:50]
k3.index <- k.index[51:75]
k4.index <- k.index[76:100]
# 4-fold CVなので行列は4つ用意します。
CV.K1_1 <- matrix(rep(NA, length(x1.vec) * length(x2.vec)),
nrow = length(x1.vec))
CV.K1_2 <- matrix(rep(NA, length(x1.vec) * length(x2.vec)),
nrow = length(x1.vec))
CV.K1_3 <- matrix(rep(NA, length(x1.vec) * length(x2.vec)),
nrow = length(x1.vec))
CV.K1_4 <- matrix(rep(NA, length(x1.vec) * length(x2.vec)),
nrow = length(x1.vec))
# Training setのみを用いてkNNを実行します。
# [-k1.index]にするとk1.indexに指定された行を除外できます。
for (i in 1:length(x1.vec)) {
for (j in 1:length(x2.vec)) {
CV.K1_1[i, j] <- CountK(df = df[-k1.index, ],
p = x1.vec[i], q = x2.vec[j],
k = 1)$result
CV.K1_2[i, j] <- CountK(df = df[-k2.index, ],
p = x1.vec[i], q = x2.vec[j],
k = 1)$result
CV.K1_3[i, j] <- CountK(df = df[-k3.index, ],
p = x1.vec[i], q = x2.vec[j],
k = 1)$result
CV.K1_4[i, j] <- CountK(df = df[-k4.index, ],
p = x1.vec[i], q = x2.vec[j],
k = 1)$result
}
}背景色はkNN分類機の予測値を、点は分類機の作成に使わなかった点のみを表示させます。
# 2 x 2の図を載せるます
par(mfrow = c(2, 2))
# これまでのプロットと同じですが、散布図はTest setのみとなります。
image(x = x1.vec, y = x2.vec, z = CV.K1_1,
col = c(rgb(0, 0, 1, 0.5), rgb(1, 0, 0, 0.5)),
xlab = "X1", ylab = "X2", main = "1st folder")
points(x = df$x1[k1.index], y = df$x2[k1.index], pch = 21,
col = "black", bg = c("blue", "red")[df$y[k1.index] + 1])
contour(x = x1.vec, y = x2.vec, z = CV.K1_1,
levels = 1, labels = "", add = TRUE)
image(x = x1.vec, y = x2.vec, z = CV.K1_2,
col = c(rgb(0, 0, 1, 0.5), rgb(1, 0, 0, 0.5)),
xlab = "X1", ylab = "X2", main = "2nd folder")
points(x = df$x1[k2.index], y = df$x2[k2.index], pch = 21,
col = "black", bg = c("blue", "red")[df$y[k2.index] + 1])
contour(x = x1.vec, y = x2.vec, z = CV.K1_2,
levels = 1, labels = "", add = TRUE)
image(x = x1.vec, y = x2.vec, z = CV.K1_3,
col = c(rgb(0, 0, 1, 0.5), rgb(1, 0, 0, 0.5)),
xlab = "X1", ylab = "X2", main = "3rd folder")
points(x = df$x1[k3.index], y = df$x2[k3.index], pch = 21,
col = "black", bg = c("blue", "red")[df$y[k3.index] + 1])
contour(x = x1.vec, y = x2.vec, z = CV.K1_3,
levels = 1, labels = "", add = TRUE)
image(x = x1.vec, y = x2.vec, z = CV.K1_4,
col = c(rgb(0, 0, 1, 0.5), rgb(1, 0, 0, 0.5)),
xlab = "X1", ylab = "X2", main = "4th folder")
points(x = df$x1[k4.index], y = df$x2[k4.index], pch = 21,
col = "black", bg = c("blue", "red")[df$y[k4.index] + 1])
contour(x = x1.vec, y = x2.vec, z = CV.K1_4,
levels = 1, labels = "", add = TRUE)
par(mfrow = c(1, 1))先ほどはエラー率0の最強分類器でしたが、今回はそこそこ誤分類が見られます。また、色塗りのパターンにかなりの変動が見られますね。
以下ではコードを見せずに、全く同じことを\(k = \{3, 5, 7, 9\}\)の場合に対して繰り返します。
まずは\(k = 3\)です。
つづいて、\(k = 5\)です。
続いて、\(k = 7\)です。
最後に、\(k = 9\)です。
いかがですか。CVで検証してみると必ずしも小さい\(k\)が高い予測精度を保証しないことが分かります。また、\(k\)が大きいほど、予測値 (プロットの色塗り)が安定していることが分かります。つまり、4つのプロットのいの塗パターンが大きく変わらないことを意味します。これが「モデルの分散が小さい」ことの意味です。
それでは、実際、どのモデルが最も予測精度が高いかを、数値で確認してみましょう。それぞれのモデルに対して、(4-fold CVなので)4つのエラー率を計算し、その平均値が最も低いモデルが良いモデルとします。
# まずは結果を格納するデータフレームから
# 各行は分類器、列はCVのIDです
CV.df <- data.frame(Fold1 = rep(NA, 5),
Fold2 = rep(NA, 5),
Fold3 = rep(NA, 5),
Fold4 = rep(NA, 5))
# k-fold CV関数の作成
# Test setを指定するindexを引数にします
kFold <- function(df, k, index) {
# Test setの予測値を格納する空ベクトル
temp.Pred <- rep(NA, length(index))
for (i in 1:nrow(df)) {
# indexに指定された点を除外して、
# indexに指定された点の予測を行います
temp.Pred[i] <- CountK(df[-index, ],
df$x1[index][i], df$x2[index][i],
k = k)$result
}
# 実際のyと予測値を比較し、エラー率を計算
temp.Error <- 1 - mean(temp.Pred == df$y[index])
return(temp.Error)
}
# 4-fold CVを実行します。
for (i in 1:5) { # Model index
CV.df[i, 1] <- kFold(df = df, k = c(1, 3, 5, 7, 8)[i],
index = k1.index)
CV.df[i, 2] <- kFold(df = df, k = c(1, 3, 5, 7, 8)[i],
index = k2.index)
CV.df[i, 3] <- kFold(df = df, k = c(1, 3, 5, 7, 8)[i],
index = k3.index)
CV.df[i, 4] <- kFold(df = df, k = c(1, 3, 5, 7, 8)[i],
index = k4.index)
}次は、横軸を\(k\)、縦軸を平均エラー率にした折れ線グラフを表示させます。データフレームの行単位平均はrowMeans()関数で計算できます4。
plot(x = c(1, 3, 5, 7, 9), y = rowMeans(CV.df),
type = "b", lwd = 2,
xaxt = "n", xlab = "k", ylab = "Error Rate")
axis(1, at = c(1, 3, 5, 7, 9), labels = c(1, 3, 5, 7, 9))
最もエラー率の低い分類器は\(k = 5\)のkNNアルゴリズムでした (エラー率は0.3)。エラー率が最も高いのは、これまで最強のモデルだった\(k = 1\)のkNNで、エラー率は約0.318です。最強のモデルかと思ったら最凶のモデルでした。エラー率は約0.02ほど改善され、これは\(k = 3\)の分類器が\(k = 1\)より約2%ほど正確に予測を行うことを意味します。今回の例で見ると、大差はないかも知れませんが、分類機の評価にCVが大いに役立つ5ことは分かって頂けたかなと思います。
kNN法のパッケージがいくつかあります。たとえば
classパッケージのknn()があります。↩前回の例だと、\(\beta_{i \in \{0, 1, 2\}}\)の3個だけです。↩
James, Gareth, Daniela Witten, Trevor Hastie, Robert Tibshirani. 2013. An Introduction to Statistical Learning. Springer.↩
apply()関数でも可能です。apply(CV.df, 1, mean)でも同じ結果が得られます。↩今回の結果はあくまでも4-fold CVの結果です。LOO CV、または10-fold CVでは\(k = 5\)以外のモデルが最良モデルである可能性も十分あり得ます。CVにおいて最適kはありません。あえていえば、k-fold CVでkは大きいほど良いです。パソコンのパワーが許容する範囲内でkを設定すれば良いでしょう。↩