ロジスティック回帰をRで実装してみた
- 作者: C. M.ビショップ,元田浩,栗田多喜夫,樋口知之,松本裕治,村田昇
- 出版社/メーカー: シュプリンガー・ジャパン株式会社
- 発売日: 2007/12/10
- メディア: 単行本
- 購入: 18人 クリック: 1,588回
- この商品を含むブログ (111件) を見る
一般的には「PRML」って呼ばれてる。かなりの良書(らしい)。
ある程度の知識とか経験がないと内容はわけわかめw
読むのにずいぶん苦労した・・・
上巻第4章の4.3.2.ロジスティック回帰〜4.3.3.反復再重み付け最小二乗(pp204-207)をRで実装してみた。
今までの分類器だと出力は2値(2クラスの場合)だったけど、今回の出力は実数、確率を出力します。
写像関数で入力を多次元に写像する。(やり方わからないので実装してません。つまり識別面は線形)
ロジスティック関数
これを使うと出力が0〜1の実数となるので、それを確率とすればよい。
予測関数
これで入力データのクラス1の確率がわかる。
クラス2の確率は次の式から求めることができる。
重みベクトルの求め方はWidow-Hoffとかと同じ感じ。
次の交差エントロピー誤差関数を最小すればいいだけ。
これを微分して勾配ベクトルとヘシアン行列を求める。
うれしいことにベクトル・行列表現された式がある。
だからRでの実装が簡単!?
勾配ベクトルとヘシアン行列がわかったのでIRLSでwを更新していく。
IRLSの中身はニュートン法をぐるぐる繰り返しているだけ。
Rでの実装
#ロジスティック回帰(2クラス) # data:学習データ, t:学習データのクラスベクトル myLoji.train <- function(data,t){ # 学習用定数の設定 N <- nrow(data) # 学習データ数 D <- ncol(data) # 学習データの次元数 MAXITS <- 500 # 学習の最大繰り返し回数 IRLS.MAXITS <- 25 # IRLSの最大繰り返し回数 TOL <- 1e-8 # IRLSの終了条件の設定 LAMDA.MIN <- 2^(-8) # IRLSの更新時の値の差の最小値 t <- replace(t, t>= 2, 0 ) # クラスラベルの設定 # クラスラベルは0か1 # 学習データに0次元目(=1)を加える(biasの計算のため) data2 <- matrix( 1, N, D+1) for( i in 1:N ){ for( j in 1:D ){ data2[i,j+1] = data[i,j] } } data <- data2 data.t <- t(data) # 転置行列 mySigmoid <- function(x){ # ロジスティックシグモイド関数 return ( 1.0/(1+exp(-x)) ) } myLoji.calc.y <- function(data.m,w.v){ # 予測値をベクトルで返す y.m <- data.m %*% w.v # N行1列の行列になるはず y.m <- apply(y.m, 2, mySigmoid) # N行1列の行列が返される return (y.m[,1]) # 1列目だけをベクトルとして返す } w.v <- numeric(D+1) # 重みベクトルwの作成(初期値 0) # IRLSの準備-------------------------------------------------- y.v <- myLoji.calc.y( data, w.v) # err.old : 最小化する関数E(w)のこと ( E(w)=PRML式(4.90) ) # err.old <- (-1) * sum( (t)*(log(y.v)) + (1-t)*(log(1-y.v)) ) # log(0)に対応させるために泣く泣く計算を分解 # log(0) = -Inf data.term1 <- (t)*(log(y.v)) # ( 1 or 0 )*log(0)=NaN data.term1 <- ifelse( is.nan(data.term1), 0, data.term1) # Nan -> 0 とする data.term2 <- (1-t)*(log(1-y.v)) data.term2 <- ifelse( is.nan(data.term2), 0, data.term2) err.old <- (-1)*( sum(data.term1 + data.term2) ) # IRLS開始(wの更新)------------------------------------------- for( IRLS.its in 1:IRLS.MAXITS ){ rn.v <- y.v*(1-y.v) R.m <- diag( (rn.v) , N, N ) e.v <- (y.v)-t # Gradient 勾配ベクトルを求める(∇E(w)=PRML式(4.96)) Gradient.v <- data.t %*% e.v Gradient.norm <- sum(Gradient.v^2) if( IRLS.its != 0 && Gradient.norm < TOL){ # 勾配の大きさが小さければ break # IRLSを終了する } # Hessian ヘシアン行列を求める(∇∇E(w)=H=PRML式(7.111)) Hessian.m <- ( (data.t)%*%(R.m) ) %*% (data) # -Δw=( (∇∇J)^(-1) * ∇J )の計算 # QR分解 連立一次方程式 Ax=b の解xを求める delta.w.v <- ( qr.solve( Hessian.m, Gradient.v)[,1] ) # ニュートン法(式の最小化):w.new = w.old + λ*(-Δw) lamda <- 1.0 while( lamda > LAMDA.MIN ){ # 最小値より大きい間繰り返す w.new.v <- (w.v)-( (lamda)*(delta.w.v) ) # 新しいwの計算 y.new.v <- myLoji.calc.y( data, w.new.v ) # 新しいyの計算 # err.new : 最小化する関数E(w)のこと ( E(w)=PRML式(4.90) ) # err.new <- (-1) * sum( (t)*(log(y.new.v)) + (1-t)*(log(1-y.new.v)) ) # log(0)に対応させるために泣く泣く計算を分解 # log(0) = -Inf data.term1 <- (t)*(log(y.new.v)) # ( 1 or 0 )*log(0)=NaN data.term1 <- ifelse( is.nan(data.term1), 0, data.term1) # Nan -> 0 とする data.term2 <- (1-t)*(log(1-y.new.v)) data.term2 <- ifelse( is.nan(data.term2), 0, data.term2) err.new <- (-1)*( sum(data.term1 + data.term2) ) # 新しいwでwを更新するべきかの判断 if( err.new > err.old ){ #式が大きくなってしまった場合 lamda <- lamda / 2 # λを小さくしてもう一度計算 }else{ #式が小さくなってくれた場合 w.v <- w.new.v # wを更新 err.old <- err.new # 更新 break # ニュートン法を終了する } } #ニュートン法終了 } # IRLS終了(wの更新)------------------------------------------- cat("err.old=",err.old) return (w.v) } #予測関数 # w:学習結果 myLoji.predict <- function(x,w){ fai.x <- c(1,x) w.fai.x <- (t(w))%*%(fai.x) return ( 1.0/(1+exp(-w.fai.x)) ) }
誤差関数の計算は関数にした方が良かったかな?
やっかいな点が2つある。
1つは誤差関数を計算するときにlog(0)が現れること。
Rの場合、log(0)の返り値は「-Inf」とかで値ではない。
このまま計算すると「nan」とかなっちゃう。
で、そのエラー処理が必要。
もっとこー、きれいな逃げ方はないものか orz
2つ目のやっかいな点は逆行列を求めること。
ヘシアン行列が特異行列とかなっちゃうと求められない。
これも何か良い逃げ方はないものか・・・
とりあえず実行してみる。
まずは簡単な「testData0.csv」で試す。
出力が0.5より大きければクラス1、0.5以下ならクラス2とする。
#プロット関数 2次元2クラス用 # w:学習結果, data:学習データ, t:学習データのクラスラベル myLoji.plot <- function(w,data,t){ a <- 200 # プロットデータの個数,大きくすると時間がかかる tx <- seq(0,9,length=a) #tx <- seq(-3,3,length=200) ty <- tx rapMyLoji <- function(x,y){ return ( myLoji.predict(c(x,y),w) ) } z <- matrix(0,length(ty),length(tx)) for( y in 1:length(tx) ){ for( x in 1:length(ty) ){ result <- rapMyLoji(tx[x],ty[y]) result <- ifelse(result>0.5,1,2) z[x,y] <- result } } # 分類結果の領域を描画 image(tx,ty,z,col=c("#00FF80","#FFBF00"),axes=TRUE,xlab="x",ylab="y") # 学習データ点をプロット points( data[,1], data[,2], col=ifelse(t==1,"red","blue"), pch=ifelse(t==1,16,17) ) title("myLoji:[testData0.csv]") #title("myLoji:[classification.csv]") return (z) }
まぁ、普通・・・
出力が確率なので、そのまま色で表現してみる。
白っぽいほどクラス1の確率が高くて、緑っぽいほどクラス2の確率が高い。
#プロット関数その2 # 確率を色で表現 # w:学習結果, data:学習データ, t:学習データのクラスベクトル myLoji.plot2 <- function(w,data,t){ tx <- seq(0,9,length=200) #tx <- seq(-3,3,length=200) ty <- tx rapMyLoji <- function(x,y){ return ( myLoji.predict(c(x,y),w) ) } z <- matrix(0,length(ty),length(tx)) for( y in 1:length(tx) ){ for( x in 1:length(ty) ){ z[x,y] <- rapMyLoji(tx[x],ty[y]) } } # 分類結果の領域を描画 image(tx,ty,z,col=terrain.colors(100),axes=TRUE,xlab="x",ylab="y") # 学習データ点をプロット points( data[,1], data[,2], col=ifelse(t==1,"blue","red"), pch=ifelse(t==1,17,16) ) title("myLoji:[testData0.csv]") #title("myLoji:[classification.csv]") return (z) }
分布が簡単すぎてつまらん!!!
つーことでちょっと複雑なデータを使うことにする。
PRMLのサポートページから「classification.txt」を頂いてきた。
「classification.csv」で試す。
まずは2値出力版
まぁ、妥当な識別面。
続いて確率出力を色で表現した版
(ノ´▽`)ノオオオォォォォ♪
こういうのを待っていたんだよ!!w
あとは非線形の仕方さえわかれば・・・