[Rの]Andrew先生の機械学習6[練習]

5回目は飛ばして6回目。講義はex8まであるけども、このシリーズはこれで終わりにする。
今回はSVMを使った分類問題。
対象のデータは画像のような二次元特徴ベクトル + ラベル

線形分離不可能なので、SVMのカーネルトリックを使う。
今回のように線形分離できないときは、特徴ベクトルを非線形写像して線形分離可能な問題に変換する。
このあたりははじパタを読んで復習した。

実装に当たっては、非線形写像する関数でなく、カーネル関数と呼ばれる関数を定義してやれば十分。
Andrew先生の講義ではガウシアンカーネルを使ったが、今回はRBFカーネルを使う。（使うライブラリのデフォルトがRBFだからであって深い理由はない）

SVMの実装は面倒だったのでe1071というRのライブラリを使う。

• SVMの設定（cost）

SVMの中でもC-SVMというやつなので、costパラメータを設定する必要がある。
costはニューラルネットワークでいう正則化項のように、過学習を抑制するパラメータ。
ただしニューラルネットワークとは逆で、costを大きくするほど学習データに対して識別性能が良くなる代わりに過学習しがちになる。つまり、識別境界に複雑な形を許すようになる。
costを小さくすると、パラメータの上限が抑えられ、識別境界がシンプルな形になる。

• SVMの設定（gamma）

RBFカーネルはパラメータgamma(はじパタではσ)を持ち、こいつもチューニングしてやる必要がある。
こいつも大きいほど複雑な識別境界を描けるようになる。

R
※ドルマークが正しく表示されないらしい
ex6<- function()
{
    require("e1071")
    require(R.matlab)
    datamat <- readMat('ex6data2.mat')
    testdata <- data.frame(cbind(datamat$y, datamat$X))
    names(testdata) <- c("y", "x1","x2")
    testdata$y <- as.factor(testdata$y)
    plotData(datamat$X, datamat$y)
    par(new = T)

    result <- svm(y ~ . , data = testdata, cost = 10, gamma = 1)
   
    summary(result)
    print(result)

    px <- seq(0.9 * min(testdata$x1), 1.1 * max(testdata$x1),length= 100)
    py <- seq(0.9 * min(testdata$x2), 1.1 * max(testdata$x2),length= 100)
    pgrid <- expand.grid(px,py)
    names(pgrid)<-c("x1","x2")
    result.plot <- predict(result, pgrid, type="vector")
    contour(px, py, array(result.plot, dim=c(length(px),length(py))),col = "red", lwd=3, drawlabels=F, c(0.9 * min(testdata$x1), 1.1 * max(testdata$x1)), ylim = c(0.9 * min(testdata$x2), 1.1 * max(testdata$x2)))
   
    y_trained <-predict(result, cbind(testdata$x1, testdata$x2), type="vector")
    cat("precision", mean(datamat$y == y_trained) * 100, "\n")
}

plotData <- function(X, y)
{
    # ----- plot -----
    plot(0, 0, type = "n", xlim = c(0.9 * min(X[, 1]), 1.1 * max(X[, 1])), ylim = c(0.9 * min(X[, 2]), 1.1 * max(X[, 2])),xlab = "x1", ylab = "x2")   
    # points(X, col = ifelse(y == 1, "black", "yellow"), pch = ifelse(y == 1, 3, 16))
    points(X, col = c("black", "orange")[y+1], pch = c(3, 16)[y+1])
    legend("topleft", legend = c("y = 1", "y = 0"), pch = c(3, 16), col = c("black", "orange"))
}

結果
デフォルト設定でSVM
 > ex6()

Call:
svm(formula = y ~ ., data = testdata)


Parameters:
   SVM-Type:  C-classification
 SVM-Kernel:  radial
       cost:  1
      gamma:  0.5

Number of Support Vectors:  406

precision 90.38239

次にcost = 10, gamma = 1で
> ex6()

Call:
svm(formula = y ~ ., data = testdata, cost = 10, gamma = 1)


Parameters:
   SVM-Type:  C-classification
 SVM-Kernel:  radial
       cost:  10
      gamma:  1

Number of Support Vectors:  151

precision 98.84125 

[Rの]Andrew先生の機械学習4[練習]

3回目後半と4回目をまとめてしまう。
3回目後半はニューラルネットワークのフィードフォワード。
4回目はニューラルネットワークの学習。

対象は前回と同様手書き数字

今回作るニューラスネットワークの構成は以下
・手書き数字の画像20*20ピクセルなので、入力素子400個
・隠れ層が1層あり、隠れ素子素子25個
・1~10の数字を識別するので、出力素子は10個

・学習させるニューラルネットワークのパラメータは、バイアス項を除くと400*25の行列と25*10の行列。
・非線形出力関数はシグモイド関数




R
学習するニューラスネットワークのパラメータはランダムに初期化しておく。 
全てが0だったりすると学習が全く進まなかったり、ローカルミニマムが大量にあるため、こうするという記憶。
randInitializeWeights <- function(L_in, L_out)
{
  epsilon_init = 0.12
　　return(matrix(runif(L_out * (1 + L_in), min = -epsilon_init, max = epsilon_init), L_out, 1 + L_in))
}

コスト関数。誤差逆伝播法で計算する勾配関数も。
nnCostFunction <- function(nn_params, input_layer_size, hidden_layer_size, num_labels, X, y, lambda)
{
  m <- length(X[, 1])
  Theta1 <- matrix(nn_params[1:(hidden_layer_size * (input_layer_size + 1))], hidden_layer_size, input_layer_size + 1)
  Theta2 <- matrix(nn_params[(length(Theta1) + 1):length(nn_params)], num_labels, hidden_layer_size + 1)

  a1 <- X
  a1 <- cbind(matrix(1, length(a1[, 1]), 1), a1)
  z2 <- a1 %*% t(Theta1)
  a2 <- sigmoid(z2)
  a2 <- cbind(matrix(1, length(a2[, 1]), 1), a2)
  z3 <- a2 %*% t(Theta2)
  a3 <- sigmoid(z3)

  yvector <- matrix(0, length(y[, 1]), num_labels)
  for (i in 1 : num_labels)
  {
    yvector[, i] <- y[, 1] == i
  }
 
  # ----- cost function -----
  J <- sum(sum(-yvector * log(a3) - (1 - yvector) * log(1 - a3))) / m
 
  # ----- regularized cost function -----
  J = J + lambda / (2 * m) * (sum(sum(Theta1[, 2:length(Theta1[1, ])] ^ 2))+ sum(sum(Theta2[, 2:length(Theta2[1, ])] ^ 2)))

  # ----- back propagation -----
  delta3 = (a3 - yvector)
  delta2 = delta3 %*% Theta2 * sigmoidGradient(cbind(matrix(1, length(z2[, 1]), 1), z2))
  Delta1 = t(delta2[, 2:length(delta2[1, ])]) %*% a1
  Delta2 = t(delta3) %*% a2

  Theta1_grad = Delta1 / m
  Theta2_grad = Delta2 / m;

  # ----- regularized derivatives -----
  Theta1[, 1] <- 0
  Theta2[, 1] <- 0
  Theta1_grad <- Theta1_grad + lambda * Theta1 / m
  Theta2_grad <- Theta2_grad + lambda * Theta2 / m
  grad = cbind(matrix(Theta1_grad, 1, length(Theta1_grad)), matrix(Theta2_grad, 1, length(Theta2_grad)))

  return(cbind(J, grad))
}

nnCostFunction.Value <- function(nn_params, input_layer_size, hidden_layer_size, num_labels, X, y, lambda)
{
  return(nnCostFunction(nn_params, input_layer_size, hidden_layer_size, num_labels, X, y, lambda)[1])
}

nnCostFunction.Gradient <- function(nn_params, input_layer_size, hidden_layer_size, num_labels, X, y, lambda)
{
  result <- nnCostFunction(nn_params, input_layer_size, hidden_layer_size, num_labels, X, y, lambda)
  return(result[2:length(result)])
}

ニューラルネットワークのフィードフォワード
nn_predict <- function(Theta1, Theta2, X)
{
    m <- length(X[, 1])
    num_labels <- length(Theta2[, 1])
   
    a1 <- cbind(matrix(1, m, 1), X)
    a2 <- sigmoid(a1 %*% t(Theta1))
    a2 <- cbind(matrix(1, length(a2[, 1]), 1), a2)
    a3 <- sigmoid(a2 %*% t(Theta2))

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

メイン関数
ex4 <- function()
{
     require(R.matlab)
    datamat <- readMat('ex4data1.mat')
    X <- datamat{dollar}X
    y <- datamat{dollar}y
    displayData(X)

    input_layer_size <- 400  # 20x20 Input Images of Digits
    hidden_layer_size <- 25 # 25 hidden units
    num_labels <- 10 # 10 labels, from 1 to 10  

    # ----- randomly initialize nn parameters ----- 
    initial_Theta1 <- randInitializeWeights(input_layer_size, hidden_layer_size);
    initial_Theta2 <- randInitializeWeights(hidden_layer_size, num_labels);
    initial_nn_params <- cbind(matrix(initial_Theta1, 1, length(initial_Theta1)), matrix(initial_Theta2, 1, length(initial_Theta2)))

    # ----- train Neural Network -----
    lambda <- 1
    result <- optim(par = initial_nn_params, fn = nnCostFunction.Value, gr = nnCostFunction.Gradient, input_layer_size = input_layer_size, hidden_layer_size = hidden_layer_size, num_labels = num_labels, X = X, y = y, lambda = lambda, method= "CG")

    # ----- extract trained parameters -----
    nn_params <- result$par
    Theta1 <- matrix(nn_params[1:(hidden_layer_size * (input_layer_size + 1))], hidden_layer_size, input_layer_size + 1)
    Theta2 <- matrix(nn_params[(length(Theta1) + 1):length(nn_params)], num_labels, hidden_layer_size + 1)
   
    # ----- feed forward with NN -----
    pred <- nn_predict(Theta1, Theta2, X)
    cat("Training Set Accuracy : ", mean(pred == y) * 100, "\n")
    return(summary(result))
}



結果
初期パラメータがランダムなので結果は毎回変わる。
> ex4()
Training Set Accuracy :  95.14
            Length Class  Mode  
par         10285  -none- numeric
value           1  -none- numeric
counts          2  -none- numeric
convergence     1  -none- numeric
message         0  -none- NULL  
>

学習にそれなりに時間がかかってしまうが、、誤差逆伝播法を使わなかった場合さらに時間がかかる。
これを見るには、上のメインのoptim関数の引数のうちgr = nnCostFunction.Gradientを消すとよい。そうすると、optim関数がパラメータひとつひとつに対して数値微分を行うため、計算量が恐ろしく増大してしまう。

[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で書いたものより識別率が上がった。
講義を受けているときには、このあたりで機械学習している実感が出てきてテンションあがった。