packages feed

learning-hmm (empty) → 0.1.0.0

raw patch · 9 files changed

+580/−0 lines, 9 filesdep +basedep +doctestdep +logfloatsetup-changed

Dependencies added: base, doctest, logfloat, random-fu, vector

Files

+ LICENSE view
@@ -0,0 +1,22 @@+The MIT License (MIT)++Copyright (c) 2014 Mitsuhiro Nakamura++Permission is hereby granted, free of charge, to any person obtaining a copy+of this software and associated documentation files (the "Software"), to deal+in the Software without restriction, including without limitation the rights+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell+copies of the Software, and to permit persons to whom the Software is+furnished to do so, subject to the following conditions:++The above copyright notice and this permission notice shall be included in all+copies or substantial portions of the Software.++THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE+SOFTWARE.+
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ learning-hmm.cabal view
@@ -0,0 +1,46 @@+name:                learning-hmm+version:             0.1.0.0+stability:           experimental++synopsis:            Yet another library for hidden Markov models+description:         This library provides functions for the maximum likelihood+                     estimation of discrete hidden Markov models. At present,+                     only Baum-Welch and Viterbi algorithms are implemented.+category:            Algorithms, Machine Learning, Statistics++author:              Mitsuhiro Nakamura+maintainer:          Mitsuhiro Nakamura <m.nacamura@gmail.com>+copyright:           Copyright (c) 2014 Mitsuhiro Nakamura+license:             MIT+license-file:        LICENSE+homepage:            https://github.com/mnacamura/learning-hmm++cabal-version:       >=1.10+build-type:          Simple++source-repository head+  type:              git+  location:          https://github.com/mnacamura/learning-hmm.git++library+  exposed-modules:   Learning.HMM+  other-modules:     Data.Random.Distribution.Categorical.Util+                   , Data.Vector.Util+                   , Data.Vector.Util.LinearAlgebra+                   , Learning.HMM.Internal+  -- other-extensions:  +  build-depends:     base >=4.7 && <4.8+                   , logfloat+                   , random-fu+                   , vector+  hs-source-dirs:    src+  default-language:  Haskell2010+  ghc-options:       -Wall++test-suite doctests+  main-is:           doctests.hs+  type:              exitcode-stdio-1.0+  build-depends:     base, doctest >= 0.9.11+  hs-source-dirs:    tests+  default-language:  Haskell2010+  ghc-options:       -threaded -Wall
+ src/Data/Random/Distribution/Categorical/Util.hs view
@@ -0,0 +1,12 @@+{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-}++module Data.Random.Distribution.Categorical.Util () where++import Data.Random.Distribution (PDF, pdf)+import Data.Random.Distribution.Categorical (Categorical, toList, totalWeight)+import Data.Tuple (swap)++instance Eq a => PDF (Categorical Double) a where+  pdf cat a = case lookup a $ map swap $ toList cat of+                Nothing -> 0+                Just p  -> p / totalWeight cat
+ src/Data/Vector/Util.hs view
@@ -0,0 +1,13 @@+-- | Miscellaneous utility functions for "Data.Vector"+module Data.Vector.Util (+    unsafeElemIndex+  ) where++import Data.Maybe (fromJust)+import Data.Vector (Vector, elemIndex)++-- | Return the index of the first occurrence of the given element or throw+--   an error if no such occurrence.+{-# INLINE unsafeElemIndex #-}+unsafeElemIndex :: Eq a => a -> Vector a -> Int+unsafeElemIndex e = fromJust . elemIndex e
+ src/Data/Vector/Util/LinearAlgebra.hs view
@@ -0,0 +1,139 @@+-- | Operators commonly used in the basic linear algebra. Note that all the+--   functions defined here do not check the dimension/length of+--   vectors/matrices.+module Data.Vector.Util.LinearAlgebra (+  -- * Pairwise operators+    (>+>)+  -- , (>->)+  , (>.>)+  , (>/>)+  , (#+#)+  -- , (#-#)+  -- , (#.#)+  -- , (#/#)++  -- * Scalar-vector/vector-scalar operators+  -- , (+>)+  -- , (->)+  , (.>)+  -- , (/>)+  -- , (>+)+  -- , (>-)+  -- , (>.)+  , (>/)++  -- * Scalar-matrix/matrix-scalar operators+  -- , (+#)+  -- , (-#)+  -- , (.#)+  -- , (/#)+  -- , (#+)+  -- , (#-)+  -- , (#.)+  , (#/)++  -- * Dot and matrix-vector/vector-matrix products+  , (<.>)+  , (#.>)+  , (<.#)++  -- * Unary operators+  , transpose+  ) where++import Prelude hiding (any, head, map, null, sum, tail, zipWith)+import Data.Vector (Vector, any, cons, empty, head, map, null, sum, tail, zipWith)++-- $setup+-- >>> :module + Data.Vector++-- | Pairwise addition between two vectors+--+-- >>> fromList [1, 2] >+> fromList [3, 4 :: Int]+-- fromList [4,6]+(>+>) :: Num a => Vector a -> Vector a -> Vector a+{-# INLINE (>+>) #-}+u >+> v = zipWith (+) u v++-- | Pairwise product between two vectors+--+-- >>> fromList [1, 2] >.> fromList [3, 4 :: Double]+-- fromList [3.0,8.0]+(>.>) :: Num a => Vector a -> Vector a -> Vector a+{-# INLINE (>.>) #-}+u >.> v = zipWith (*) u v++-- | Pairwise division between two vectors+--+-- >>> fromList [1, 2] >/> fromList [3, 4 :: Double]+-- fromList [0.3333333333333333,0.5]+(>/>) :: Fractional a => Vector a -> Vector a -> Vector a+{-# INLINE (>/>) #-}+u >/> v = zipWith (/) u v++-- | Pairwise addition between two matrices+--+-- >>> fromList [fromList [1, 2], fromList [3, 4]] #+# fromList [fromList [5, 6], fromList [7, 8 :: Int]]+-- fromList [fromList [6,8],fromList [10,12]]+(#+#) :: Num a => Vector (Vector a) -> Vector (Vector a) -> Vector (Vector a)+{-# INLINE (#+#) #-}+m #+# n = zipWith (>+>) m n++-- | Scalar-vector product+--+-- >>> 2 .> fromList [1, 2 :: Integer]+-- fromList [2,4]+(.>) :: Num a => a -> Vector a -> Vector a+{-# INLINE (.>) #-}+s .> v = map (s *) v++-- | Vector-scalar division+--+-- >>> fromList [1, 2 :: Double] >/ 2+-- fromList [0.5,1.0]+(>/) :: Fractional a => Vector a -> a -> Vector a+{-# INLINE (>/) #-}+v >/ s = map (/ s) v++-- | Matrix-scalar division+--+-- >>> fromList [fromList [1, 2], fromList [3, 4 :: Double]] #/ 2+-- fromList [fromList [0.5,1.0],fromList [1.5,2.0]]+(#/) :: Fractional a => Vector (Vector a) -> a -> Vector (Vector a)+{-# INLINE (#/) #-}+m #/ s = map (>/ s) m++-- | Dot product+--+-- >>> fromList [1, 2] <.> fromList [3, 4 :: Int]+-- 11+(<.>) :: Num a => Vector a -> Vector a -> a+{-# INLINE (<.>) #-}+u <.> v = sum $ u >.> v++-- | Matrix-vector product+--+-- >>> fromList [fromList [1, 2], fromList [3, 4]] #.> fromList [1, 2 :: Double]+-- fromList [5.0,11.0]+(#.>) :: Num a => Vector (Vector a) -> Vector a -> Vector a+{-# INLINE (#.>) #-}+m #.> v = map (<.> v) m++-- | Vector-matrix product+--+-- >>> fromList [1, 2 :: Double] <.# fromList [fromList [1, 2], fromList [3, 4]]+-- fromList [7.0,10.0]+(<.#) :: Num a => Vector a -> Vector (Vector a) -> Vector a+{-# INLINE (<.#) #-}+v <.# m | any null m = empty+        | otherwise  = (v <.> map head m) `cons` (v <.# map tail m)++-- | Matrix transpose+--+-- >>> transpose $ fromList [fromList "ab", fromList "cd"]+-- fromList [fromList "ac",fromList "bd"]+transpose :: Vector (Vector a) -> Vector (Vector a)+{-# INLINE transpose #-}+transpose m+  | any null m = empty+  | otherwise  = map head m `cons` transpose (map tail m)
+ src/Learning/HMM.hs view
@@ -0,0 +1,166 @@+module Learning.HMM (+    HMM (..)+  , LogLikelihood+  , new+  , viterbi+  , baumWelch+  ) where++import Control.Arrow ((***), first)+import Data.Random.Distribution (pdf)+import Data.Random.Distribution.Categorical (Categorical)+import qualified Data.Random.Distribution.Categorical as C (+    fromList, fromWeightedList, normalizeCategoricalPs+  )+import Data.Random.Distribution.Categorical.Util ()+import Data.List (genericLength)+import Data.Number.LogFloat (fromLogFloat, logFloat, logFromLogFloat)+import Data.Vector ((!))+import qualified Data.Vector as V (elemIndex, fromList, map, toList, zip)+import qualified Data.Vector.Util.LinearAlgebra as V (transpose)+import Learning.HMM.Internal++type LogLikelihood = Double++-- | Parameter set of the hidden Markov model. Direct use of the+--   constructor is not recommended. Instead, call 'new'.+data HMM s o = HMM { states  :: [s] -- ^ Hidden states+                   , outputs :: [o] -- ^ Observed outputs+                   , initialStateDist :: Categorical Double s+                     -- ^ Categorical distribusion of initial states+                   , transitionDist :: s -> Categorical Double s+                     -- ^ Categorical distribution of next states+                     --   conditioned by the previous states+                   , emissionDist :: s -> Categorical Double o+                     -- ^ Categorical distribution of outputs conditioned+                     --   by the hidden states+                   }++instance (Show s, Show o) => Show (HMM s o) where+  show = showHMM++showHMM :: (Show s, Show o) => HMM s o -> String+showHMM hmm = "HMM {states = "           ++ show ss+              ++ ", outputs = "          ++ show os+              ++ ", initialStateDist = " ++ show pi0+              ++ ", transitionDist = "   ++ show [(w s, s) | s <- ss]+              ++ ", emissionDist = "     ++ show [(phi s, s) | s <- ss]+              ++ "}"+  where+    ss  = states hmm+    os  = outputs hmm+    pi0 = initialStateDist hmm+    w   = transitionDist hmm+    phi = emissionDist hmm++-- | Construct a 'HMM' from the given states and outputs. The+--   'initialStateDist' and 'emissionDist' are set to be uniform+--   distributions. The 'transitionDist' is specified as follows: with+--   probability 1/2, move to the same state, otherwise, move to a random+--   state (which might be the same state).+new :: (Ord s, Ord o) => [s] -> [o] -> HMM s o+new ss os = HMM { states           = ss+                , outputs          = os+                , initialStateDist = pi0+                , transitionDist   = w+                , emissionDist     = phi+                }+  where+    pi0 = C.fromWeightedList [(1, s) | s <- ss]+    w s | s `elem` ss = C.fromList [(p s', s') | s' <- ss]+        | otherwise   = C.fromList []+      where+        k = genericLength ss+        p s' | s' == s   = 1/2 * (1 + 1/k)+             | otherwise = 1/2 / k+    phi s | s `elem` ss = C.fromWeightedList [(1, o) | o <- os]+          | otherwise   = C.fromList []++-- | Perform the Viterbi algorithm and return the most likely state path+--   and its log likelihood.+viterbi :: (Eq s, Eq o) => HMM s o -> [o] -> ([s], LogLikelihood)+viterbi model xs =+  checkModelIn "viterbi" model `seq`+  checkDataIn "viterbi" model xs `seq`+  (V.toList *** logFromLogFloat) $ viterbi' model' xs'+  where+    model' = toHMM' model+    xs'    = V.fromList xs++-- | Perform the Baum-Welch algorithm steps iteratively and return+--   a list of updated 'HMM' parameters and their corresponding log+--   likelihoods.+baumWelch :: (Eq s, Eq o) => HMM s o -> [o] -> [(HMM s o, LogLikelihood)]+baumWelch model xs =+  checkModelIn "baumWelch" model `seq`+  checkDataIn "baumWelch" model xs `seq`+  map (fromHMM' *** logFromLogFloat) $ baumWelch' model' xs'+  where+    model' = toHMM' model+    xs'    = V.fromList xs++-- | Check if the 'HMM' is valid in the sense of whether the 'states' and+--   'outputs' are not empty.+checkModelIn :: String -> HMM s o -> ()+checkModelIn fun hmm+  | null ss   = err "empty states"+  | null os   = err "empty outputs"+  | otherwise = ()+  where+    ss = states hmm+    os = outputs hmm+    err = errorIn fun++-- | Check if all the elements of the data are contained in the 'outputs'+--   of the 'HMM'.+checkDataIn :: Eq o => String -> HMM s o -> [o] -> ()+checkDataIn fun hmm xs+  | all (`elem` os) xs = ()+  | otherwise          = err "illegal data"+  where+    os = outputs hmm+    err = errorIn fun++-- | Convert 'HMM'' to 'HMM'.+fromHMM' :: (Eq s, Eq o) => HMM' s o -> HMM s o+fromHMM' hmm' = HMM { states           = V.toList ss+                    , outputs          = V.toList os+                    , initialStateDist = C.fromList pi0'+                    , transitionDist   = \s -> case V.elemIndex s ss of+                                                 Nothing -> C.fromList []+                                                 Just i  -> C.fromList $ w' i+                    , emissionDist     = \s -> case V.elemIndex s ss of+                                                 Nothing -> C.fromList []+                                                 Just i  -> C.fromList $ phi' i+                    }+  where+    ss  = states' hmm'+    os  = outputs' hmm'+    pi0 = initialStateDist' hmm'+    w   = transitionDist' hmm'+    phi = V.transpose $ emissionDistT' hmm'+    pi0'   = V.toList $ V.map (first fromLogFloat) $ V.zip pi0 ss+    w' i   = V.toList $ V.map (first fromLogFloat) $ V.zip (w ! i) ss+    phi' i = V.toList $ V.map (first fromLogFloat) $ V.zip (phi ! i) os++-- | Convert 'HMM' to 'HMM''. The 'initialStateDist'', 'transisionDist'',+--   and 'emissionDistT'' are normalized.+toHMM' :: (Eq s, Eq o) => HMM s o -> HMM' s o+toHMM' hmm = HMM' { states'           = V.fromList ss+                  , outputs'          = V.fromList os+                  , initialStateDist' = V.fromList pi0'+                  , transitionDist'   = V.fromList w'+                  , emissionDistT'    = V.fromList phi'+                  }+  where+    ss  = states hmm+    os  = outputs hmm+    pi0 = C.normalizeCategoricalPs $ initialStateDist hmm+    w   = C.normalizeCategoricalPs . transitionDist hmm+    phi = C.normalizeCategoricalPs . emissionDist hmm+    pi0' = [logFloat $ pdf pi0 s | s <- ss]+    w'   = [V.fromList [logFloat $ pdf (w s) s' | s' <- ss] | s <- ss]+    phi' = [V.fromList [logFloat $ pdf (phi s) o | s <- ss] | o <- os]++errorIn :: String -> String -> a+errorIn fun msg = error $ "Learning.HMM." ++ fun ++ ": " ++ msg
+ src/Learning/HMM/Internal.hs view
@@ -0,0 +1,174 @@+module Learning.HMM.Internal (+    HMM' (..)+  , Likelihood+  , Probability+  , viterbi'+  , baumWelch'+  -- , baumWelch1'+  -- , forward'+  -- , backward'+  ) where++import Control.Monad (forM_)+import Control.Monad.ST (runST)+import Data.Number.LogFloat (LogFloat, logFloat)+import Data.Vector (Vector, (!))+import qualified Data.Vector as V (+    filter, foldl1', freeze, last, length, map, maximum, maxIndex+  , replicate, sum , tail, zip , zipWith, zipWith3, zipWith4+  )+import qualified Data.Vector.Mutable as M (new, read, write)+import qualified Data.Vector.Util as V (unsafeElemIndex)+import Data.Vector.Util.LinearAlgebra (+    (>+>), (>.>), (>/>), (#+#), (.>), (>/), (#/), (<.>), (#.>), (<.#)+  )+import qualified Data.Vector.Util.LinearAlgebra as V (transpose)++type Likelihood  = LogFloat+type Probability = LogFloat++-- | More efficient data structure of the HMM parameters. This should be+--   only used internally. The 'emissionDistT'' is a transposed matrix in+--   order to simplify the calculation.+data HMM' s o = HMM' { states'           :: Vector s+                     , outputs'          :: Vector o+                     , initialStateDist' :: Vector Probability+                     , transitionDist'   :: Vector (Vector Probability)+                     , emissionDistT'    :: Vector (Vector Probability)+                     }++-- | Perform the Viterbi algorithm and return the most likely state path+--   and its likelihood.+viterbi' :: Eq o => HMM' s o -> Vector o -> (Vector s, Likelihood)+viterbi' model xs = (path, likelihood)+  where+    -- The following procedure is based on+    -- http://ibisforest.org/index.php?cmd=read&page=Viterbi%E3%82%A2%E3%83%AB%E3%82%B4%E3%83%AA%E3%82%BA%E3%83%A0&word=Viterbi+    path = V.map (ss !) $ runST $ do+      ix <- M.new n+      ix `M.write` (n-1) $ V.maxIndex $ deltas ! (n-1)+      forM_ (reverse [0..(n-2)]) $ \i -> do+        j <- ix `M.read` (i+1)+        ix `M.write` i $ psis ! (i+1) ! j+      V.freeze ix+      where+        ss = states' model+    likelihood = V.maximum $ deltas ! (n-1)++    deltas :: Vector (Vector Probability)+    psis   :: Vector (Vector Int)+    (deltas, psis) = runST $ do+      ds <- M.new n+      ps <- M.new n+      ds `M.write` 0 $ (phi' ! x 0) >.> pi0+      ps `M.write` 0 $ V.replicate k (0 :: Int)+      forM_ [1..(n-1)] $ \i -> do+        d <- ds `M.read` (i-1)+        let dws = V.map (d >.>) w'+        ds `M.write` i $ phi' ! x i >.> V.map V.maximum dws+        ps `M.write` i $ V.map V.maxIndex dws+      ds' <- V.freeze ds+      ps' <- V.freeze ps+      return (ds', ps')+      where+        k   = V.length $ states' model+        x i = let os  = outputs' model+                  xs' = V.map (`V.unsafeElemIndex` os) xs+              in xs' ! i+        pi0  = initialStateDist' model+        w'   = V.transpose $ transitionDist' model+        phi' = emissionDistT' model++    -- Here we assumed that+    n = V.length xs++-- | Perform the Baum-Welch algorithm steps iteratively and return+--   a list of updated 'HMM'' parameters and their corresponding likelihoods.+baumWelch' :: (Eq s, Eq o) => HMM' s o -> Vector o -> [(HMM' s o, Likelihood)]+baumWelch' model xs = zip ms $ tail ells+  where+    (ms, ells) = unzip $ iterate ((`baumWelch1'` xs) . fst) (model, undefined)++-- | Perform one step of the Baum-Welch algorithm and return the updated+--   'HMM'' parameters and the likelihood of the old parameters.+baumWelch1' :: (Eq s, Eq o) => HMM' s o -> Vector o -> (HMM' s o, Likelihood)+baumWelch1' model xs = (model', likelihood)+  where+    model' = model { initialStateDist' = pi0+                   , transitionDist'   = w+                   , emissionDistT'    = phi'+                   }+    likelihood = V.last ells++    -- First, we calculate the alpha and beta values using the+    -- forward-backward algorithm.+    alphas = forward' model xs+    betas  = backward' model xs++    -- Then, we obtain the likelihoods for each time. This should be+    -- constant over time.+    ells = V.zipWith (<.>) alphas betas++    -- Based on the alpha, beta, and likelihood values, we calculate the+    -- gamma and xi values.+    gammas = V.zipWith3 (\a b l -> a >.> b >/ l) alphas betas ells+    xis    = V.zipWith4 (\a b l x -> let w1 = V.zipWith (.>) a w0+                                         w2 = V.map (phi0 ! x >.> b >.>) w1+                                     in w2 #/ l)+               alphas (V.tail betas) (V.tail ells) (V.tail xs')+      where+        xs'  = V.map (`V.unsafeElemIndex` os) xs+        w0   = transitionDist' model+        phi0 = emissionDistT' model++    -- Using the gamma and xi values, we finally obtain the optimal initial+    -- state probability vector, transition probability matrix, and+    -- emission probability matrix.+    pi0  = let gs = gammas ! 0+           in gs >/ V.sum gs+    w    = let ws = V.foldl1' (#+#) xis+               zs = V.map V.sum ws+           in V.zipWith (>/) ws zs+    phi' = let gs' o = V.map snd $ V.filter ((== o) . fst) $ V.zip xs gammas+               phis  = V.foldl1' (>+>) . gs'+               zs    = V.foldl1' (>+>) gammas+           in V.map (\o -> phis o >/> zs) os++    -- Here we assumed that+    os = outputs' model++-- | Baum-Welch forward algorithm that generates α values+forward' :: Eq o => HMM' s o -> Vector o -> Vector (Vector Probability)+forward' model xs = runST $ do+  v <- M.new n+  v `M.write` 0 $ (phi' ! x 0) >.> pi0+  forM_ [1..(n-1)] $ \i -> do+    a <- v `M.read` (i-1)+    v `M.write` i $ (phi' ! x i) >.> (a <.# w)+  V.freeze v+  where+    n   = V.length xs+    x i = let os  = outputs' model+              xs' = V.map (`V.unsafeElemIndex` os) xs+          in xs' ! i+    pi0  = initialStateDist' model+    w    = transitionDist' model+    phi' = emissionDistT' model++-- | Baum-Welch backward algorithm that generates β values+backward' :: Eq o => HMM' s o -> Vector o -> Vector (Vector Probability)+backward' model xs = runST $ do+  v <- M.new n+  v `M.write` (n-1) $ V.replicate k $ logFloat (1 :: Double)+  forM_ (reverse [0..(n-2)]) $ \i -> do+    b <- v `M.read` (i+1)+    v `M.write` i $ w #.> ((phi' ! x (i+1)) >.> b)+  V.freeze v+  where+    n   = V.length xs+    k   = V.length $ states' model+    x i = let os  = outputs' model+              xs' = V.map (`V.unsafeElemIndex` os) xs+          in xs' ! i+    w    = transitionDist' model+    phi' = emissionDistT' model
+ tests/doctests.hs view
@@ -0,0 +1,6 @@+import Test.DocTest++main :: IO ()+main = doctest [ "-isrc"+               , "src/Data/Vector/Util/LinearAlgebra.hs"+               ]