分からんこと多すぎ

何故か就職できたので技術磨く。一人前の多重線形代数使いを目指しつつ、機械学習、データマイニングをやる

ゼロから始めるDeepLearning_その2_ホップフィールドネットワークとは

MachineLearning R

はじめに

学部四年生向け。

ゼロから始めるDeepLearning_その1_ニューラルネットとは - 分からんこと多すぎ
→（Auto Encoder）
→（Denoising AutoEncoder）
→ホップフィールドネットワーク（この記事）
→ゼロから始めるDeepLearning_その3_ボルツマンマシンとは - 分からんこと多すぎ
→Restricted ボルツマンマシン（後日）
→（Gaussian Binary - Restricted Boltzmann Machines）
→（Deep Belief Networks）
→（Deep Neural Networks）
→畳み込みニューラルネット（後日）

太線以外は読み飛ばしてOK

知識に乏しく深く踏み込めなかったので、間違ってたらごめんなさい。

参考図書

An Introduction to Restricted Boltzmann Machines

CiNii 論文 -  統計的モデルとしてのボルツマンマシン

ニューラルネットワーク情報処理―コネクショニズム入門、あるいは柔らかな記号に向けて

ニューラルネットワーク情報処理―コネクショニズム入門、あるいは柔らかな記号に向けて

やっぱり麻生先生はすごい。長岡先生も説明が分かりやすい。

この記事と次の記事の概要

ボルツマンマシンには、ホップフィールドネットワークという前身がある。

ホップフィールドネットワークとは、1つの行列に複数の画像を埋め込む技術である。

ホップフィールドネットワークは、ノイズありの画像から元々のノイズなしの画像を復元できる。

ホップフィールドネットワークには、古典的な評価関数が存在し、これはエネルギー関数と呼ばれている。

エネルギー関数の微分をすると、次の時刻での素子の出力（要するにノイズなしの画素値）を求めることができる。

この更新は、マルコフ連鎖モンテカルロ法（MCMC）である。

個人的には面白かったが、この記事を読まなくてもRBMが理解できる説はある。
（エネルギー関数だけは重要）

今回できるようになること

まず画像データをもらってくる。標準画像

今回はこの画像データを、脳的なもの（ネットワーク：マルコフ確率場）に記憶させるアルゴリズムの話をする。

ここでは、2つの画像データを1つの脳的なものに記憶させる例を上げる。

記憶させたイメージその１（LENAの一部111*111画素）

記憶させたイメージその２（Cameramanの一部111*111画素）

この脳は連想記憶というものをもっていて、
適当なデータを投げ込むと、記憶の中にあるものを勝手に思い出す。

脳に入ったデータ（いわば視覚）

この情報（脳に入ったデータ）が、脳の中で記憶と照らし合わされることを、想起と呼ぶことにする。
（何か覚えがあるけど、あれってなんだっけ…と頭を捻っている感じ。
　想起するごとにイメージが鮮明になり、最終的に思い出す（元画像が復元される））
（注・勝手に想起とか言っているだけで、全然厳密な意味ではない）

5,000回想起後

10,000回想起後　　　　　　　

100,000回想起後

適当な視覚データから、記憶させた画像が取り出せたことが分かる。
（枯れ木が幽霊に見える原理はたぶんこんな感じ）

こういうことができると、ノイズの除去なんかに役に立つ。
（砂嵐が混ざった画像から、綺麗な画像を取り出せる）

今回のメインは、どうやってネットワークに記憶を埋め込むのかという話。

ちなみに、画像は本当に実験した結果。

同じ脳（ネットワーク）は、もう1つの画像（カメラマン）も同時に覚えているので、
こういう画像を見せると

こういう画像を思い出してくれる。

概要・ボルツマンマシンってなに？

1985 年にヒントンが作ったアルゴリズム。

ボルツマンマシンは、グラフ構造上での状態遷移を確率的に扱うためのものである。

ホップフィールドネットワークというものを下敷きにしている。

意味がわからないので、ホップフィールドネットワークの説明する。

ホップフィールドネットワーク

上図のような、交差点に相当するもの（ノード：今回はｉやｊ）を、
道路に相当するもの（エッジ：今回はwij）がつないだものを、グラフ構造と呼ぶ。

今回はw_{ij}=w_{ji}なので、無向グラフと呼ばれる。（対義語は有向グラフw_{ij} \neq w_{ji}）

グラフ理論とかは、このへんでちょっと触れた。
グラフ（ネットワーク状）構造のクラスタリング（グループ分け）基礎：The Markov Cluster Algorithm（MCL）まとめた - 分からんこと多すぎ


ノードｉには観測可能な出力x_{i}={-1,1}と、観測不可能な内部状態yiが存在する。

結合係数wは、ノードｉの観測可能な出力x_{i}={-1,1}と
記憶させるイメージ名ｓを用いて、次の式で設定する。

f:id:rishida:20140301030330p:plain

安定なネットワークを作るために、自己想起w_{ii}=0


ホップフィールドネットワークとは、以上のような想定の下で、
ネットワークに情報を記憶させるもののことである。


ホップフィールドネットワークをフローチャート的に書くと、以下のようになる。


0.ネットワークに記憶させるデータにしたがって、結合係数w_{ij}を設定する

1.ランダムにノードｉを選ぶ

2.ノードの状態y_{i} = \sum_{j=1}^{N} w_{ij} x_{j} - \theta_{i} を計算する

3.ノードの状態y_{i} < 0なら、出力x_{i} = -1
　ノードの状態y_{i} = 0なら、出力x_{i}はそのまま
　ノードの状態y_{i} > 0なら、出力x_{i} = 1

4.以上を繰り返すと、最終的に記憶させた状態（ある定常状態）に至る
　（最初に記憶させたデータが取り出せる＝人間の記憶の仕組みっぽい）
　（一種のセルオートマトンとも言える）


図的に説明する。

まず、25個のノードに、以下の2つのイメージを記憶させる。

この時、重みｗは以下のように設計される（25*25の行列）
このｗは計算式から一意に定まる。

この重みwが定められている時、
ランダムな初期状態（下図の一番上の状態）からフローチャートの2.と3.の更新を繰り返すと、一番下の状態に至る。

一番下の状態は、最初に記憶させたイメージの片方と一致する。
（完全に一致させるのは、何か色々難しいらしい）

ネットワークに乱数を突っ込むと、予め定めた定常状態に至る ＝ ネットワークが情報を記憶している。

このように、相互に結合のあるネットワークに情報を記憶させるのが、
ホップフィールドネットワークである。
（Rのコードは最後に載せる）


ちなみに、最終的に定常状態に至るのは、
更新によってネットワークのエネルギー関数が常に減少する（安定な状態になる）ためである。

　　

エネルギー関数

　　エネルギー関数とは、要するに評価関数のことで、以下の様なものである。
　　（エネルギー関数はボルツマンマシンでも出てくるので重要）

　　E = - \frac{1}{2} \sum_{ij} w_{ij} x_{i} x_{j} + \sum_{i} \theta_{i} x_{i}

　　この関数が小さな値を取る時、ネットワークが安定となることから、この関数はエネルギー関数と呼ばれる。

　　ホップフィールドネットワークでは、この関数の鞍点のうちの１つを求めることができるが、
　　必ずしも最下点であるとは言えない。


　　以下、エネルギー関数が更新によって常に減少することの説明

　　　　2.ノードの状態y_{i} = \sum_{j=1}^{N} w_{ij} x_{j} - \theta_{i} を計算する

　　　　3.ノードの状態y_{i} < 0なら、出力x_{i} = -1
　　　　　ノードの状態y_{i} = 0なら、出力x_{i}はそのまま
　　　　　ノードの状態y_{i} > 0なら、出力x_{i} = 1

　　という更新がホップフィールドネットワークの中では行われる。

　　実は、2番の式はエネルギー関数のx_{i}による微分と等しい。

　　y_{i} = \sum_{j=1}^{N} w_{ij} x_{j} - \theta_{i} = - \frac{\partial E}{\partial x_{i}}

　　1個の変数x_{1}を持つある関数f (x_{1})のテイラー展開は、
　　関数f (x_{1})の微分を用いて、以下のように表される。

　　f (x_{1} + \Delta x_{1} ) = f (x_{1} ) + \frac{\partial \hat{f}}{\partial x_{1}} ( x_{1} ) \Delta x_{1} + \cdots

　　したがって、エネルギー関数は以下のように一次近似できる。
　　（エネルギー関数の二次微分以降は0）

　　E (x_{i} + \Delta x_{i} ) = E (x_{i} ) - y_{i} \Delta x_{i}

　　今、y_{i} = \sum_{j=1}^{N} w_{ij} x_{j} - \theta_{i} = - \frac{\partial E}{\partial x_{i}} < 0だったとする。

　　y_{i} < 0 なので x_{i} は -1 の値に変更される。

　　①.x_{i} = 1からy_{i} < 0でx_{i} = -1に変更された場合 \Delta x_{i} = -2
　　　したがって、E (x_{i} + \Delta x_{i} ) = E (x_{i} ) - y_{i} \Delta x_{i} = E (x_{i} ) - ( - ) (-2) < E (x_{i} )

　　②.x_{i} = -1からy_{i} < 0でx_{i} = -1に変更された場合 \Delta x_{i} = 0
　　　したがって、E (x_{i} + \Delta x_{i} ) = E (x_{i} ) - y_{i} \Delta x_{i} = E (x_{i} ) - ( - ) (0) = E (x_{i} )

　　①,②より、E (x_{i} + \Delta x_{i} ) = E (x_{i} ) - y_{i} \Delta x_{i} \leq E (x_{i} )

　　y_{i} > 0の場合も同様である。

　　以上のことより、出力ｘの更新によって（想起を重ねることによって）
　　エネルギー関数が常に減少することが示された。



思いがけず重たくなったので、今回はここで終了。

次回はボルツマンマシン。

実装

ホップフィールドネットワークをRで実装する。

#使い方
if(FALSE){
  library('bmp')
  data1 <- read.bmp('Downloads/mono/LENNA.bmp')
  data2 <- read.bmp('Downloads/mono/Cameraman.bmp')
  #実際には大きすぎてメモリが足りないので
  #dataの一部を小さく切り出してください
  Mem <- list(data1, data2)
  huga <- hop.net(memory = Mem)
  image(huga$state[[100]])
}

#Hopfield networks
  #必須項目
  #non
    #オプション
    #dims.vec:記憶させるイメージの大きさ
    #p:記憶させるイメージの1の割合
    #theta:発火のしきい値の値
    #mnum:記憶させるものの数(ニューロン数の15%あたりが限界らしい)
    #memory:記憶させるイメージのリスト（0と1だけで表現される同じサイズの行列のリスト）
hop.net <- function(dims.vec = c(5,5), p = 0.3, theta = 0.5, mnum = 2, memory = NULL){
  #set memory
  if(is.null(memory)){
    memory <- list(NULL)
    #layout(matrix(c(1:mnum),nrow=mnum))
    for(i in 1:mnum){
      memory[[i]] <- matrix(rbinom(n = prod(dims.vec), prob = p, size = 1), nrow = dims.vec[1], ncol = dims.vec[2])
      memory[[i]][memory[[i]] < 1] <- -1
      #image_text(memory[[i]])
    }
  }else{
    mnum <- length(memory)
    for(i in 1:mnum){
      memory[[i]] <- round(memory[[i]])
      memory[[i]][memory[[i]] < 1] <- -1
    }
    dims.vec <- dim(memory[[1]])
  }

  #set weight
  w <- matrix(0, nrow = prod(dims.vec), ncol = prod(dims.vec))
  for(i in 1:mnum){
    w <- w + (as.vector(memory[[i]]) ) %*% t(as.vector(memory[[i]]))
  }
  diag(w) <- 0
  #image_text(w)
  browser()

  #remember memory
  state <- list(NULL)
  #state[[1]] <- matrix(rbinom(n = prod(dims.vec), prob = p, size = 1), nrow = dims.vec[1], ncol = dims.vec[2])
  state[[1]] <- memory[[1]] + rnorm(mean=0,sd=0.3,n=prod(dims.vec))
  state[[1]][state[[1]] > 1] <- 1
  state[[1]][state[[1]] < 1] <- -1
  huga <- list(NULL)
  for(i in 1:100000){
    print(paste(i, 'th iteration'))
    #ランダムにノードを選択
    id <- sample(x = c(1:prod(dims.vec)), size = 1)
    state[[i+1]] <- state[[i]]
    #ノードの出力を更新
    hoge <- w[id,] %*% as.vector(state[[i]]) - theta
    hoge <- sign(hoge) * sign(ceiling(abs(hoge)))
    #hoge <- ceiling(sign(hoge))
    state[[i+1]][id] <- hoge
    for(j in 1:mnum){
      #記憶が連想できたら終了
      if( sum(memory[[j]] == state[[i+1]]) == prod(dims.vec)){
        #dev.new()
        print('image is remembered')
        #layout(matrix(c(1:mnum),ncol=1))
        #for(k in 1:mnum){
          #image_text(memory[[k]])
        #}
        #dev.new()
        #layout(matrix(c(1:4),ncol=1))
        #for(k in c(1,(i-1),(i),(i+1))){
          #image_text(state[[k]])
        #}
        ans <- list(memory,state)
        names(ans) <- c('memory','state')
        return(ans)
      }
    }
    if( (i %% 1000) == 1){
      huga[[((i %/% 1000) + 1)]] <- state[[i]]
    }
    state[1:i] <- 1
  }

  #返り値
  #ans <- list(memory,state)
  ans <- list(memory,huga)
  names(ans) <- c('memory','state')
  return(ans)
}

#描画関数
#ネット上で公開されていたはずだが元コードが見つからない
image_text <- function(x,list.I=0,list.J=0,TEXT=TRUE,...){
  at_v <- 1/(2*(ncol(x)-1))
  at_v <- seq(0,1+2*at_v,by=at_v)-at_v
  at_h <- 1/(2*(nrow(x)-1))
  at_h <- seq(0,1+2*at_h,by=at_h)-at_h
  xy <- do.call("rbind",lapply(at_v[seq(2,length(at_v),by=2)],function(x,y) cbind(x,y),y=rev(at_h[seq(2,length(at_h),by=2)])))
  image(1-t(x[nrow(x):1,ncol(x):1])[ncol(x):1,])
  if(TEXT==TRUE){
        text(round(xy,2),labels=x)
  }
  abline(v=at_v[seq(3,length(at_v)-1,by=2)],h=at_h[seq(3,length(at_h)-1,by=2)],col=8)
  if(length(list.I)!=1){
        tmp<-unlist(lapply(length(list.I):1,function(i){length(list.I[[i]])}))
        abline(h=at_h[c(1,2*cumsum(tmp)+1,2*sum(tmp)+1)],col=1,lwd=3,lty=2)
  }
  if(length(list.J)!=1){
        tmp<-unlist(lapply(1:length(list.J),function(i){length(list.J[[i]])}))
        abline(v=at_v[c(1,2*cumsum(tmp)+1,2*sum(tmp)+1)],col=1,lwd=3,lty=2)
  }
  #return(list(v=at_v,h=at_h))
  return(1)
}

疑問

ホップフィールドネットワークのエネルギー関数は、なぜあの形になったのか。

w_{ij}とx_{i},x_{j}の相関を一致させる一番簡単な式だから？

ゼロから始めるDeepLearning_その1_ニューラルネットとは

対象とする人

ディープラーニングすごい！　←聞き飽きた

チュートリアルあるよ！　←ふわっとしすぎて具体的なところが分からん

こういう論文あるよ！　←読めるわけないだろ

そういう人向け。（たぶん学部四年程度向け）

ニューラルネット初学者が、書ききるまで怪しいところ満載でも突っ走ります。

ニューラルネット（この記事）
→（AutoEncoder）
→（DenoisingAutoEncoder）
→ホップフィールドネットワーク
→ボルツマンマシン
→Restrictedボルツマンマシン
→（Gaussian Binary - Restricted Boltzmann Machines）
→（DeepBeliefNetwork）
→（DeepNeuralNetworks）
→畳み込みニューラルネット（後日）

までやる。

太線以外は読み飛ばしてOK

本文中では怖い式は使わない。（Appendixに書く）

分からない人がいたら、説明を追加する。

正直算数が覚束ないので、間違ってたらごめんなさい。

ディープラーニングとは

①.ディープラーニングとは、”ニューラルネット（っぽいもの）”を”多層に”積み上げたものによる学習器の総称である。

②.特徴量を自動で作ってくれる。

③.画像認識とかの記録を驚くほど塗り替えた、新進気鋭の技術。

今回は、まずニューラルネットワークの話をする。

ニューラルネットをとりあえず使ってみる

ニューラルネットを使うと、任意の関数が近似できる。

クラス判別の場合は、ある関数が境界面をつくっていると考える。

境界面をつくっているある関数を推定する問題と見れば、任意の関数を近似する問題と同様に解ける。

関数近似

図的にはこんな感じ。

横軸が入力 x 、縦軸が出力 y .

図中の黒線は推定したい関数y=sin(x)

図中の赤点は推定された関数（初期のニューラルネットによる関数）

状況：赤点の横軸上の値 x と黒線上の値 y=sin(x) + noise が与えられている

目的：赤点の縦軸上の値（出力値 y' ）が黒線に近づくように未知の関数（今回はsin関数）を求めたい

ニューラルネットの役割：未知の関数そのものを、近似する関数になること

更新を重ねると、以下のように推定した関数が変化していく。

ニューラルネットとは

ニューラルネットとは、多層パーセプトロンのこと。

そもそもパーセプトロンってなに？

　　1962年にローゼンブラットが作ったアルゴリズム。

　　y \left( x \right) = f \left( w^{T} \phi \left( x \right) \right)

　　で表される線形識別モデルのこと。

　　ここで非線形活性化関数fは

　　 f \left( a \right) = \begin{cases} +1, \; a \geq 0 \\ -1, \; a < 0 \end{cases}

　　というステップ関数である。

　　つまり、

　　入力データxを、

　　関数Φによって特徴量に変換して、

　　重みwで足しあわせ、

　　ステップ関数で二値に分類するものである。

尤度関数がパラメータの凸関数にならないため、勾配法を何回も解いて解に近づけることになる。

当然過学習もする。（EarlyStoppingという、途中で学習を打ち切る方法で対処したりする）

良いところは、あらゆる関数を近似できることと、
学習し終わったモデル自体は軽量なので、新しいデータを高速に処理できること。

ニューラルネットを図的に書くと、こんな感じになる。

入力変数ベクトルxが与えられた時、線形変換wをいじることで、出力変数ベクトルyを目標ベクトルtに近づけるのが目的となる。

したがって、以下の誤差関数を最小化すれば良い。

一個ずつ説明していく。

ニューラルネットの仕組み

①.まず、d次元のデータベクトルがやってくる

あと、関数の切片とかをいじるために、バイアスパラメータx0=1をくっつける。

②.入力ベクトルが線形変換されて、各活性aに入る

つまり，

a_{j} = \sum_{i=0}^{d} w_{ji}^{(1)} x_{i}

③.各活性がしきい値を超えたら発火する（1の値を取る）

つまり，

z_{j} = h( a_{j} )

非線形関数h(・)はロジスティックシグモイド関数や、tanh関数のようなシグモイド関数を用いる。

例えば，標準シグモイド関数は，こんな関数．

\sigma \left( a \right) = \frac{1}{1 + \exp \left( - a \right)}

非線形変換がないと任意の関数を近似することができないので、何らかの関数が必要とされる。

特にこれらのシグモイド関数が用いられるのは、微分した時の値の計算が簡単なためである。

指数関数族は微分にて最強。

④.隠れ層zが線形変換されて、出力ユニット活性aに入る

ついでに隠れ層にバイアスz0=1をくっつけて、

a_{k} = \sum_{j=0}^{m} w_{kj}^{(2)} z_{j}

⑤.（クラス判別問題の場合）各活性がしきい値を超えたら発火する

y_{k} \left( x,w \right) = \sigma \left( a_{k} \right)

（関数近似の場合は、ここでシグモイド関数を用いない）

したがって、①から⑤までの式をくっつけると、

図からも見て取れるように、ニューラルネットワークは、パーセプトロンが二段重なった構造になっている。

パーセプトロンとの違いは、ステップ関数ではなく、連続で非線形なシグモイド関数を用いる点である。

連続関数を用いているため、パラメータに関して微分可能であり、高速な学習を可能としている。

また、シグモイド関数の非線形領域を用いた場合のみ、ニューラルネットは万能の関数近似器になる。

（シグモイド関数の真ん中ら辺は、線形領域なので、ただの線形変換しかしない。
　ただの線形変換によってニューラルネットを作った時、
　モデルの表現力が著しく損なわれることを、ローゼンブラットが証明した）

隠れ層の素子数を、入出力次元数よりも小さくすると、
主成分分析のように特徴量が圧縮される（AutoEncoder）。

逆に、隠れ層の素子数を、入出力次元数よりも大きくすると、
スパースコーディングのように特徴量がスパースになる（SparseAutoEncoder）。

どうやって学習するの？

誤差関数は、次の式で与えられる。

これをwを変化させることで最小化する。

関数の最小値は、微分して勾配が0になる点にあるので、

\nabla E(w) = 0

最急降下法を用いると、wの更新式は以下のようになる。

w_{(n+1)} = w_{(n)} - \eta \nabla E(w_{(n)})

これはすべてのデータを一度に使うバッチ学習というもので、
ニューラルネットの学習にはふさわしくない。
（計算オーダーが、後に紹介する最適化よりもはるかに大きい）

そこで今回は、得られたデータごとに逐次、勾配降下法を行うこととする。
（オンライン学習）

n番目のデータが来るごとに、誤差E_{n}を計算する。

その和が最終的な誤差関数となる。

E ( w) = \sum_{n=1}^{N} E_{n} (w)

この誤差関数を最小化していくのが、オンライン勾配降下法である。

オンライン勾配降下法は、確率的勾配降下法（StochasticGradientDescent）と呼ばれている。

SGDの話は、昔ちょっと触れた。

確率的勾配降下法（SGD）の並列化について - 分からんこと多すぎ

w_{(n+1)} = w_{(n)} - \eta \nabla E_{n} (w_{(n)})

要するに、上式のようなパラメータの更新をすれば良い。

Back propagation

上図のような関係を考える。

y^{(l+1)}は、y^{(l)}を重みw^{(l)}で足しあわせて、関数f^(l+1)をかけたものである。

y_{k}^{(l+1)} = f^{(l+1)} ( \sum_{j} w_{kj}^{(l)} y_{j}^{(l)} )

この時、n番目のデータに関する誤差関数E_{n}^{(l)}の重みパラメータw^{(l)}による微分は、以下のようになる。

（(l)層での出力ベクトル（シグモイド関数後の値）をy^{(l)}と書くこととする）


（導出は長いので最後に書いた）

f:id:rishida:20140218211520p:plain

ちなみにNは層の総数。

したがってパラメータの更新式は、以下のようになる。

w_{(n+1)}^{(l)} = w_{(n)}^{(l)} - \eta \nabla E_{n} (w_{(n)}^{(l)}) = w_{(n)}^{(l)} - \eta \delta^{(l+1)} {y^{(l)}}^{T}

実装

Rで実装する。

#使い方例
if(FALSE){
  w <- NULL
  for(i in 1:500){
    x <- runif(min = -3, max = 3, n = 100)
    y <- sin(x)
    if( (i < 10) || (490 < i)){
      w <- neural.net(x = matrix(x,nrow = 1), y = matrix(y, nrow=1) , w = w, graph = 'y = sin(x)', noise = 0.1)
    }
    w <- neural.net(x = matrix(x,nrow = 1), y = matrix(y, nrow=1) , w = w, graph = NULL, noise = 0.1)
  }
}

#ニューラルネット関数
#必須項目：
#入力行列x（入力の次元数=行数）,出力行列y,重みベクトルのリストw
#オプション：
#中間層素子数nnum,更新重みeta,クラス判別問題か否かのフラグCLASS,
#グラフを書くときのタイトルgraph,ノイズの分散noise
neural.net <- function(x,y,w = NULL, nnum = 5, eta = 0.1, CLASS = FALSE, graph = NULL, noise = 0){
  sigmoid <- function(x){
    return(( 1 / ( 1 + exp( -x))))
  }

  #create weight
  if(is.null(w)){
    w <- list(NULL)
    w[[1]] <- matrix(rnorm(n=nnum * (dim(x)[1]+1),mean=0,sd=(1/sqrt(dim(x)[1]))),nrow=nnum,ncol=(dim(x)[1]+1) )
    w[[2]] <- matrix(rnorm(n=dim(y)[1] * (nnum+1),mean=0,sd=(1/sqrt(dim(x)[1]))),nrow=dim(y)[1],ncol=(nnum+1))
  }

  #mix data
  idx <- sample(c(1:dim(x)[2]),size=dim(x)[2],replace=FALSE)
  x <- matrix(x[,idx],ncol=length(idx))
  y <- matrix(y[,idx],ncol=length(idx))

  #add noise
  y <- y + rnorm(mean=0,sd=noise,n=prod(dim(y)))

  for(i in 1:dim(x)[2]){
    #compute neural.net
    out <- list(NULL)
    out[[1]] <- rbind(1,x[,i])
    out[[2]] <- w[[1]] %*% out[[1]]
    out[[2]] <- rbind(1,sigmoid(out[[2]]))
    out[[3]] <- w[[2]] %*% out[[2]]
    if(CLASS == TRUE){
      out[[3]] <- sigmoid(out[[3]])
    }

    #Back Propagate
    #BP1
    delta <- -1 * (y[,i] - out[[3]]) * (sigmoid(out[[3]])) #* (1 - sigmoid(out[[3]])))
    wdet <- t(out[[2]] %*% t(delta))
    hoge <- y[,i] + (y[,i] - out[[3]])
    #eta <- (matrix(c(w[[2]] * wdet) * c(out[[2]]),nrow=dim(y)[1]) - (wdet * t(out[[2]] %*% out[[3]]))) / (wdet^2)
    w[[2]] <- w[[2]] - eta * wdet
    #BP2
    delta <- ((out[[2]][-1] * (1 - out[[2]][-1])) * as.vector(w[[2]][,-1] * delta))
    w[[1]] <- w[[1]] - eta * t(out[[1]] %*% t(delta))

    #plot
    if(!is.null(graph)){
      if( ((i - 1) %% 10) == 0 ){
        sim.neural(x,y,w,CLASS,graph)
      }
    }
  }

  print(paste('err',out[[3]] - y[,dim(x)[2]]))
  return(w)
}

#推定結果を図示する関数
sim.neural <- function(x,y,w,CLASS = FALSE,graph){
  sigmoid <- function(x){
    return( 1 / ( 1 + exp( -x)))
  }

  #sort
  idx <- order(x)
  x <- matrix(x[,idx],ncol=length(idx))
  y <- matrix(y[,idx],ncol=length(idx))

  #simurate
  out <- list(NULL)
  out[[1]] <- w[[1]] %*% rbind(1,x)
  out[[1]] <- sigmoid(out[[1]])
  out[[2]] <- w[[2]] %*% rbind(1,out[[1]])
  if(CLASS == TRUE){
    out[[2]] <- sigmoid(out[[2]])
  }

  #plot
  plot(x,y,xlim=c(min(x)-0.1,max(x)+0.1),ylim=c(min(y,out[[2]])-0.1,max(y,out[[2]])+0.1),xlab="",ylab="",type='l')
  par(new=TRUE)
  plot(x,out[[2]],col='red',xlim=c(min(x)-0.1,max(x)+0.1),ylim=c(min(y,out[[2]])-0.1,max(y,out[[2]])+0.1),xlab="x",ylab="y",main=graph)
  return(0)
}

飛ばしたBackPropagationの計算。


次回はボルツマンマシンの話。

Rのxtableにおける浮動小数点表記（floating point expression, scientific notation）

R

Rでは，行列をlatexの表形式にして出力してくれるライブラリが存在する．

この表に出力する時に，浮動小数点表記（floating point expression）をする方法が日本語で書かれていなかったため，記事にする．

library('xtable')
data <- matrix(0.0000001,3,2)
xtable(data, display=c('e','e','e'))

% latex table generated in R 2.15.3 by xtable 1.7-1 package
% Wed Jan 29 21:39:10 2014
\begin{table}[ht]
\centering
\begin{tabular}{rrr}
  \hline
 & 1 & 2 \\
  \hline
1 & 1.00e-07 & 1.00e-07 \\
  2 & 1.00e-07 & 1.00e-07 \\
  3 & 1.00e-07 & 1.00e-07 \\
   \hline
\end{tabular}
\end{table}

\begin{tabular}{rrr} \hline & 1 & 2 \\ \hline 1 & 1.00e-07 & 1.00e-07 \\ 2 & 1.00e-07 & 1.00e-07 \\ 3 & 1.00e-07 & 1.00e-07 \\ \hline \end{tabular}

ちなみに，行番号や列名を消すには，

library('xtable')
data <- matrix(0.0000001,3,2)
> print(xtable(data, display=c('e','e','e')),include.rownames=FALSE,include.colnames=FALSE)

% latex table generated in R 2.15.3 by xtable 1.7-1 package
% Wed Jan 29 21:40:35 2014
\begin{table}[ht]
\centering
\begin{tabular}{rr}
  \hline
  \hline
1.00e-07 & 1.00e-07 \\
  1.00e-07 & 1.00e-07 \\
  1.00e-07 & 1.00e-07 \\
   \hline
\end{tabular}
\end{table}

\begin{tabular}{rr} \hline \\ \hline 1.00e-07 & 1.00e-07 \\ 1.00e-07 & 1.00e-07 \\ 1.00e-07 & 1.00e-07 \\ \hline \end{tabular}