Rで実装する多層パーセプトロン
あっ, どーも僕です.
修士論文書きました. 公聴会まで少し時間があるのでブログ書きます.
本題に入る前に
Rのdplyrパッケージが爆速なのですね.
デフォルトパッケージ推しのわたしもさすがに乗り換えました.
ところで, dplyrで外生的に分解軸を設定するにはどうしたらよいのでしょう？
formulaが使えるものはformulaで, plyrならas.quotedできます.
誰かご存知でしたら教えてください.
多層パーセプトロン
はじめに
「Rは最遅」だとあちらこちらで聞きますが, パーセプトロンぐらいなら十分にRで楽しめるようです. Rcppは勉強不足です.
sylvan5 (sylvan5) on TwitterさんがPythonでパーセプトロンを実装しています.
Pythonは読めない, かつ, パーセプトロンに挫折したことあるわたしでもわかりやすいと感じるほど丁寧に記述されています.
多層パーセプトロンで手書き数字認識 - 人工知能に関する断創録
ブックマーク先では, 活性化関数の一階微分を使用することで学習の効率がよくなる, がなぜかわからないとしていました.
追記をみるとどうやらPRMLとは目的関数が違っていたためだそうです.
スクリプト
多層パーセプトロンのクラス
# 3層パーセプトロン # クラスに関しては逆引きハンドブックを参照 # 記号とか間違っていたらごめんな MultiPerceptronBP <- setRefClass ( Class = "mpbp", #付きははじめに設定して... fields = list (trainData = "matrix", # trainSignal = "matrix", # NUM_INPUT = "integer", NUM_HIDDEN = "integer", # NUM_OUTPUT = "integer", WeightKJ = "matrix", WeightJI = "matrix", isInitialize = "logical", act1 = "character", # 出力層の活性化関数 Act1 = "function", Act1.deriv = "function", act2 = "character", # 中間層の活性化関数 Act2 = "function", Act2.deriv = "function", sse = "numeric"), methods = list ( Initialize = function () { if (!all (trainData[,1,drop = TRUE] == 1)) { trainData <<- cbind (rep (1, nrow (trainData)), trainData) } # trainData <<- sweep (trainData, 2, colMaxs (trainData), "/") 正規化はどうするよ... NUM_INPUT <<- ncol (trainData) NUM_OUTPUT <<- ncol (trainSignal) WeightKJ <<- matrix (runif (NUM_OUTPUT * NUM_HIDDEN, -1, 1), NUM_OUTPUT, NUM_HIDDEN) WeightJI <<- matrix (runif (NUM_HIDDEN * NUM_INPUT, -1, 1), NUM_HIDDEN, NUM_INPUT) setAct1 () setAct2 () sse <<- SSE () isInitialize <<- TRUE }, setAct1 = function () { if (act1 == "tanh") { Act1 <<- tanh Act1.deriv <<- function (x) 1 - Act1 (x) ** 2 } else if (act1 == "sigmoid") { Act1 <<- function (x) 1 / (1 + exp (-x)) Act1.deriv <<- function (x) Act1 (x) * (1 - Act1 (x)) } else { Act1 <<- function (x) x Act1.deriv <<- function (x) { if (is.matrix (x)) { d <- attr (x, "dim") return (matrix (1, d[1], d[2])) } else { rep (1, length (x)) } } } }, setAct2 = function () { if (act2 == "tanh") { Act2 <<- tanh Act2.deriv <<- function (x) 1 - Act1 (x) ** 2 } else if (act1 == "sigmoid") { Act2 <<- function (x) 1 / (1 + exp (-x)) Act2.deriv <<- function (x) Act2 (x) * (1 - Act2 (x)) } else { Act2 <<- function (x) x Act2.deriv <<- function (x) { if (is.matrix (x)) { d <- attr (x, "dim") return (matrix (1, d[1], d[2])) } else { rep (1, length (x)) } } } }, SSE = function () { Z <- Act2 (tcrossprod (WeightJI, trainData)) Z[1, ] <- 1 Y <- Act1 (WeightKJ %*% Z) sse <<- c (crossprod (c ((trainSignal - t (Y))))) * .5 return (sse) }, ForwardPropagation = function (x.biased) { z <- Act2 (WeightJI %*% x.biased) z[1] <- 1 y <- Act1 (WeightKJ %*% z) return (list (y = y, z = z)) }, BackPropagation = function (y, z, signal) { delta.output <- c (y - signal) * c (Act1.deriv (WeightKJ %*% z)) delta.hidden <- c (Act2.deriv (z)) * c (crossprod (delta.output , WeightKJ)) return (list (o = delta.output, h = c (delta.hidden))) }, fit = function (eta = .02, epochs = 1) { cat ("更新前二乗誤差 = ", sse, "\n") for (loop in seq_len (epochs)) { for (m in 1:nrow (trainData)) { # 入力データ x.biased <- trainData[m, , drop = TRUE] signal <- trainSignal[m, , drop = TRUE] # 順伝播 fp <- ForwardPropagation (x.biased = x.biased) z <- fp$z y <- fp$y # 逆伝播 bp <- BackPropagation (z = z, y = y, signal = signal) delta.output <- bp$o delta.hidden <- bp$h # 重み更新量の計算 WeightJI.delta <- - eta * tcrossprod (delta.hidden, x.biased) WeightKJ.delta <- - eta * tcrossprod (delta.output, z) # 重みの更新 WeightJI <<- WeightJI + WeightJI.delta WeightKJ <<- WeightKJ + WeightKJ.delta } } # 二乗誤差 sse <<- SSE () cat ("更新後二乗誤差 = ", sse, "\n") }, predict = function (X) { X <- cbind (rep (1, nrow (X)), X) Z <- Act2 (tcrossprod (WeightJI, X)) # [J, N] Z[1, ] <- 1 Y <- Act1 (WeightKJ %*% Z) # [K, N] return (list (y = t (Y), z = t (Z))) }) )
グラフ
# ************************************************************** # example : 関数近似, sin # ************************************************************** # 訓練データ X <- matrix (c (seq (-pi, pi, length.out = 30)), 30, 1) Y <- sin (X) # fit perSIN <- MultiPerceptronBP$new (trainData = X, trainSignal = Y, act1 = "linear", act2 = "tanh", NUM_HIDDEN = 5L) perSIN$Initialize () perSIN$fit (eta = .01, epochs = 5000) #テストデータ test <- matrix (c (seq (-pi, pi, length.out = 100)), 100) y <- perSIN$predict (test)$y #グラフ化 plot (cbind (X, Y), pch = 19, col = "dodgerblue1", cex = 1.5, xlab = "x", ylab = "y", tcl = .5, cex.lab = 1.5) points (test, y, type = "l", lwd = 2, col = "deeppink") legend ("topleft", bty = "n", legend = c ("訓練データ", "予測値"), pch = c(19, -1), lty = c(-1, 1), col = c ("dodgerblue1", "deeppink"), cex = 2, lwd = 3, pt.cex = 2)
おまけ
# ************************************************************** # example : パターン認識, iris # ************************************************************** myiris <- iris[(1:75) * 2, ] X <- as.matrix (myiris[,-5]) Y <- as.matrix (model.matrix (~ 0 + Species, data = myiris)) perIRIS <- MultiPerceptronBP$new (trainData = X, trainSignal = Y, act1 = "sigmoid", act2 = "tanh", NUM_HIDDEN = 10L) perIRIS$Initialize () perIRIS$fit (eta = .01, epochs = 3000) XX <- as.matrix (iris[-(1:75) * 2, -5]) table (c (iris[-(1:75) * 2, 5]), max.col (perIRIS$predict (X = XX)$y)) 1 2 3 1 25 0 0 2 1 22 2 3 0 0 25
おわりに
最近どーも, 自分で新しい勉強ネタを見つけずに人に乗っかてしまいます.
2014年はその辺を変えていけたらなと思います.
