hmm (empty) → 0.1

raw patch · 5 files changed

+270/−0 lines, 5 filesdep +basedep +containerssetup-changed

Dependencies added: base, containers

Files

• Data/HMM.hs +152/−0
• Data/Lognum.hs +50/−0
• LICENSE +30/−0
• Setup.lhs +4/−0
• hmm.cabal +34/−0

@@ -0,0 +1,152 @@+module Data.HMM+    (Prob, HMM, train, bestSequence, sequenceProb)+    where++import qualified Data.Map as M+import Data.List (sort, groupBy, maximumBy)+import Data.Maybe (fromMaybe, fromJust)+import System.IO.Unsafe (unsafeInterleaveIO)+import System.Environment (getArgs)+import Control.Monad+import qualified Data.Foldable+import Debug.Trace+import Data.Lognum++type Prob = Lognum Double++-- | The type of Hidden Markov Models.+data HMM state observation = HMM [state] [Prob] [[Prob]] (observation -> [Prob])++-- | Perform a single step in the Viterbi algorithm.+--  +--   Takes a list of path probabilities, and an observation, and returns the updated+--   list of (surviving) paths with probabilities.+viterbi :: Ord observation =>+               HMM state observation+            -> [(Prob, [state])]+            -> observation+            -> [(Prob, [state])]+viterbi (HMM states _ state_transitions observations) prev x = +    [maximumBy (compare `on` fst)+            [(transition_prob * prev_prob * observation_prob,+               new_state:path)+                    | transition_prob <- transition_probs+                    | (prev_prob, path) <- prev+                    | observation_prob <- observation_probs]+        | transition_probs <- state_transitions+        | new_state <- states]+    where+        observation_probs = observations x++-- | The initial value for the Viterbi algorithm+viterbi_init :: HMM state observation -> [(Prob, [state])]+viterbi_init (HMM states state_probs _ _) = zip state_probs (map (:[]) states)++-- | Perform a single step of the forward algorithm+-- +--   Each item in the input and output list is the probability that the system+--   ended in the respective state.+forward :: Ord observation =>+               HMM state observation+            -> [Prob]+            -> observation+            -> [Prob]+forward (HMM _ _ state_transitions observations) prev x =+    [sum [transition_prob * prev_prob * observation_prob+                | transition_prob <- transition_probs+                | prev_prob <- prev+                | observation_prob <- observation_probs]+        | transition_probs <- state_transitions]+    where+        observation_probs = observations x++-- | The initial value for the forward algorithm+forward_init :: HMM state observation -> [Prob]+forward_init (HMM _ state_probs _ _) = state_probs++learn_states :: (Ord state, Fractional prob) => [(observation, state)] -> M.Map state prob+learn_states xs = histogram $ map snd xs++learn_transitions :: (Ord state, Fractional prob) => [(observation, state)] -> M.Map (state, state) prob+learn_transitions xs = let xs' = map snd xs in+                        histogram $ zip xs' (tail xs')++learn_observations ::  (Ord state, Ord observation, Fractional prob) =>+                       M.Map state prob+                    -> [(observation, state)]+                    -> M.Map (observation, state) prob+learn_observations state_prob = M.mapWithKey (\ (observation, state) prob -> prob / (fromJust $ M.lookup state state_prob))+                            . histogram++histogram :: (Ord a, Fractional prob) => [a] -> M.Map a prob+histogram xs = let hist = foldr (flip (M.insertWith (+)) 1) M.empty xs in+                M.map (/ M.fold (+) 0 hist) hist++readBrownFile :: FilePath -> IO [(String, String)]+readBrownFile = (liftM (map split . words)) . readFile+    where+        split [] = ([], [])+        split ('/':xs) = ([], xs)+        split (x:xs) = let (first, snd) = split xs in+                        (x:first, snd)++-- | Calculate the parameters of an HMM from a list of observations+--   and the corresponding states.+train :: (Ord observation, Ord state) =>+            [(observation, state)]+         -> HMM state observation+train sample = model+    where+        states = learn_states sample+        state_list = M.keys states+        +        transitions = learn_transitions sample+        trans_prob_mtx = [[fromMaybe 1e-10 $ M.lookup (old_state, new_state) transitions+                                | old_state <- state_list]+                                | new_state <- state_list]++        observations = learn_observations states sample+        observation_probs = fromMaybe (fill state_list []) . (flip M.lookup $+                            M.fromList $ map (\ (e, xs) -> (e, fill state_list xs)) $+                                map (\ xs -> (fst $ head xs, map snd xs)) $+                                groupBy     ((==) `on` fst)+                                            [(observation, (state, prob))+                                                | ((observation, state), prob) <- M.toAscList observations])++        initial = map (\ state -> (fromJust $ M.lookup state states, [state])) state_list++        model = HMM state_list (fill state_list $ M.toAscList states) trans_prob_mtx observation_probs++        fill :: Eq state => [state] -> [(state, Prob)] -> [Prob]+        fill states [] = map (const 1e-10) states+        fill (s:states) xs@((s', p):xs') = if s /= s' then+                                            1e-10 : fill states xs+                                           else+                                            p : fill states xs'++-- | Test Viterbi's algorithm on an HMM by comparing the predicted states+--   against known states for the observations.+testViterbi :: (Ord observation, Ord state) =>+                   HMM state observation+                -> [(observation, state)]+                -> Rational+testViterbi hmm testData = (fromIntegral $ length $ filter id $ zipWith (==) (bestSequence hmm observations) states) +                            / (fromIntegral $ length testData)+    where+        observations = map fst testData+        states = map snd testData++train2 :: FilePath -> IO (HMM String String)+train2 = (liftM train) . readBrownFile++-- | Calculate the most likely sequence of states for a given sequence of observations+--   using Viterbi's algorithm+bestSequence :: (Ord observation) => HMM state observation -> [observation] -> [state]+bestSequence hmm = (reverse . tail . snd . (maximumBy (compare `on` fst))) . (foldl (viterbi hmm) (viterbi_init hmm))++-- | Calculate the probability of a given sequence of observations+--    using the forward algorithm.+sequenceProb :: (Ord observation) => HMM state observation -> [observation] -> Prob+sequenceProb hmm = sum . (foldl (forward hmm) (forward_init hmm))++on f g a b = f (g a) (g b)

@@ -0,0 +1,50 @@+module Data.Lognum (Lognum)+    where++import Data.Ratio++data Lognum t = L !Int !t deriving (Eq, Ord)++fromFloating :: (Floating t, Ord t) => t -> Lognum t+fromFloating a = case a `compare` 0 of+                LT -> L (-1) (log $ -a)+                EQ -> L 0 0+                GT -> L 1 (log a)++toFloating :: (Floating t, Ord t) => Lognum t -> t+toFloating (L s m) = case s of+                        -1 -> negate $ exp m+                        0  -> 0+                        1  -> exp m++logify2 :: (Floating t, Ord t) => (t -> t -> t) -> (Lognum t -> Lognum t -> Lognum t)+logify2 (*) a b = fromFloating $ toFloating a * toFloating b++instance (Floating t, Ord t) => Show (Lognum t) where+    show = show . toFloating++instance (Floating t, Ord t) => Num (Lognum t) where+    (+) = logify2 (+)++    (L s m) * (L s' m') =   if s == 0 || s' == 0 then+                                L 0 0+                            else+                                L (s*s') (m+m')++    (-) = logify2 (-)++    negate (L s m) = L (negate s) m++    abs (L s m) = L (abs s) m++    signum (L s m) = (L s 0)++    fromInteger = fromFloating . fromInteger++instance (Floating t, Ord t) => Fractional (Lognum t) where+    _ / (L 0 _) = error "division by zero"+    (L s m) / (L s' m') = L (s*s') (m-m')+    +    recip (L s m) = L s (-m)+    +    fromRational x = (fromInteger $ numerator x) / (fromInteger $ denominator x)

@@ -0,0 +1,30 @@+Copyright (c) Max Rabkin++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions+are met:++1. Redistributions of source code must retain the above copyright+   notice, this list of conditions and the following disclaimer.++2. Redistributions in binary form must reproduce the above copyright+   notice, this list of conditions and the following disclaimer in the+   documentation and/or other materials provided with the distribution.++3. Neither the name of the author nor the names of his contributors+   may be used to endorse or promote products derived from this software+   without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE CONTRIBUTORS ``AS IS'' AND ANY EXPRESS+OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE+DISCLAIMED.  IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR+ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS+OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)+HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,+STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN+ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE+POSSIBILITY OF SUCH DAMAGE.

@@ -0,0 +1,4 @@+#! /usr/bin/env runhaskell++> import Distribution.Simple+> main = defaultMain

@@ -0,0 +1,34 @@+Name:               hmm+Version:            0.1+Description:+    A simple library for working with Hidden Markov Models.+    Should be usable even by people who are not familiar with+    HMMs. Includes implementations of Viterbi's algorithm and+    the forward algorithm.+Category:           algorithms, natural language processing, data mining+Synopsis:           Hidden Markov Model algorithms+License:            BSD3+License-file:       LICENSE+Author:             Max Rabkin+Maintainer:         max.rabkin@gmail.com+Stability:          Alpha+Build-Type:         Simple+Cabal-Version:      >= 1.2++Flag small_base+    Description: Choose the new smaller, split-up base package.++Library+    if flag(small_base)+        Build-Depends:  base >= 3, containers+    else+        Build-Depends:  base < 3+        +    Exposed-Modules:+        Data.HMM++    Other-Modules:+        Data.Lognum++    Extensions:+        ParallelListComp