learning-hmm 0.2.0.0 → 0.2.1.0
raw patch · 6 files changed
+174/−264 lines, 6 filesdep +hmatrixPVP ok
version bump matches the API change (PVP)
Dependencies added: hmatrix
API changes (from Hackage documentation)
Files
- CHANGES.md +4/−0
- learning-hmm.cabal +2/−2
- src/Data/Vector/Generic/Util/LinearAlgebra.hs +0/−141
- src/Learning/HMM.hs +35/−29
- src/Learning/HMM/Internal.hs +133/−91
- tests/doctests.hs +0/−1
CHANGES.md view
@@ -1,6 +1,10 @@ Revision history for Haskell package learning-hmm === +## Version 0.2.1.0+- Performance improvements with employing 'hmatrix'+- `withEmission` now does iterative re-estimation based on the Viterbi path+ ## Version 0.2.0.0 - Remove dependency on the 'logfloat' package - Performance improvements
learning-hmm.cabal view
@@ -1,5 +1,5 @@ name: learning-hmm-version: 0.2.0.0+version: 0.2.1.0 stability: experimental synopsis: Yet another library for hidden Markov models@@ -28,12 +28,12 @@ other-modules: Data.Random.Distribution.Categorical.Util , Data.Random.Distribution.Simplex , Data.Vector.Generic.Util- , Data.Vector.Generic.Util.LinearAlgebra , Learning.HMM.Internal -- other-extensions: build-depends: base >=4.7 && <4.8 , containers , deepseq+ , hmatrix < 0.16.1 , random-fu , random-source , vector
− src/Data/Vector/Generic/Util/LinearAlgebra.hs
@@ -1,141 +0,0 @@--- | 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.Generic.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.Generic (- Vector, any, cons, convert, 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 v a) => v a -> v a -> v 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 v a) => v a -> v a -> v 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 v a) => v a -> v a -> v 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 v a, Vector w (v a)) => w (v a) -> w (v a) -> w (v a)-{-# INLINE (#+#) #-}-m #+# n = zipWith (>+>) m n---- | Scalar-vector product------ >>> 2 .> fromList [1, 2 :: Integer]--- fromList [2,4]-(.>) :: (Num a, Vector v a) => a -> v a -> v a-{-# INLINE (.>) #-}-s .> v = map (s *) v---- | Vector-scalar division------ >>> fromList [1, 2 :: Double] >/ 2--- fromList [0.5,1.0]-(>/) :: (Fractional a, Vector v a) => v a -> a -> v 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 v a, Vector w (v a)) => w (v a) -> a -> w (v a)-{-# INLINE (#/) #-}-m #/ s = map (>/ s) m---- | Dot product------ >>> fromList [1, 2] <.> fromList [3, 4 :: Int]--- 11-(<.>) :: (Num a, Vector v a) => v a -> v 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 v a, Vector w (v a), Vector w a) => w (v a) -> v a -> v a-{-# INLINE (#.>) #-}-m #.> v = convert $ 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 v a, Vector w (v a), Vector w a) => v a -> w (v a) -> v a-{-# INLINE (<.#) #-}-v <.# m | any null m = empty- | otherwise = (v <.> convert (map head m)) `cons` (v <.# map tail m)---- | Matrix transpose------ >>> transpose $ fromList [fromList "ab", fromList "cd"]--- fromList [fromList "ac",fromList "bd"]-transpose :: (Vector v a, Vector w (v a), Vector w a) => w (v a) -> w (v a)-{-# INLINE transpose #-}-transpose m- | any null m = empty- | otherwise = convert (map head m) `cons` transpose (map tail m)
src/Learning/HMM.hs view
@@ -22,10 +22,15 @@ import Data.Random.Distribution.Categorical.Util () import Data.Random.RVar (RVar) import Data.Random.Sample (sample)-import qualified Data.Vector as V ((!), elemIndex, fromList, map, toList)+import qualified Data.Vector as V (+ elemIndex, fromList, map, toList, unsafeIndex+ ) import qualified Data.Vector.Generic as G (convert)-import qualified Data.Vector.Generic.Util.LinearAlgebra as G (transpose)-import qualified Data.Vector.Unboxed as U (fromList, toList)+import qualified Data.Vector.Unboxed as U (fromList)+import qualified Numeric.LinearAlgebra.Data as H (+ (!), fromList, fromLists, toList+ )+import qualified Numeric.LinearAlgebra.HMatrix as H (tr) import Learning.HMM.Internal -- | Parameter set of the hidden Markov model. Direct use of the@@ -33,13 +38,13 @@ data HMM s o = HMM { states :: [s] -- ^ Hidden states , outputs :: [o] -- ^ Outputs , initialStateDist :: Categorical Double s- -- ^ Categorical distribution of initial states+ -- ^ Categorical distribution of initial state , transitionDist :: s -> Categorical Double s- -- ^ Categorical distribution of next states- -- conditioned by the previous states+ -- ^ Categorical distribution of next state+ -- conditioned by the previous state , emissionDist :: s -> Categorical Double o- -- ^ Categorical distribution of outputs conditioned- -- by the hidden states+ -- ^ Categorical distribution of output conditioned+ -- by the hidden state } instance (Show s, Show o) => Show (HMM s o) where@@ -92,10 +97,10 @@ init ss os = fromHMM' ss os <$> init' (length ss) (length os) -- | @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.+-- 'emissionDist' is updated by re-estimations using the observed outputs+-- @xs@. The 'emissionDist' is set to be normalized histograms each of+-- which is calculated from segumentations of @xs@ based on the Viterbi+-- state path. withEmission :: (Eq s, Eq o) => HMM s o -> [o] -> HMM s o withEmission model xs = fromHMM' ss os $ withEmission' model' xs' where@@ -112,12 +117,13 @@ viterbi model xs = checkModelIn "viterbi" model `seq` checkDataIn "viterbi" model xs `seq`- first (V.toList . V.map (ss V.!) . G.convert) $ viterbi' model' xs'+ first toStates $ viterbi' model' xs' where- ss = V.fromList $ states model+ ss' = V.fromList $ states model os' = V.fromList $ outputs model model' = toHMM' model xs' = U.fromList $ fromJust $ mapM (`V.elemIndex` os') xs+ toStates = V.toList . V.map (V.unsafeIndex ss') . G.convert -- | @baumWelch model xs@ iteratively performs the Baum-Welch algorithm -- using the observed outputs @xs@, and returns a list of updated models@@ -188,29 +194,29 @@ where pi0 = initialStateDist' hmm' w = transitionDist' hmm'- phi = G.transpose $ emissionDistT' hmm'- pi0' = zip (U.toList pi0) ss- w' i = zip (U.toList $ w V.! i) ss- phi' i = zip (U.toList $ phi V.! i) os+ phi = H.tr $ emissionDistT' hmm'+ pi0' = zip (H.toList pi0) ss+ w' i = zip (H.toList $ w H.! i) ss+ phi' i = zip (H.toList $ phi H.! i) os -- | Convert 'HMM' to 'HMM''. The 'initialStateDist'', 'transitionDist'', -- and 'emissionDistT'' are normalized. toHMM' :: (Eq s, Eq o) => HMM s o -> HMM' toHMM' hmm = HMM' { nStates' = length ss , nOutputs' = length os- , initialStateDist' = U.fromList pi0'- , transitionDist' = V.fromList w'- , emissionDistT' = V.fromList phi'+ , initialStateDist' = pi0+ , transitionDist' = w+ , emissionDistT' = phi' } where- ss = states hmm- os = outputs hmm- pi0 = C.normalizeCategoricalPs $ initialStateDist hmm- w = C.normalizeCategoricalPs . transitionDist hmm- phi = C.normalizeCategoricalPs . emissionDist hmm- pi0' = [pdf pi0 s | s <- ss]- w' = [U.fromList [pdf (w s) s' | s' <- ss] | s <- ss]- phi' = [U.fromList [pdf (phi s) o | s <- ss] | o <- os]+ ss = states hmm+ os = outputs hmm+ pi0_ = C.normalizeCategoricalPs $ initialStateDist hmm+ w_ = C.normalizeCategoricalPs . transitionDist hmm+ phi_ = C.normalizeCategoricalPs . emissionDist hmm+ pi0 = H.fromList [pdf pi0_ s | s <- ss]+ w = H.fromLists [[pdf (w_ s) s' | s' <- ss] | s <- ss]+ phi' = H.fromLists [[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
@@ -18,21 +18,28 @@ import Data.Random.RVar (RVar) import Data.Random.Distribution.Simplex (stdSimplex) import qualified Data.Vector as V (- Vector, (!), filter, foldl1', freeze, fromList, generate, map , tail- , zip, zipWith, zipWith3+ Vector, filter, foldl1', map, unsafeFreeze, unsafeIndex, unsafeTail+ , zip, zipWith3 ) import qualified Data.Vector.Generic as G (convert) import qualified Data.Vector.Generic.Util as G (frequencies)-import Data.Vector.Generic.Util.LinearAlgebra (- (>+>), (>.>), (>/>), (#+#), (.>), (>/), (#.>), (<.#)+import qualified Data.Vector.Mutable as MV (+ unsafeNew, unsafeRead, unsafeWrite )-import qualified Data.Vector.Generic.Util.LinearAlgebra as G (transpose)-import qualified Data.Vector.Mutable as MV (new, read, write) import qualified Data.Vector.Unboxed as U (- Vector, (!), freeze, fromList, generate, length, map, maxIndex, maximum- , replicate, sum, tail, zip+ Vector, fromList, length, map, sum, unsafeFreeze, unsafeIndex+ , unsafeTail, zip )-import qualified Data.Vector.Unboxed.Mutable as MU (new, read, write)+import qualified Data.Vector.Unboxed.Mutable as MU (+ unsafeNew, unsafeRead, unsafeWrite+ )+import qualified Numeric.LinearAlgebra.Data as H (+ (!), Matrix, Vector, diag, fromColumns, fromList, fromLists+ , fromRows, konst, maxElement, maxIndex, toColumns, tr+ )+import qualified Numeric.LinearAlgebra.HMatrix as H (+ (<>), (#>), sumElements+ ) type LogLikelihood = Double @@ -43,43 +50,69 @@ -- in order to simplify the calculation. data HMM' = HMM' { nStates' :: Int -- ^ Number of states , nOutputs' :: Int -- ^ Number of outputs- , initialStateDist' :: U.Vector Double- , transitionDist' :: V.Vector (U.Vector Double)- , emissionDistT' :: V.Vector (U.Vector Double)+ , initialStateDist' :: H.Vector Double+ , transitionDist' :: H.Matrix Double+ , emissionDistT' :: H.Matrix Double } instance NFData HMM' where- rnf hmm' = rnf n `seq` rnf m `seq` rnf pi0 `seq` rnf w `seq` rnf phi'+ rnf hmm' = rnf k `seq` rnf l `seq` rnf pi0 `seq` rnf w `seq` rnf phi' where- n = nStates' hmm'- m = nOutputs' hmm'+ k = nStates' hmm'+ l = nOutputs' hmm' pi0 = initialStateDist' hmm' w = transitionDist' hmm' phi' = emissionDistT' hmm' init' :: Int -> Int -> RVar HMM'-init' n m = do- pi0 <- U.fromList <$> stdSimplex (n-1)- w <- V.fromList <$> replicateM n (U.fromList <$> stdSimplex (n-1))- phi <- V.fromList <$> replicateM n (U.fromList <$> stdSimplex (m-1))- return HMM' { nStates' = n- , nOutputs' = m+init' k l = do+ pi0 <- H.fromList <$> stdSimplex (k-1)+ w <- H.fromLists <$> replicateM k (stdSimplex (k-1))+ phi <- H.fromLists <$> replicateM k (stdSimplex (l-1))+ return HMM' { nStates' = k+ , nOutputs' = l , initialStateDist' = pi0 , transitionDist' = w- , emissionDistT' = G.transpose phi+ , emissionDistT' = H.tr phi } withEmission' :: HMM' -> U.Vector Int -> HMM'-withEmission' model xs = model { emissionDistT' = phi' }+withEmission' model xs = model' where- ss = V.generate (nStates' model) id- os = U.generate (nOutputs' model) id- phi' = let (path, _) = viterbi' model xs- freqs = G.frequencies $ U.zip path xs- hists = V.map (\s -> U.map (\o ->- M.findWithDefault 0 (s, o) freqs) os) ss- in V.map (\f -> f >/ U.sum f) hists+ n = U.length xs+ k = nStates' model+ l = nOutputs' model+ ss = [0..(k-1)]+ os = [0..(l-1)] + step m = fst $ baumWelch1' (m { emissionDistT' = H.tr phi }) n xs+ where+ phi :: H.Matrix Double+ phi = let zs = fst $ viterbi' m xs+ fs = G.frequencies $ U.zip zs xs+ hs = H.fromLists $ map (\s -> map (\o ->+ M.findWithDefault 0 (s, o) fs) os) ss+ -- hs' is needed to not yield NaN vectors+ hs' = hs + H.konst 1e-9 (k, l)+ ns = hs' H.#> H.konst 1 k+ in hs' / H.fromColumns (replicate l ns)++ ms = iterate step model+ ms' = tail ms+ ds = zipWith euclideanDistance ms ms'++ model' = fst $ head $ dropWhile ((> 1e-9) . snd) $ zip ms' ds++-- | Return the Euclidean distance between two models.+euclideanDistance :: HMM' -> HMM' -> Double+euclideanDistance model model' =+ sqrt $ (H.sumElements $ (w - w') ** 2) + (H.sumElements $ (phi - phi') ** 2)+ where+ w = transitionDist' model+ w' = transitionDist' model'+ phi = emissionDistT' model+ phi' = emissionDistT' model'+ viterbi' :: HMM' -> U.Vector Int -> (U.Vector Int, LogLikelihood) viterbi' model xs = (path, logL) where@@ -87,36 +120,39 @@ -- First, we calculate the value function and the state maximizing it -- for each time.- deltas :: V.Vector (U.Vector Double)+ deltas :: V.Vector (H.Vector Double) psis :: V.Vector (U.Vector Int) (deltas, psis) = runST $ do- ds <- MV.new n- ps <- MV.new n- MV.write ds 0 $ U.map log (phi' V.! (xs U.! 0)) >+> U.map log pi0- MV.write ps 0 $ U.replicate k 0+ ds <- MV.unsafeNew n+ ps <- MV.unsafeNew n+ let x0 = U.unsafeIndex xs 0+ MV.unsafeWrite ds 0 $ log (phi' H.! x0) + log pi0 forM_ [1..(n-1)] $ \t -> do- d <- MV.read ds (t-1)- let dws = V.map (\wj -> d >+> U.map log wj) w'- MV.write ds t $ U.map log (phi' V.! (xs U.! t)) >+> G.convert (V.map U.maximum dws)- MV.write ps t $ G.convert (V.map U.maxIndex dws)- ds' <- V.freeze ds- ps' <- V.freeze ps+ d <- MV.unsafeRead ds (t-1)+ let x = U.unsafeIndex xs t+ dws = map (\wj -> d + log wj) w'+ MV.unsafeWrite ds t $ log (phi' H.! x) + H.fromList (map H.maxElement dws)+ MV.unsafeWrite ps t $ U.fromList (map H.maxIndex dws)+ ds' <- V.unsafeFreeze ds+ ps' <- V.unsafeFreeze ps return (ds', ps') where- k = nStates' model pi0 = initialStateDist' model- w' = G.transpose $ transitionDist' model+ w' = H.toColumns $ transitionDist' model phi' = emissionDistT' model - -- The most likely path and corresponding log likelihood is as follows.+ deltaE = V.unsafeIndex deltas (n-1)++ -- The most likely path and corresponding log likelihood are as follows. path = runST $ do- ix <- MU.new n- MU.write ix (n-1) $ U.maxIndex (deltas V.! (n-1))- forM_ (reverse [0..(n-2)]) $ \t -> do- i <- MU.read ix (t+1)- MU.write ix t $ psis V.! (t+1) U.! i- U.freeze ix- logL = U.maximum $ deltas V.! (n-1)+ ix <- MU.unsafeNew n+ MU.unsafeWrite ix (n-1) $ H.maxIndex deltaE+ forM_ [n-l | l <- [1..(n-1)]] $ \t -> do+ i <- MU.unsafeRead ix t+ let psi = V.unsafeIndex psis t+ MU.unsafeWrite ix (t-1) $ U.unsafeIndex psi i+ U.unsafeFreeze ix+ logL = H.maxElement deltaE baumWelch' :: HMM' -> U.Vector Int -> [(HMM', LogLikelihood)] baumWelch' model xs = zip models (tail logLs)@@ -130,6 +166,9 @@ baumWelch1' :: HMM' -> Int -> U.Vector Int -> (HMM', LogLikelihood) baumWelch1' model n xs = force (model', logL) where+ k = nStates' model+ l = nOutputs' model+ -- First, we calculate the alpha, beta, and scaling values using the -- forward-backward algorithm. (alphas, cs) = forward' model n xs@@ -142,16 +181,14 @@ -- Using the gamma and xi values, we obtain the optimal initial state -- probability vector, transition probability matrix, and emission -- probability matrix.- pi0 = gammas V.! 0- w = let ds = V.foldl1' (#+#) xis -- denominators- ns = V.map U.sum ds -- numerators- in V.zipWith (>/) ds ns+ pi0 = V.unsafeIndex gammas 0+ w = let ds = V.foldl1' (+) xis -- denominators+ ns = ds H.#> H.konst 1 k -- numerators+ in H.diag (H.konst 1 k / ns) H.<> ds phi' = let gs' o = V.map snd $ V.filter ((== o) . fst) $ V.zip (G.convert xs) gammas- ds = V.foldl1' (>+>) . gs' -- denominators- ns = V.foldl1' (>+>) gammas -- numerators- in V.map (\o -> ds o >/> ns) os- where- os = V.generate (nOutputs' model) id+ ds = V.foldl1' (+) . gs' -- denominators+ ns = V.foldl1' (+) gammas -- numerators+ in H.fromRows $ map (\o -> ds o / ns) [0..(l-1)] -- We finally obtain the new model and the likelihood for the old model. model' = model { initialStateDist' = pi0@@ -161,56 +198,61 @@ logL = - (U.sum $ U.map log cs) -- | Return alphas and scaling variables.-forward' :: HMM' -> Int -> U.Vector Int -> (V.Vector (U.Vector Double), U.Vector Double)+forward' :: HMM' -> Int -> U.Vector Int -> (V.Vector (H.Vector Double), U.Vector Double) {-# INLINE forward' #-} forward' model n xs = runST $ do- as <- MV.new n- cs <- MU.new n- let a0 = (phi' V.! (xs U.! 0)) >.> pi0- c0 = 1 / U.sum a0- MV.write as 0 (c0 .> a0)- MU.write cs 0 c0+ as <- MV.unsafeNew n+ cs <- MU.unsafeNew n+ let x0 = U.unsafeIndex xs 0+ a0 = (phi' H.! x0) * pi0+ c0 = 1 / H.sumElements a0+ MV.unsafeWrite as 0 (H.konst c0 k * a0)+ MU.unsafeWrite cs 0 c0 forM_ [1..(n-1)] $ \t -> do- a <- MV.read as (t-1)- let a' = (phi' V.! (xs U.! t)) >.> (a <.# w)- c' = 1 / U.sum a'- MV.write as t (c' .> a')- MU.write cs t c'- as' <- V.freeze as- cs' <- U.freeze cs+ a <- MV.unsafeRead as (t-1)+ let x = U.unsafeIndex xs t+ a' = (phi' H.! x) * (w' H.#> a)+ c' = 1 / H.sumElements a'+ MV.unsafeWrite as t (H.konst c' k * a')+ MU.unsafeWrite cs t c'+ as' <- V.unsafeFreeze as+ cs' <- U.unsafeFreeze cs return (as', cs') where+ k = nStates' model pi0 = initialStateDist' model- w = transitionDist' model+ w' = H.tr $ transitionDist' model phi' = emissionDistT' model -- | Return betas using scaling variables.-backward' :: HMM' -> Int -> U.Vector Int -> U.Vector Double -> V.Vector (U.Vector Double)+backward' :: HMM' -> Int -> U.Vector Int -> U.Vector Double -> V.Vector (H.Vector Double) {-# INLINE backward' #-} backward' model n xs cs = runST $ do- bs <- MV.new n- let bE = U.replicate k 1- cE = cs U.! (n-1)- MV.write bs (n-1) $ cE .> bE- forM_ (reverse [0..(n-2)]) $ \t -> do- b <- MV.read bs (t+1)- let b' = w #.> ((phi' V.! (xs U.! (t+1))) >.> b)- c' = cs U.! t- MV.write bs t $ c' .> b'- V.freeze bs+ bs <- MV.unsafeNew n+ let bE = H.konst 1 k+ cE = U.unsafeIndex cs (n-1)+ MV.unsafeWrite bs (n-1) (H.konst cE k * bE)+ forM_ [n-l | l <- [1..(n-1)]] $ \t -> do+ b <- MV.unsafeRead bs t+ let x = U.unsafeIndex xs t+ b' = w H.#> ((phi' H.! x) * b)+ c' = U.unsafeIndex cs (t-1)+ MV.unsafeWrite bs (t-1) (H.konst c' k * b')+ V.unsafeFreeze bs where k = nStates' model w = transitionDist' model phi' = emissionDistT' model -- | Return the posterior distribution.-posterior' :: HMM' -> Int -> U.Vector Int -> V.Vector (U.Vector Double) -> V.Vector (U.Vector Double) -> U.Vector Double -> (V.Vector (U.Vector Double), V.Vector (V.Vector (U.Vector Double)))+posterior' :: HMM' -> Int -> U.Vector Int -> V.Vector (H.Vector Double) -> V.Vector (H.Vector Double) -> U.Vector Double -> (V.Vector (H.Vector Double), V.Vector (H.Matrix Double)) {-# INLINE posterior' #-} posterior' model _ xs alphas betas cs = (gammas, xis) where- gammas = V.zipWith3 (\a b c -> a >.> b >/ c) alphas betas (G.convert cs)- xis = V.zipWith3 (\a b x -> let w' = V.zipWith (.>) (G.convert a) w- in V.map ((phi' V.! x) >.> b >.>) w')- alphas (V.tail betas) (G.convert $ U.tail xs)+ gammas = V.zipWith3 (\a b c -> a * b / H.konst c k)+ alphas betas (G.convert cs)+ xis = V.zipWith3 (\a b x -> H.diag a H.<> w H.<> H.diag (b * (phi' H.! x)))+ alphas (V.unsafeTail betas) (G.convert $ U.unsafeTail xs)+ k = nStates' model w = transitionDist' model phi' = emissionDistT' model
tests/doctests.hs view
@@ -4,6 +4,5 @@ main :: IO () main = doctest [ "-isrc"- , "src/Data/Vector/Generic/Util/LinearAlgebra.hs" , "src/Learning/HMM.hs" ]