ロジスティック回帰をRで実装してみた

作者: C. M.ビショップ,元田浩,栗田多喜夫,樋口知之,松本裕治,村田昇
出版社/メーカー: シュプリンガー・ジャパン株式会社
発売日: 2007/12/10
メディア: 単行本
購入: 18人クリック: 1,588回
この商品を含むブログ (111件) を見る

この本。
一般的には「PRML」って呼ばれてる。かなりの良書（らしい）。
ある程度の知識とか経験がないと内容はわけわかめｗ
読むのにずいぶん苦労した・・・
上巻第4章の4.3.2.ロジスティック回帰〜4.3.3.反復再重み付け最小二乗(pp204-207)をRで実装してみた。

今までの分類器だと出力は2値(2クラスの場合)だったけど、今回の出力は実数、確率を出力します。
写像関数 $\phi$ で入力を多次元に写像する。（やり方わからないので実装してません。つまり識別面は線形）
ロジスティック関数
$\sigma (a) = \frac{1}{1+exp(-a)}$
これを使うと出力が0〜1の実数となるので、それを確率とすればよい。
予測関数
$p ( C_{1} | \phi ) = y ( \phi ) = \sigma ( w^{T} \phi )$
これで入力データのクラス1の確率がわかる。
クラス２の確率は次の式から求めることができる。
$p ( C_{2} | \phi ) = 1 - y( \phi )$
重みベクトル $w$ の求め方はWidow-Hoffとかと同じ感じ。
次の交差エントロピー誤差関数を最小すればいいだけ。
$E(w) = - \ln p(t|w)= - \sum_{n=1}^{N} \{ t_{n} \ln y_{n} + (1-t_{n}) \ln (1-y_{n}) \}$
これを微分して勾配ベクトルとヘシアン行列を求める。
うれしいことにベクトル・行列表現された式がある。
だからRでの実装が簡単！?
$\nabla E(w) = \sum_{n=1}^{N} (y_{n} - t_{n})\phi_{n} = \Phi^{T}( {\bf y} - {\bf t})$
$H=\nabla \nabla E(w) = \sum_{n=1}^{N} y_{n} ( 1 - y_{n} ) \phi_{n} \phi_{n}^{T} = \Phi^{T} R \Phi$
$R = diag( y_{n}(1-y_{n}) )$
勾配ベクトルとヘシアン行列がわかったので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)) )
}

誤差関数の計算は関数にした方が良かったかな？
やっかいな点が２つある。
１つは誤差関数を計算するときにlog(0)が現れること。
Rの場合、log(0)の返り値は「-Inf」とかで値ではない。
このまま計算すると「nan」とかなっちゃう。
で、そのエラー処理が必要。
もっとこー、きれいな逃げ方はないものか orz

２つ目のやっかいな点は逆行列を求めること。
ヘシアン行列が特異行列とかなっちゃうと求められない。
これも何か良い逃げ方はないものか・・・

とりあえず実行してみる。
まずは簡単な「testData0.csv」で試す。
出力が0.5より大きければクラス１、0.5以下ならクラス２とする。

#プロット関数 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)
}

まぁ、普通・・・
出力が確率なので、そのまま色で表現してみる。
白っぽいほどクラス１の確率が高くて、緑っぽいほどクラス２の確率が高い。

#プロット関数その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値出力版

まぁ、妥当な識別面。
続いて確率出力を色で表現した版

(ﾉ´▽`)ノオオオォォォォ♪
こういうのを待っていたんだよ！！ｗ

あとは非線形の仕方さえわかれば・・・