[Rの]Andrew先生の機械学習3-1[練習]

3回目前半はロジスティック回帰を使った文字認識

識別するのは1~10の数字。
20 X 20ピクセルの手書き数字の画像サンプルが5000個

講義の問題ではMATLABファイル形式で与えられるので、readMat関数で読む。

やりたいのは書かれている数字は1~10のどれなのか識別することなので、10クラスの識別問題となる。
こういう場合は1度に多クラス識別をするのでなく、1対多の識別モデルをクラスの数だけ作成し。最もあてはまるクラスに識別するのが常套手段らしい。

サンプルデータの可視化
実はロジスティック回帰の部分よりも可視化に手間がかかった。
20 X 20の画像データ が5000個あり、これが5000 X 400の配列として保存されているため、これを画像データとして復元している。

displayData <- function(X)
{
    N = 10
    m = length(X[, 1])
    n = length(X[1, ])
    width = sqrt(n)
    mnist = matrix(0, N * width, N * width)

    for (i in 1 : N) {
        for (j in 1 : N) {
            rIndexStart = width * (i - 1) + 1
            rIndexEnd = width * i
            cIndexStart = width * (j - 1) + 1
            cIndexEnd = width * j
            mnist[rIndexStart : rIndexEnd, cIndexStart : cIndexEnd] <- makeImage(X[floor(runif(1, min = 1, max = m + 1)), ])
        }
    }
    image(mnist)
}
makeImage <- function(X)
{
    n <- length(X)
    width <- sqrt(n)
    mnist <- matrix(0, width, width)
    for (i in 1 : width) {
        for (j in 1 : width) {
            mnist[i , width - j + 1] <- X[width * (i - 1) + j]
        }
    }
    return(mnist)
}


※Rでは":"より"*"や"+"の方が優先される。

ロジスティック回帰パラメータのトレーニング 
ロジスティック回帰のパラメータを、forループを回して10個分トレーニングしているだけで、やっていることは前回と変わらない。
ベクトル化しているので、できあがるのは10 X 401の行列。(画像データが400次元のベクトルで表現されており、それにバイアス項が加わり401次元)

oneVsAll <- function(X, y, num_labels, lambda)
{
    m = length(X[, 1])
    n = length(X[1, ])

    # add bias
    initial_theta = rep(0, n + 1)
    X = cbind(matrix(1, m, 1), X)
   
    all_theta = matrix(0, num_labels, n + 1)
    print(" ---started to train paramters---")
    for (i in 1 : num_labels)
    {   
        result <- optim(par = initial_theta, fn = costFunctionReg, gr = costGr,method = "BFGS" , X = X, y = y == i, lambda = lambda)
        all_theta[i, ] = result$par
        cat(i * 100 / num_labels, "% optimized\n")
    }
    return(all_theta)
}
※コスト関数は前回と同一

訓練したパラメータでロジスティック回帰
インプットに対して回帰を行い、10個のモデルのうち最も適合したインデックスを返す。
必要であれば閾値を設けて、すべてのモデルで閾値を下回ったら エラーか何かを返すこともできる。

predictOneVsAll <- function(all_theta, X)
{
    m <- length(X[, 1])
    num_labels <- length(all_theta[, 1])
    X <- cbind(matrix(1, m, 1), X)

    predicts <- matrix(0, m, num_labels)
    for (i in 1 : num_labels)
    {
        predicts[, i] <- X %*% t(all_theta)[, i]
    }

    p <- matrix(0, m, 1)
    for (i in 1 : m)
    {
        p[i, 1] <- which.max(predicts[i, ])
    }
    return(p)
}




実行してみる
ex3 <- function()
{
    require(R.matlab)
    datamat <- readMat('ex3data1.mat')
    X <- datamat$X
    y <- datama$y
    displayData(X)

    lambda = 0.1
    num_labels = 10
    all_theta = oneVsAll(X, y, num_labels, lambda)
    pred = predictOneVsAll(all_theta, X)
    print("Training Set Accuracy:")
    print(mean(pred == y) * 100)
}

結果
> ex3()
[1] " ---started to train paramters---"
10 % optimized
20 % optimized
30 % optimized
40 % optimized
50 % optimized
60 % optimized
70 % optimized
80 % optimized
90 % optimized
100 % optimized
[1] "Training Set Accuracy:"
[1] 95.04

なぜかMATLABで書いたものより識別率が上がった。
講義を受けているときには、このあたりで機械学習している実感が出てきてテンションあがった。