日々の流転

λ. 自動微分でロジスティック回帰

「自動微分を使って線形回帰をしてみる」と「nonlinear-optimizationパッケージで遊ぶ」で線形回帰をして遊んでみたので、次は予定通り自動微分でロジスティック回帰を試してみることに。 例として何を使おうか悩んでいたのだけれど、たまたま ロジスティック回帰で分類を試す - Negative/Positive Thinking (2013-11-29) という記事を（@chezouさんのtweetで）見かけたので、これと同じデータで試してみることにした。

使用するデータ

というわけで、使用するデータはこれ。

• LIBSVMのページにあるUCIデータセットのa9a
  • 学習データ : a9a
  • テストデータ : a9a.t

コード

実装にはnonlinear-optimizationパッケージを、自動微分を行うadパッケージと組み合わせて便利に使うための、拙作のnonlinear-optimization-adパッケージを使用。

コスト関数としては Courseraの Machine Learning のコースで扱った形の、正則化なしのコスト関数と、L2正則化を行うコスト関数の場合の両方を試す。

import Control.Monad
import Data.Map (Map)
import qualified Data.Map as Map
import Data.Set (Set)
import qualified Data.Set as Set
import qualified Numeric.Optimization.Algorithms.HagerZhang05.AD as CG
import Text.Printf

main :: IO ()
main = do
  trainingData <- liftM parse $ readFile "a9a"
  testData     <- liftM parse $ readFile "a9a.t"

  let features :: Set String
      features = Set.unions [Map.keysSet fs | (_, fs) <- trainingData]

      cost, costReg :: RealFloat a => Map String a -> a
      cost theta    = sum [ if y then - log p else - log (1 - p)
                          | (y,x) <- trainingData, let p = h theta (fmap realToFrac x) ]
                      / fromIntegral (length trainingData)
      costReg theta = cost theta + (lambda / (2 * fromIntegral (length trainingData))) * sum [x**2 | (feat,x) <- Map.toList theta, feat/=""]

      params = CG.defaultParameters
               { CG.printFinal  = True
               , CG.printParams = True
               , CG.verbose     = CG.Verbose
               , CG.maxItersFac = 50 / fromIntegral (1 + Set.size features)
               }
      grad_tol = 0.00001
      theta0 = Map.fromAscList [(feat,0) | feat <- Set.toAscList (Set.insert "" features)]

  (theta1, _result1, _stat1) <- CG.optimize params grad_tol theta0 cost
  (theta2, _result2, _stat2) <- CG.optimize params grad_tol theta0 costReg

  forM_ [("without regularization", theta1), ("with regularization", theta2)] $ \(s, theta) -> do
    printf "[%s]\n" s
    forM_ [("training set", trainingData), ("test set", testData)] $ \(name, samples) -> do
      let correct  = filter (\(y,x) -> predict theta x == y) samples
          accuracy :: Double
          accuracy = fromIntegral (length correct) / fromIntegral (length samples)
      printf "  %s accuracy: %d/%d = %f %%\n" name (length correct) (length samples) (100*accuracy)

-- regularization parameter
lambda :: Num a => a
lambda = 1

-- hypothesis
h :: RealFloat a => Map String a -> (Map String a -> a)
h theta x = sigmoid $ sum [(theta Map.! feat) * val | (feat, val) <- ("", 1) : Map.toList x]

predict :: Map String Double -> (Map String Double -> Bool)
predict theta x = h theta x > 0.5

sigmoid :: RealFloat a => a -> a
sigmoid x = 1 / (1 + exp (-x))

parse :: String -> [(Bool, Map String Double)]
parse = map f . lines
  where
    f l = (readInt w > 0, Map.fromList [(name, read val) | s <- ws, let (name,_:val) = break (':'==) s])
      where
        (w:ws) = words l
        readInt :: String -> Int
        readInt ('+':s) = read s
        readInt s = read s

ちょっと試してみたところなかなか収束しなかったので、参照したブログと同じくイテレーション回数50回で打ち切ることに。これは、nonlinear-optimizationパッケージのoptimize関数は、イテレーション回数の指定がパラメータ数に対する倍率になっていて、ちょっと面倒だった。

結果

結果は以下のとおりで、参照したブログの数字ともほぼ一緒だし、まあこんなものだろう。

[without regularization]
  training set accuracy: 27644/32561 = 84.89911243512176 %
  test set accuracy: 13846/16281 = 85.04391622136232 %
[with regularization]
  training set accuracy: 27648/32561 = 84.91139707011456 %
  test set accuracy: 13850/16281 = 85.06848473680978 %