読者です 読者をやめる 読者になる 読者になる

Rで実装する多層パーセプトロン

R

あっ, どーも僕です.

修士論文書きました. 公聴会まで少し時間があるのでブログ書きます.

本題に入る前に

Rのdplyrパッケージが爆速なのですね.
デフォルトパッケージ推しのわたしもさすがに乗り換えました.

ところで, dplyrで外生的に分解軸を設定するにはどうしたらよいのでしょう?
formulaが使えるものはformulaで, plyrならas.quotedできます.
誰かご存知でしたら教えてください.


多層パーセプトロン

はじめに

「Rは最遅」だとあちらこちらで聞きますが, パーセプトロンぐらいなら十分にRで楽しめるようです. Rcppは勉強不足です.

sylvan5 (sylvan5) on TwitterさんがPythonパーセプトロンを実装しています.
Pythonは読めない, かつ, パーセプトロンに挫折したことあるわたしでもわかりやすいと感じるほど丁寧に記述されています.

多層パーセプトロンで手書き数字認識 - 人工知能に関する断創録


ブックマーク先では, 活性化関数の一階微分を使用することで学習の効率がよくなる, がなぜかわからないとしていました.
追記をみるとどうやらPRMLとは目的関数が違っていたためだそうです.

二乗誤差を微分

PRMLは持っていないのですが, 二乗誤差の場合ならsylvan5さんのブログに掲載されていました.
そこで, 二乗誤差の場合に本当に出力層の活性化関数の一次微分が出てくるのかを確認してみました.

自分用に作ったメモを載せておきます. 一応それらしい式が導出されますが, 自信ないので間違っていたら教えてください.

Rで実装

Rで多層パーセプトロンを実装しました. sylvan5さんのPythonのプログラムを参考にしました.

はじめてリファレンスクラスを使ったのですがなんとか動いてよかったです. 継承とかを理解していればもっとわかりやすくつくれたのかな?

次の図はsin関数の近似結果です.
ちゃんと動作してくれそうなことがわかります.


f:id:aaaazzzz036:20140202161850p:plain

実装の注意点を挙げるとすると3つあります.

  • 恒等関数はちゃんと同じ次元の配列で返すようにしなければいけない
  • integerは数字にL付きでわたす
  • crossprodが返すのはmatrix

スクリプト

多層パーセプトロンのクラス
# 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年はその辺を変えていけたらなと思います.