learning-hmm 0.1.0.0 → 0.1.1.0

raw patch · 7 files changed

+198/−47 lines, 7 filesdep +containersdep +random-sourcePVP ok

version bump matches the API change (PVP)

Dependencies added: containers, random-source

API changes (from Hackage documentation)

Files

• CHANGES.md +9/−0
• learning-hmm.cabal +6/−1
• src/Data/Random/Distribution/Simplex.hs +50/−0
• src/Data/Random/Distribution/Uniform/Util.hs +17/−0
• src/Learning/HMM.hs +57/−15
• src/Learning/HMM/Internal.hs +56/−31
• tests/doctests.hs +3/−0

@@ -0,0 +1,9 @@+Revision history for Haskell package learning-hmm+===++## Version 0.1.1.0+- Add function `init` for random initialization+- Add function `simulate` for running a Markov process++## Version 0.1.0.0+- Original version

@@ -1,5 +1,5 @@ name:                learning-hmm-version:             0.1.0.0+version:             0.1.1.0 stability:           experimental  synopsis:            Yet another library for hidden Markov models@@ -17,6 +17,7 @@  cabal-version:       >=1.10 build-type:          Simple+extra-source-files:  CHANGES.md  source-repository head   type:              git@@ -25,13 +26,17 @@ library   exposed-modules:   Learning.HMM   other-modules:     Data.Random.Distribution.Categorical.Util+                   , Data.Random.Distribution.Simplex+                   , Data.Random.Distribution.Uniform.Util                    , Data.Vector.Util                    , Data.Vector.Util.LinearAlgebra                    , Learning.HMM.Internal   -- other-extensions:     build-depends:     base >=4.7 && <4.8+                   , containers                    , logfloat                    , random-fu+                   , random-source                    , vector   hs-source-dirs:    src   default-language:  Haskell2010

@@ -0,0 +1,50 @@+{-# LANGUAGE+    MultiParamTypeClasses,+    FlexibleContexts, FlexibleInstances,+    UndecidableInstances, GADTs+  #-}++module Data.Random.Distribution.Simplex+    ( StdSimplex(..)+    , stdSimplex+    , stdSimplexT+    , fractionalStdSimplex+    ) where++import Control.Applicative+import Control.Monad+import Data.List+import Data.Random.RVar+import Data.Random.Distribution+import Data.Random.Distribution.Uniform++-- |Uniform distribution over a standard simplex.+newtype StdSimplex as =+    -- | @StdSimplex k@ constructs a standard simplex of dimension @k@+    -- (standard /k/-simplex).+    -- An element of the simplex represents a vector variable @as = (a_0,+    -- a_1, ..., a_k)@. The elements of @as@ are more than or equal to @0@+    -- and @sum as@ is always equal to @1@.+    StdSimplex Int+    deriving (Eq, Show)++instance (Ord a, Fractional a, Distribution StdUniform a) => Distribution StdSimplex [a] where+    rvar (StdSimplex k) = fractionalStdSimplex k++-- |@stdSimplex k@ returns a random variable being uniformly distributed over+-- a standard simplex of dimension @k@.+stdSimplex :: Distribution StdSimplex [a] => Int -> RVar [a]+stdSimplex k = rvar (StdSimplex k)++stdSimplexT :: Distribution StdSimplex [a] => Int -> RVarT m [a]+stdSimplexT k = rvarT (StdSimplex k)++-- |An algorithm proposed by Rubinstein & Melamed (1998).+-- See, /e.g./, S. Onn, I. Weissman.+-- Generating uniform random vectors over a simplex with implications to+-- the volume of a certain polytope and to multivariate extremes.+-- /Ann Oper Res/ (2011) __189__:331-342.+fractionalStdSimplex :: (Ord a, Fractional a, Distribution StdUniform a) => Int -> RVar [a]+fractionalStdSimplex k = do us <- replicateM k stdUniform+                            let us' = sort us ++ [1]+                            return $ zipWith (-) us' (0 : us')

@@ -0,0 +1,17 @@+{-# LANGUAGE FlexibleContexts, FlexibleInstances, MultiParamTypeClasses #-}++module Data.Random.Distribution.Uniform.Util () where++import Control.Applicative ((<$>))+import Data.Number.LogFloat (LogFloat, logFloat, fromLogFloat)+import Data.Random.Distribution (Distribution)+import Data.Random.Distribution.Uniform -- (StdUniform(..), Uniform(..), doubleUniform)+import Data.Random (rvarT)+import Data.Random.Source (getRandomDouble)++instance Distribution Uniform LogFloat where+  rvarT (Uniform a b) = do x <- doubleUniform (fromLogFloat a) (fromLogFloat b)+                           return $ logFloat x++instance Distribution StdUniform LogFloat where+  rvarT _ = logFloat <$> getRandomDouble

@@ -2,17 +2,24 @@     HMM (..)   , LogLikelihood   , new+  , init+  , withEmission   , viterbi   , baumWelch+  , simulate   ) where +import Prelude hiding (init)+import Control.Applicative ((<$>)) import Control.Arrow ((***), first)-import Data.Random.Distribution (pdf)+import Data.Random.Distribution (pdf, rvar) 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.Random.RVar (RVar)+import Data.Random.Sample (sample) import Data.List (genericLength) import Data.Number.LogFloat (fromLogFloat, logFloat, logFromLogFloat) import Data.Vector ((!))@@ -23,11 +30,11 @@ type LogLikelihood = Double  -- | Parameter set of the hidden Markov model. Direct use of the---   constructor is not recommended. Instead, call 'new'.+--   constructor is not recommended. Instead, call 'new' or 'init'. data HMM s o = HMM { states  :: [s] -- ^ Hidden states-                   , outputs :: [o] -- ^ Observed outputs+                   , outputs :: [o] -- ^ Outputs                    , initialStateDist :: Categorical Double s-                     -- ^ Categorical distribusion of initial states+                     -- ^ Categorical distribution of initial states                    , transitionDist :: s -> Categorical Double s                      -- ^ Categorical distribution of next states                      --   conditioned by the previous states@@ -53,11 +60,14 @@     w   = transitionDist hmm     phi = emissionDist hmm --- | Construct a 'HMM' from the given states and outputs. The---   'initialStateDist' and 'emissionDist' are set to be uniform+-- | @new states outputs@ returns a model from the @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 [1, 2 :: Int] ['C', 'D']+--   HMM {states = [1,2], outputs = "CD", initialStateDist = fromList [(0.5,1),(0.5,2)], transitionDist = [(fromList [(0.75,1),(0.25,2)],1),(fromList [(0.25,1),(0.75,2)],2)], emissionDist = [(fromList [(0.5,'C'),(0.5,'D')],1),(fromList [(0.5,'C'),(0.5,'D')],2)]} new :: (Ord s, Ord o) => [s] -> [o] -> HMM s o new ss os = HMM { states           = ss                 , outputs          = os@@ -76,8 +86,24 @@     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.+-- | @init states outputs@ returns a random variable of the model with+--   @states@ and @outputs@, wherein parameters are sampled from uniform+--   distributions.+init :: (Ord s, Ord o) => [s] -> [o] -> RVar (HMM s o)+init ss os = do hmm' <- init' (V.fromList ss) (V.fromList os)+                return $ fromHMM' hmm'++-- | @model \`withEmission\` xs@ returns a model in which the+--   'emissionDist' is updated by using the observed outputs @xs@. The+--   'emissionDist' is set to be normalized histograms each of which is+--   calculated from a partial set of @xs@ for each state. The partition is+--   based on the most likely state path obtained by the Viterbi algorithm.+withEmission :: (Ord s, Ord o) => HMM s o -> [o] -> HMM s o+withEmission model xs = fromHMM' $ withEmission' (toHMM' model) (V.fromList xs)++-- | @viterbi model xs@ performs the Viterbi algorithm using the observed+--   outputs @xs@, and returns 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`@@ -87,9 +113,9 @@     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 model xs@ performs the Baum-Welch algorithm using the+--   observed outputs @xs@, and iteratively returns a list of updated+--   models 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`@@ -99,7 +125,23 @@     model' = toHMM' model     xs'    = V.fromList xs --- | Check if the 'HMM' is valid in the sense of whether the 'states' and+-- | @simulate model t@ generates a Markov process of length @t@ using the+--   @model@, and returns its state path and observed outputs.+simulate :: HMM s o -> Int -> RVar ([s], [o])+simulate model step | step < 1  = return ([], [])+                    | otherwise = do s0 <- sample $ rvar pi0+                                     x0 <- sample $ rvar $ phi s0+                                     unzip . ((s0, x0) :) <$> sim s0 (step - 1)+  where+    sim _ 0 = return []+    sim s t = do s' <- sample $ rvar $ w s+                 x' <- sample $ rvar $ phi s+                 ((s', x') :) <$> sim s' (t - 1)+    pi0 = initialStateDist model+    w   = transitionDist model+    phi = emissionDist model++-- | Check if the model is valid in the sense of whether the 'states' and --   'outputs' are not empty. checkModelIn :: String -> HMM s o -> () checkModelIn fun hmm@@ -111,8 +153,8 @@     os = outputs hmm     err = errorIn fun --- | Check if all the elements of the data are contained in the 'outputs'---   of the 'HMM'.+-- | Check if all the elements of the observed outputs are contained in the+--   'outputs' of the model. checkDataIn :: Eq o => String -> HMM s o -> [o] -> () checkDataIn fun hmm xs   | all (`elem` os) xs = ()@@ -143,7 +185,7 @@     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'',+-- | Convert 'HMM' to 'HMM''. The 'initialStateDist'', 'transitionDist'', --   and 'emissionDistT'' are normalized. toHMM' :: (Eq s, Eq o) => HMM s o -> HMM' s o toHMM' hmm = HMM' { states'           = V.fromList ss

@@ -2,6 +2,8 @@     HMM' (..)   , Likelihood   , Probability+  , init'+  , withEmission'   , viterbi'   , baumWelch'   -- , baumWelch1'@@ -9,15 +11,20 @@   -- , backward'   ) where -import Control.Monad (forM_)+import Control.Applicative ((<$>))+import Control.Monad (forM_, replicateM) import Control.Monad.ST (runST)+import qualified Data.Map.Strict as M (empty, insertWith, findWithDefault) import Data.Number.LogFloat (LogFloat, logFloat)+import Data.Random.RVar (RVar)+import Data.Random.Distribution.Simplex (stdSimplex)+import Data.Random.Distribution.Uniform.Util () 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+    filter, foldl', foldl1', freeze, fromList, 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.Mutable as MV (new, read, write) import qualified Data.Vector.Util as V (unsafeElemIndex) import Data.Vector.Util.LinearAlgebra (     (>+>), (>.>), (>/>), (#+#), (.>), (>/), (#/), (<.>), (#.>), (<.#)@@ -27,7 +34,7 @@ type Likelihood  = LogFloat type Probability = LogFloat --- | More efficient data structure of the HMM parameters. This should be+-- | More efficient data structure of the 'HMM' model. 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@@ -37,19 +44,41 @@                      , emissionDistT'    :: Vector (Vector Probability)                      } --- | Perform the Viterbi algorithm and return the most likely state path---   and its likelihood.+init' :: Vector s -> Vector o -> RVar (HMM' s o)+init' ss os = do+  let n = V.length ss+      m = V.length os+  pi0 <- V.fromList <$> stdSimplex (n-1)+  w   <- V.fromList <$> replicateM n (V.fromList <$> stdSimplex (n-1))+  phi <- V.fromList <$> replicateM n (V.fromList <$> stdSimplex (m-1))+  return HMM' { states'           = ss+              , outputs'          = os+              , initialStateDist' = pi0+              , transitionDist'   = w+              , emissionDistT'    = V.transpose phi+              }++withEmission' :: (Ord s, Ord o) => HMM' s o -> Vector o -> HMM' s o+withEmission' model xs = model { emissionDistT' = phi' }+  where+    ss = states' model+    os = outputs' model+    (path, _) = viterbi' model xs+    mp    = V.foldl' (\m k -> M.insertWith (+) k 1 m) M.empty $ V.zip path xs+    hists = V.map (\s -> V.map (\o -> M.findWithDefault 0 (s, o) mp) os) ss+    phi'  = V.transpose $ V.map (\h -> h >/ V.sum h) hists+ 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)+      ix <- MV.new n+      ix `MV.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+        j <- ix `MV.read` (i+1)+        ix `MV.write` i $ psis ! (i+1) ! j       V.freeze ix       where         ss = states' model@@ -58,15 +87,15 @@     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)+      ds <- MV.new n+      ps <- MV.new n+      ds `MV.write` 0 $ (phi' ! x 0) >.> pi0+      ps `MV.write` 0 $ V.replicate k (0 :: Int)       forM_ [1..(n-1)] $ \i -> do-        d <- ds `M.read` (i-1)+        d <- ds `MV.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 `MV.write` i $ phi' ! x i >.> V.map V.maximum dws+        ps `MV.write` i $ V.map V.maxIndex dws       ds' <- V.freeze ds       ps' <- V.freeze ps       return (ds', ps')@@ -82,15 +111,13 @@     -- 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.+--   model and the likelihood of the old model. baumWelch1' :: (Eq s, Eq o) => HMM' s o -> Vector o -> (HMM' s o, Likelihood) baumWelch1' model xs = (model', likelihood)   where@@ -137,14 +164,13 @@     -- 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+  v <- MV.new n+  v `MV.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)+    a <- v `MV.read` (i-1)+    v `MV.write` i $ (phi' ! x i) >.> (a <.# w)   V.freeze v   where     n   = V.length xs@@ -155,14 +181,13 @@     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)+  v <- MV.new n+  v `MV.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)+    b <- v `MV.read` (i+1)+    v `MV.write` i $ w #.> ((phi' ! x (i+1)) >.> b)   V.freeze v   where     n   = V.length xs

@@ -1,6 +1,9 @@+module Main (main) where+ import Test.DocTest  main :: IO () main = doctest [ "-isrc"                , "src/Data/Vector/Util/LinearAlgebra.hs"+               , "src/Learning/HMM.hs"                ]