diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright (c) 2015, Henning Thielemann + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * 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. + + * Neither the name of Henning Thielemann nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND 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 COPYRIGHT +OWNER 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. diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,3 @@ +#! /usr/bin/env runhaskell +> import Distribution.Simple +> main = defaultMain diff --git a/hmm-hmatrix.cabal b/hmm-hmatrix.cabal new file mode 100644 --- /dev/null +++ b/hmm-hmatrix.cabal @@ -0,0 +1,79 @@ +Name: hmm-hmatrix +Version: 0.0 +Synopsis: Hidden Markov Models using HMatrix primitives +Description: + Hidden Markov Models implemented using HMatrix data types and operations. + + . + It implements: + . + * generation of samples of emission sequences, + . + * computation of the likelihood of an observed sequence of emissions, + . + * construction of most likely state sequence + that produces an observed sequence of emissions, + . + * supervised and unsupervised training of the model by Baum-Welch algorithm. + . + It supports any kind of emission distribution, + where discrete and multivariate Gaussian distributions + are implemented as examples. + . + For an introduction please refer to the examples: + . + * "Math.HiddenMarkovModel.Example.TrafficLight" + . + * "Math.HiddenMarkovModel.Example.SineWave" + . + * "Math.HiddenMarkovModel.Example.Circle" + . + An alternative package without foreign calls is @hmm@. +Homepage: http://code.haskell.org/~thielema/hmm-hmatrix +License: BSD3 +License-File: LICENSE +Author: Henning Thielemann +Maintainer: haskell@henning-thielemann.de +Category: Math +Build-Type: Simple +Cabal-Version: >=1.10 + +Source-Repository this + Tag: 0.0 + Type: darcs + Location: http://code.haskell.org/~thielema/hmm-hmatrix + +Source-Repository head + Type: darcs + Location: http://code.haskell.org/~thielema/hmm-hmatrix + +Library + Exposed-Modules: + Math.HiddenMarkovModel + Math.HiddenMarkovModel.Named + Math.HiddenMarkovModel.Distribution + Math.HiddenMarkovModel.Pattern + Math.HiddenMarkovModel.Example.TrafficLight + Math.HiddenMarkovModel.Example.SineWave + Math.HiddenMarkovModel.Example.Circle + Other-Modules: + Math.HiddenMarkovModel.Normalized + Math.HiddenMarkovModel.Private + Math.HiddenMarkovModel.Utility + Math.HiddenMarkovModel.CSV + Math.HiddenMarkovModel.Test + Build-Depends: + hmatrix >=0.15 && <0.16, + explicit-exception >=0.1.7 && <0.2, + lazy-csv >=0.5 && <0.6, + random >=1.0 && <1.1, + transformers >= 0.2 && <0.5, + non-empty >=0.2.1 && <0.3, + semigroups >=0.8.4.1 && <0.17, + containers >=0.4.2 && <0.6, + array >=0.4 && <0.6, + utility-ht >=0.0.10 && <0.1, + base >=4.5 && <4.8 + Hs-Source-Dirs: src + Default-Language: Haskell2010 + GHC-Options: -Wall diff --git a/src/Math/HiddenMarkovModel.hs b/src/Math/HiddenMarkovModel.hs new file mode 100644 --- /dev/null +++ b/src/Math/HiddenMarkovModel.hs @@ -0,0 +1,178 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +module Math.HiddenMarkovModel ( + T(..), Distr.State, state, + Discrete, DiscreteTrained, + Gaussian, GaussianTrained, + uniform, + generate, + Normalized.logLikelihood, + Normalized.reveal, + + Trained(..), + trainSupervised, + Normalized.trainUnsupervised, + mergeTrained, finishTraining, trainMany, + deviation, + + toCSV, + fromCSV, + ) where + +import qualified Math.HiddenMarkovModel.Distribution as Distr +import qualified Math.HiddenMarkovModel.Normalized as Normalized +import qualified Math.HiddenMarkovModel.CSV as HMMCSV +import Math.HiddenMarkovModel.Private + (T(..), Trained(..), mergeTrained, toCells, parseCSV) +import Math.HiddenMarkovModel.Distribution (State(State)) +import Math.HiddenMarkovModel.Utility + (randomItemProp, normalizeProb, attachOnes) + +import qualified Numeric.LinearAlgebra.Algorithms as Algo +import qualified Numeric.Container as NC +import qualified Data.Packed.Matrix as Matrix +import qualified Data.Packed.Vector as Vector +import Data.Packed.Matrix (Matrix) +import Data.Packed.Vector (Vector) + +import qualified Text.CSV.Lazy.String as CSV + +import qualified System.Random as Rnd + +import qualified Control.Monad.Exception.Synchronous as ME +import qualified Control.Monad.Trans.State as MS +import qualified Control.Monad.HT as Monad + +import qualified Data.NonEmpty as NonEmpty +import qualified Data.Array as Array +import Data.Foldable (Foldable) +import Data.Array (accumArray) + + + +state :: Int -> State +state = State + + +type DiscreteTrained prob symbol = Trained (Distr.DiscreteTrained prob symbol) prob +type Discrete prob symbol = T (Distr.Discrete prob symbol) prob + +type GaussianTrained a = Trained (Distr.GaussianTrained a) a +type Gaussian a = T (Distr.Gaussian a) a + + +{- | +Create a model with uniform probabilities +for initial vector and transition matrix +given a distribution for the emissions. +You can use this as a starting point for 'Normalized.trainUnsupervised'. +-} +uniform :: + (Distr.Info distr, Distr.Probability distr ~ prob) => + distr -> T distr prob +uniform distr = + let n = Distr.numberOfStates distr + c = recip $ fromIntegral n + in Cons { + initial = NC.constant c n, + transition = NC.konst c (n,n), + distribution = distr + } + + +generate :: + (Rnd.RandomGen g, Ord prob, Rnd.Random prob, + Distr.Generate distr, Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) => + T distr prob -> g -> [emission] +generate hmm = + MS.evalState $ + flip MS.evalStateT (initial hmm) $ + Monad.repeat $ MS.StateT $ \v0 -> do + s <- randomItemProp $ zip [0..] (Vector.toList v0) + x <- Distr.generate (distribution hmm) (State s) + return (x, takeColumn s $ transition hmm) + +takeColumn :: (Matrix.Element a) => Int -> Matrix a -> Vector a +takeColumn n = Matrix.flatten . Matrix.extractRows [n] . Matrix.trans + + + +{- | +Contribute a manually labeled emission sequence to a HMM training. +-} +trainSupervised :: + (Distr.Estimate tdistr, Distr.Distribution tdistr ~ distr, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) => + Int -> NonEmpty.T [] (State, emission) -> Trained tdistr prob +trainSupervised n xs = + let getState (State s, _x) = s + in Trained { + trainedInitial = NC.assoc n 0 [(getState (NonEmpty.head xs), 1)], + trainedTransition = + Matrix.trans $ NC.accum (NC.konst 0 (n,n)) (+) $ + attachOnes $ NonEmpty.mapAdjacent (,) $ fmap getState xs, + trainedDistribution = + Distr.accumulateEmissions $ map attachOnes $ Array.elems $ + accumArray (flip (:)) [] (State 0, State (n-1)) $ NonEmpty.flatten xs + } + +finishTraining :: + (Distr.Estimate tdistr, Distr.Distribution tdistr ~ distr, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) => + Trained tdistr prob -> T distr prob +finishTraining hmm = + Cons { + initial = normalizeProb $ trainedInitial hmm, + transition = + Matrix.fromColumns $ map normalizeProb $ + Matrix.toColumns $ trainedTransition hmm, + distribution = Distr.normalize $ trainedDistribution hmm + } + +trainMany :: + (Distr.Estimate tdistr, Distr.Distribution tdistr ~ distr, + Distr.Probability distr ~ prob, + Foldable f) => + (trainingData -> Trained tdistr prob) -> + NonEmpty.T f trainingData -> T distr prob +trainMany train = + finishTraining . NonEmpty.foldl1Map mergeTrained train + + + + + +{- | +Compute maximum deviation between initial and transition probabilities. +You can use this as abort criterion for unsupervised training. +We omit computation of differences between the emission probabilities. +This simplifies matters a lot and +should suffice for defining an abort criterion. +-} +deviation :: + (Algo.Field prob, Ord prob) => T distr prob -> T distr prob -> prob +deviation hmm0 hmm1 = + deviationVec (initial hmm0) (initial hmm1) + `max` + deviationVec (transition hmm0) (transition hmm1) + +deviationVec :: + (Ord a, NC.Container c a) => + c a -> c a -> a +deviationVec x y = + let d = NC.sub x y + in NC.maxElement d `max` negate (NC.minElement d) + + +toCSV :: + (Distr.CSV distr, Algo.Field prob, Show prob) => + T distr prob -> String +toCSV hmm = + CSV.ppCSVTable $ snd $ CSV.toCSVTable $ HMMCSV.padTable "" $ + toCells hmm + +fromCSV :: + (Distr.CSV distr, Algo.Field prob, Read prob) => + String -> ME.Exceptional String (T distr prob) +fromCSV = + MS.evalStateT parseCSV . map HMMCSV.fixShortRow . CSV.parseCSV diff --git a/src/Math/HiddenMarkovModel/CSV.hs b/src/Math/HiddenMarkovModel/CSV.hs new file mode 100644 --- /dev/null +++ b/src/Math/HiddenMarkovModel/CSV.hs @@ -0,0 +1,151 @@ +module Math.HiddenMarkovModel.CSV where + +import qualified Numeric.LinearAlgebra.Algorithms as Algo +import qualified Data.Packed.Matrix as Matrix +import qualified Data.Packed.Vector as Vector +import Data.Packed.Matrix (Matrix) +import Data.Packed.Vector (Vector) + +import qualified Text.CSV.Lazy.String as CSV +import Text.Read.HT (maybeRead) +import Text.Printf (printf) + +import qualified Control.Monad.Exception.Synchronous as ME +import qualified Control.Monad.Trans.Class as MT +import qualified Control.Monad.Trans.State as MS +import Control.Monad.Exception.Synchronous (Exceptional) +import Control.Monad (liftM2, replicateM, unless) + +import qualified Data.List.HT as ListHT + + +cellsFromVector :: + (Show a, Algo.Field a) => + Vector a -> [String] +cellsFromVector = map show . Vector.toList + +cellsFromMatrix :: + (Show a, Matrix.Element a) => + Matrix.Matrix a -> [[String]] +cellsFromMatrix = map (map show) . Matrix.toLists + +padTable :: a -> [[a]] -> [[a]] +padTable x xs = + let width = maximum (map length xs) + in map (ListHT.padRight x width) xs + + +type CSVParser = MS.StateT CSV.CSVResult (Exceptional String) + +assert :: Bool -> String -> CSVParser () +assert cond msg = + unless cond $ MT.lift $ ME.throw msg + +retrieveShortRow :: CSV.CSVError -> Maybe CSV.CSVRow +retrieveShortRow err = + case err of + CSV.IncorrectRow {CSV.csvFields = row} -> Just row + _ -> Nothing + +fixShortRow :: + Either [CSV.CSVError] CSV.CSVRow -> Either [CSV.CSVError] CSV.CSVRow +fixShortRow erow = + case erow of + Left errs -> + case ListHT.partitionMaybe retrieveShortRow errs of + ([row], []) -> Right row + _ -> Left errs + _ -> erow + +maybeGetRow :: CSVParser (Maybe CSV.CSVRow) +maybeGetRow = do + csv0 <- MS.get + case csv0 of + [] -> return Nothing + item : csv1 -> do + MS.put csv1 + case item of + Right row -> return (Just row) + Left errors -> + MT.lift $ ME.throw $ unlines $ map CSV.ppCSVError errors + +getRow :: CSVParser CSV.CSVRow +getRow = + MT.lift . ME.fromMaybe "unexpected end of file" =<< maybeGetRow + +checkEmptyRow :: CSV.CSVRow -> Exceptional String () +checkEmptyRow row = + case filter (not . null . CSV.csvFieldContent) row of + [] -> return () + cell:_ -> ME.throw $ printf "%d: expected empty row" (CSV.csvRowNum cell) + +skipEmptyRow :: CSVParser () +skipEmptyRow = MT.lift . checkEmptyRow =<< getRow + +manySepUntilEnd :: CSVParser a -> CSVParser [a] +manySepUntilEnd p = + let go = liftM2 (:) p $ do + mrow <- maybeGetRow + case mrow of + Nothing -> return [] + Just row -> do + MT.lift $ checkEmptyRow row + go + in go + +manyRowsUntilEnd :: (CSV.CSVRow -> CSVParser a) -> CSVParser [a] +manyRowsUntilEnd p = + let go = do + mrow <- maybeGetRow + case mrow of + Nothing -> return [] + Just row -> liftM2 (:) (p row) go + in go + +parseVectorCells :: + (Read a, Algo.Field a) => + CSVParser (Vector a) +parseVectorCells = + parseVectorFields =<< getRow + +parseVectorFields :: + (Read a, Algo.Field a) => + CSV.CSVRow -> CSVParser (Vector a) +parseVectorFields = + MT.lift . fmap Vector.fromList . mapM parseNumberCell . + ListHT.dropWhileRev (null . CSV.csvFieldContent) + +parseNonEmptyVectorCells :: + (Read a, Algo.Field a) => + CSVParser (Vector a) +parseNonEmptyVectorCells = do + v <- parseVectorCells + assert (Vector.dim v > 0) "no data for vector" + return v + +cellContent :: CSV.CSVField -> Exceptional String String +cellContent field = + case field of + CSV.CSVFieldError {} -> ME.throw $ CSV.ppCSVField field + CSV.CSVField { CSV.csvFieldContent = str } -> return str + +parseNumberCell :: (Read a) => CSV.CSVField -> Exceptional String a +parseNumberCell field = do + str <- cellContent field + ME.fromMaybe (printf "field content \"%s\" is not a number" str) $ + maybeRead str + +parseSquareMatrixCells :: + (Read a, Algo.Field a) => + Int -> CSVParser (Matrix a) +parseSquareMatrixCells n = do + rows <- replicateM n parseVectorCells + assert (not $ null rows) "no rows" + assert (all ((n==) . Vector.dim) rows) "inconsistent matrix dimensions" + return $ Matrix.fromRows rows + +parseStringList :: + CSV.CSVRow -> CSVParser [String] +parseStringList = + MT.lift . mapM cellContent . + ListHT.dropWhileRev (null . CSV.csvFieldContent) diff --git a/src/Math/HiddenMarkovModel/Distribution.hs b/src/Math/HiddenMarkovModel/Distribution.hs new file mode 100644 --- /dev/null +++ b/src/Math/HiddenMarkovModel/Distribution.hs @@ -0,0 +1,327 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +module Math.HiddenMarkovModel.Distribution ( + State(..), + Emission, Probability, Trained, + Info(..), Generate(..), EmissionProb(..), Estimate(..), + + Discrete(..), DiscreteTrained(..), + Gaussian(..), GaussianTrained(..), gaussian, + + CSV(..), HMMCSV.CSVParser, CSVSymbol(..), + ) where + +import qualified Math.HiddenMarkovModel.CSV as HMMCSV +import Math.HiddenMarkovModel.Utility (randomItemProp, normalizeProb) + +import qualified Numeric.LinearAlgebra.Algorithms as Algo +import qualified Numeric.Container as NC +import qualified Data.Packed.Matrix as Matrix +import qualified Data.Packed.Vector as Vector +import Numeric.Container ((<>)) +import Data.Packed.Matrix (Matrix) +import Data.Packed.Vector (Vector) + +import qualified System.Random as Rnd + +import qualified Text.CSV.Lazy.String as CSV +import Text.Read.HT (maybeRead) +import Text.Printf (printf) + +import qualified Control.Monad.Exception.Synchronous as ME +import qualified Control.Monad.Trans.Class as MT +import qualified Control.Monad.Trans.State as MS +import Control.Monad (liftM2) + +import qualified Data.NonEmpty as NonEmpty +import qualified Data.Foldable as Fold +import qualified Data.Map as Map +import qualified Data.Set as Set +import qualified Data.Array as Array +import qualified Data.List as List +import Data.Foldable (foldMap) +import Data.Tuple.HT (mapFst) +import Data.Array (Array, Ix, listArray, (!)) +import Data.Map (Map) +import Data.Maybe (listToMaybe) + + +newtype State = State Int + deriving (Eq, Ord, Show, Read, Ix) + +instance Enum State where + toEnum = State + fromEnum (State n) = n + + +type family Probability distr +type family Emission distr +type family Trained distr + + +class + (NC.Container Vector (Probability distr), NC.Product (Probability distr)) => + Info distr where + numberOfStates :: distr -> Int + +class + (NC.Container Vector (Probability distr), NC.Product (Probability distr)) => + Generate distr where + generate :: + (Rnd.RandomGen g, Probability distr ~ prob, Emission distr ~ emission) => + distr -> State -> MS.State g emission + +class + (NC.Container Vector (Probability distr), NC.Product (Probability distr)) => + EmissionProb distr where + emissionProb :: distr -> Emission distr -> Vector (Probability distr) + +class + (EmissionProb (Distribution tdistr), + Trained (Distribution tdistr) ~ tdistr) => + Estimate tdistr where + type Distribution tdistr + accumulateEmissions :: + (Distribution tdistr ~ distr, Probability distr ~ prob) => + [[(Emission distr, prob)]] -> tdistr + -- could as well be in Semigroup class + combine :: tdistr -> tdistr -> tdistr + normalize :: (Distribution tdistr ~ distr) => tdistr -> distr + + + +newtype Discrete prob symbol = Discrete (Map symbol (Vector prob)) + deriving (Show) + +newtype DiscreteTrained prob symbol = DiscreteTrained (Map symbol (Vector prob)) + deriving (Show) + +type instance Probability (Discrete prob symbol) = prob +type instance Emission (Discrete prob symbol) = symbol + +type instance Trained (Discrete prob symbol) = DiscreteTrained prob symbol + +instance + (NC.Container Vector prob, NC.Product prob, Ord symbol) => + Info (Discrete prob symbol) where + numberOfStates (Discrete m) = Vector.dim $ snd $ Map.findMin m + +instance + (NC.Container Vector prob, NC.Product prob, Ord symbol, + Ord prob, Rnd.Random prob) => + Generate (Discrete prob symbol) where + generate (Discrete m) (State state) = + randomItemProp $ Map.toAscList $ fmap (flip NC.atIndex state) m + +instance + (NC.Container Vector prob, NC.Product prob, Ord symbol) => + EmissionProb (Discrete prob symbol) where + emissionProb (Discrete m) = + mapLookup "emitDiscrete: unknown emission symbol" m + +instance + (NC.Container Vector prob, NC.Product prob, Ord symbol) => + Estimate (DiscreteTrained prob symbol) where + type Distribution (DiscreteTrained prob symbol) = Discrete prob symbol + accumulateEmissions grouped = + let set = Set.toAscList $ foldMap (Set.fromList . map fst) grouped + emi = Map.fromAscList $ zip set [0..] + in DiscreteTrained $ Map.fromAscList $ zip set $ + transposeVectorList $ + map + (NC.accum (NC.konst 0 (length set)) (+) . + map (mapFst + (mapLookup "estimateDiscrete: unknown emission symbol" emi))) + grouped + combine (DiscreteTrained distr0) (DiscreteTrained distr1) = + DiscreteTrained $ Map.unionWith NC.add distr0 distr1 + normalize (DiscreteTrained distr) = + Discrete $ Map.fromAscList $ zip (Map.keys distr) $ + transposeVectorList $ map normalizeProb $ + transposeVectorList $ Map.elems distr + +transposeVectorList :: (Matrix.Element a) => [Vector a] -> [Vector a] +transposeVectorList = Matrix.toRows . Matrix.fromColumns + +mapLookup :: (Ord k) => String -> Map.Map k a -> k -> a +mapLookup msg dict x = + Map.findWithDefault (error msg) x dict + + +newtype Gaussian a = Gaussian (Array State (Vector a, Matrix a, a)) + deriving (Show) + +newtype GaussianTrained a = GaussianTrained (Map State (Vector a, Matrix a, a)) + deriving (Show) + +type instance Probability (Gaussian a) = a +type instance Emission (Gaussian a) = Vector a + +type instance Trained (Gaussian a) = GaussianTrained a + +instance (Algo.Field a) => Info (Gaussian a) where + numberOfStates (Gaussian params) = Array.rangeSize $ Array.bounds params + +instance (Algo.Field a, Ord a, Rnd.Random a) => Generate (Gaussian a) where + generate (Gaussian allParams) state = do + let (center, covarianceChol, _c) = allParams ! state + seed <- MS.state Rnd.random + return $ + NC.add center $ + NC.cmap realToFrac + (NC.randomVector seed NC.Gaussian (Vector.dim center)) + <> covarianceChol + +instance (Algo.Field a) => EmissionProb (Gaussian a) where + emissionProb (Gaussian allParams) = + let cholSolve m x = + Matrix.flatten $ Algo.cholSolve m $ Matrix.asColumn x + prob x (center, covarianceChol, c) = + let x0 = NC.sub x center + in c * exp ((-1/2) * NC.dot x0 (cholSolve covarianceChol x0)) + in \x -> Vector.fromList $ map (prob x) $ Array.elems allParams + +instance (Algo.Field a) => Estimate (GaussianTrained a) where + type Distribution (GaussianTrained a) = Gaussian a + accumulateEmissions = + let params xs = + let center = + NonEmpty.foldl1Map NC.add (uncurry $ flip NC.scale) xs + covariance = + NonEmpty.foldl1Map NC.add (\(x,c) -> NC.scale c $ NC.outer x x) xs + in (center, covariance, Fold.sum $ fmap snd xs) + in GaussianTrained . fmap params . Map.mapMaybe NonEmpty.fetch . + Map.fromList . zip [State 0 ..] + combine (GaussianTrained distr0) (GaussianTrained distr1) = + let comb (center0, covariance0, weight0) + (center1, covariance1, weight1) = + (NC.add center0 center1, + NC.add covariance0 covariance1, + weight0 + weight1) + in GaussianTrained $ Map.unionWith comb distr0 distr1 + {- + Sum_i (xi-mi) * (xi-mi)^T + = Sum_i xi*xi^T + Sum_i mi*mi^T - Sum_i xi*mi^T - Sum_i mi*xi^T + = Sum_i xi*xi^T - Sum_i mi*mi^T + = Sum_i xi*xi^T - n * mi*mi^T + -} + normalize (GaussianTrained distr) = + let params (centerSum, covarianceSum, weight) = + let c = recip weight + center = NC.scale c centerSum + in (center, + NC.sub (NC.scale c covarianceSum) (NC.outer center center)) + in Gaussian $ + Array.array (fst $ Map.findMin distr, fst $ Map.findMax distr) $ + Map.toList $ fmap (gaussianParameters . params) distr + +gaussian :: + (Algo.Field prob) => + [(Vector prob, Matrix prob)] -> Gaussian prob +gaussian = + consGaussian . map gaussianParameters + +gaussianParameters :: + (Algo.Field prob) => + (Vector prob, Matrix prob) -> (Vector prob, Matrix prob, prob) +gaussianParameters (center, covariance) = + gaussianFromCholesky center $ Algo.chol covariance + +consGaussian :: [(Vector a, Matrix a, a)] -> Gaussian a +consGaussian xs = + Gaussian $ listArray (State 0, State $ length xs - 1) xs + +gaussianFromCholesky :: + (Algo.Field prob) => + Vector prob -> Matrix prob -> (Vector prob, Matrix prob, prob) +gaussianFromCholesky center covarianceChol = + let covarianceSqrtDet = NC.prodElements $ Matrix.takeDiag covarianceChol + in (center, covarianceChol, + recip (sqrt (2*pi) ^ Vector.dim center * covarianceSqrtDet)) + + + +class CSV distr where + toCells :: distr -> [[String]] + parseCells :: Int -> HMMCSV.CSVParser distr + +class (Ord symbol) => CSVSymbol symbol where + cellFromSymbol :: symbol -> String + symbolFromCell :: String -> Maybe symbol + +instance CSVSymbol Char where + cellFromSymbol = (:[]) + symbolFromCell = listToMaybe + +instance CSVSymbol Int where + cellFromSymbol = show + symbolFromCell = maybeRead + + +instance + (Algo.Field prob, Show prob, Read prob, CSVSymbol symbol) => + CSV (Discrete prob symbol) where + toCells (Discrete m) = + map + (\(symbol, probs) -> + cellFromSymbol symbol : HMMCSV.cellsFromVector probs) $ + Map.toAscList m + parseCells n = + fmap (Discrete . Map.fromList) $ + HMMCSV.manyRowsUntilEnd $ parseSymbolProb n + +parseSymbolProb :: + (Algo.Field prob, Read prob, CSVSymbol symbol) => + Int -> CSV.CSVRow -> HMMCSV.CSVParser (symbol, Vector prob) +parseSymbolProb n row = + case row of + [] -> MT.lift $ ME.throw "missing symbol" + c:cs -> + liftM2 (,) + (let str = CSV.csvFieldContent c + in MT.lift $ ME.fromMaybe (printf "unknown symbol %s" str) $ + symbolFromCell str) + (do v <- HMMCSV.parseVectorFields cs + HMMCSV.assert (n == Vector.dim v) + (printf "number of states (%d) and size of probability vector (%d) mismatch" + n (Vector.dim v)) + return v) + + +instance (Algo.Field a, Eq a, Show a, Read a) => CSV (Gaussian a) where + toCells (Gaussian params) = + List.intercalate [[]] $ + map + (\(center, covarianceChol, _) -> + HMMCSV.cellsFromVector center : + HMMCSV.cellsFromMatrix covarianceChol) $ + Array.elems params + parseCells n = do + gs <- HMMCSV.manySepUntilEnd parseSingleGaussian + HMMCSV.assert (length gs == n) $ + printf "number of states (%d) and number of Gaussians (%d) mismatch" + n (length gs) + return $ consGaussian gs + +parseSingleGaussian :: + (Algo.Field prob, Eq prob, Read prob) => + HMMCSV.CSVParser (Vector prob, Matrix prob, prob) +parseSingleGaussian = do + center <- HMMCSV.parseNonEmptyVectorCells + covarianceChol <- HMMCSV.parseSquareMatrixCells $ Vector.dim center + HMMCSV.assert (isUpperTriang covarianceChol) $ + "matrices must be upper triangular" + return $ gaussianFromCholesky center covarianceChol + + +{- +Maybe this test is too strict. +It would also be ok, and certainly more intuitive +to use an orthogonal but not normalized matrix. +We could get such a matrix from the eigensystem. +-} +isUpperTriang :: (Algo.Field a, Eq a) => Matrix a -> Bool +isUpperTriang m = + let cleared = Matrix.mapMatrixWithIndex (\(i,j) x -> if i>j then x else 0) m + in NC.minElement cleared == 0 && + NC.maxElement cleared == 0 diff --git a/src/Math/HiddenMarkovModel/Example/Circle.hs b/src/Math/HiddenMarkovModel/Example/Circle.hs new file mode 100644 --- /dev/null +++ b/src/Math/HiddenMarkovModel/Example/Circle.hs @@ -0,0 +1,93 @@ +{- | +Example of an HMM with continuous emissions with two-dimensional observations. +We train a model to accept a parametric curve of a circle with a certain speed. +This is like "Math.HiddenMarkovModel.Example.SineWave" but in two dimensions. + +The four hidden states correspond to the four quadrants. +-} +module Math.HiddenMarkovModel.Example.Circle +{-# WARNING "do not import that module, it is only intended for demonstration" #-} + where + +import qualified Math.HiddenMarkovModel as HMM +import qualified Math.HiddenMarkovModel.Distribution as Distr + +import qualified Data.Packed.Matrix as Matrix +import qualified Data.Packed.Vector as Vector +import Data.Packed.Vector (Vector) + +import qualified System.Random as Rnd + +import qualified Control.Monad.Trans.State as MS +import Control.Monad (liftM2, replicateM) + +import qualified Data.NonEmpty.Class as NonEmptyC +import qualified Data.NonEmpty as NonEmpty +import Data.Function.HT (nest) +import Data.NonEmpty ((!:)) + + + +hmm :: HMM.Gaussian Double +hmm = + HMM.Cons { + HMM.initial = Vector.fromList [1/4, 1/4, 1/4, 1/4], + HMM.transition = + Matrix.fromLists $ + [0.9, 0.0, 0.0, 0.1] : + [0.1, 0.9, 0.0, 0.0] : + [0.0, 0.1, 0.9, 0.0] : + [0.0, 0.0, 0.1, 0.9] : + [], + HMM.distribution = + let cov0 = Matrix.fromLists [[0.10, -0.09],[-0.09, 0.10]] + cov1 = Matrix.fromLists [[0.10, 0.09],[ 0.09, 0.10]] + in Distr.gaussian $ + (Vector.fromList [ 0.5, 0.5], cov0) : + (Vector.fromList [-0.5, 0.5], cov1) : + (Vector.fromList [-0.5, -0.5], cov0) : + (Vector.fromList [ 0.5, -0.5], cov1) : + [] + } + +circleLabeled :: NonEmpty.T [] (HMM.State, Vector Double) +circleLabeled = + NonEmpty.mapTail (take 200) $ + fmap + (\x -> + (HMM.state $ mod (floor (x*2/pi)) 4, + Vector.fromList [cos x, sin x])) $ + NonEmptyC.iterate (0.1+) 0 + +circle :: NonEmpty.T [] (Vector Double) +circle = fmap snd circleLabeled + +revealed :: NonEmpty.T [] HMM.State +revealed = HMM.reveal hmm circle + +{- | +Sample multivariate normal distribution and reconstruct it from the samples. +You should obtain the same parameters. +-} +reconstructDistribution :: HMM.Gaussian Double +reconstructDistribution = + let s0 = HMM.state 0 + gen = Distr.generate (HMM.distribution hmm) s0 + in HMM.finishTraining $ HMM.trainSupervised 1 $ fmap ((,) s0) $ + flip MS.evalState (Rnd.mkStdGen 23) $ + liftM2 (!:) gen $ replicateM 1000 gen + + +hmmTrainedSupervised :: HMM.Gaussian Double +hmmTrainedSupervised = + HMM.finishTraining $ HMM.trainSupervised 4 circleLabeled + +hmmTrainedUnsupervised :: HMM.Gaussian Double +hmmTrainedUnsupervised = + HMM.finishTraining $ HMM.trainUnsupervised hmm circle + +hmmIterativelyTrained :: HMM.Gaussian Double +hmmIterativelyTrained = + nest 100 + (HMM.finishTraining . flip HMM.trainUnsupervised circle) + hmm diff --git a/src/Math/HiddenMarkovModel/Example/SineWave.hs b/src/Math/HiddenMarkovModel/Example/SineWave.hs new file mode 100644 --- /dev/null +++ b/src/Math/HiddenMarkovModel/Example/SineWave.hs @@ -0,0 +1,76 @@ +{- | +Example of an HMM with continuous emissions. +We train a model to accept sine waves of a certain frequency. + +There are four hidden states: +0 - rising, +1 - high, +2 - falling, +3 - low. +-} +module Math.HiddenMarkovModel.Example.SineWave +{-# WARNING "do not import that module, it is only intended for demonstration" #-} + where + +import qualified Math.HiddenMarkovModel as HMM +import qualified Math.HiddenMarkovModel.Distribution as Distr + +import qualified Numeric.Container as NC +import qualified Data.Packed.Matrix as Matrix +import qualified Data.Packed.Vector as Vector + +import qualified Data.NonEmpty.Class as NonEmptyC +import qualified Data.NonEmpty as NonEmpty +import Data.Function.HT (nest) +import Data.Tuple.HT (mapSnd) + + + +hmm :: HMM.Gaussian Double +hmm = + HMM.Cons { + HMM.initial = Vector.fromList [1/4, 1/4, 1/4, 1/4], + HMM.transition = + Matrix.fromLists $ + [0.9, 0.0, 0.0, 0.1] : + [0.1, 0.9, 0.0, 0.0] : + [0.0, 0.1, 0.9, 0.0] : + [0.0, 0.0, 0.1, 0.9] : + [], + HMM.distribution = + Distr.gaussian $ + (Vector.fromList [ 0], Matrix.fromLists [[1]]) : + (Vector.fromList [ 1], Matrix.fromLists [[1]]) : + (Vector.fromList [ 0], Matrix.fromLists [[1]]) : + (Vector.fromList [-1], Matrix.fromLists [[1]]) : + [] + } + +sineWaveLabeled :: NonEmpty.T [] (HMM.State, Double) +sineWaveLabeled = + NonEmpty.mapTail (take 200) $ + fmap (\x -> (HMM.state $ mod (floor (x*2/pi+0.5)) 4, sin x)) $ + NonEmptyC.iterate (0.1+) 0 + +sineWave :: NonEmpty.T [] Double +sineWave = fmap snd sineWaveLabeled + +revealed :: NonEmpty.T [] HMM.State +revealed = HMM.reveal hmmTrainedSupervised $ fmap NC.scalar sineWave + +hmmTrainedSupervised :: HMM.Gaussian Double +hmmTrainedSupervised = + HMM.finishTraining $ HMM.trainSupervised 4 $ + fmap (mapSnd NC.scalar) sineWaveLabeled + +hmmTrainedUnsupervised :: HMM.Gaussian Double +hmmTrainedUnsupervised = + HMM.finishTraining $ HMM.trainUnsupervised hmm $ fmap NC.scalar sineWave + +hmmIterativelyTrained :: HMM.Gaussian Double +hmmIterativelyTrained = + nest 100 + (\model -> + HMM.finishTraining $ HMM.trainUnsupervised model $ + fmap NC.scalar sineWave) + hmm diff --git a/src/Math/HiddenMarkovModel/Example/TrafficLight.hs b/src/Math/HiddenMarkovModel/Example/TrafficLight.hs new file mode 100644 --- /dev/null +++ b/src/Math/HiddenMarkovModel/Example/TrafficLight.hs @@ -0,0 +1,172 @@ +{- | +This is an example of an HMM with discrete emissions. +We model a traffic light consisting of the colors red, yellow, green, +where only one lamp can be switched on at every point in time. +This way, when it is yellow you cannot tell immediately +whether it will switch to green or red. +We can only infer this from the light seen before. + +There are four hidden states: +0 emits red, 1 emits yellow between red and green, +2 emits green, 3 emits yellow between green and red. + +We quantise time in time steps. +The transition matrix of the model 'hmm' encodes +the expected duration of every state counted in time steps +and what states follow after each other. +E.g. transition probability of 0.8 of a state to itself means +that the expected duration of the state is 5 time steps (1/(1-0.8)). +However, it is a geometric distribution, +that is, shorter durations are always more probable. + +The distribution of 'hmm' encodes which lights a state activates. +In our case everything is deterministic: +Every state can switch exactly one light on. + +Given a sequence of observed lights +the function 'HMM.reveal' tells us the most likely sequence of states. +We test this with the light sequences in 'stateSequences' +where we already know the hidden states +as they are stored in 'labeledSequences'. +'verifyRevelation' compares the computed state sequence with the given one. + +We also try some trainings in 'hmmTrainedSupervised' et.al. +-} +module Math.HiddenMarkovModel.Example.TrafficLight +{-# WARNING "do not import that module, it is only intended for demonstration" #-} + where + +import qualified Math.HiddenMarkovModel as HMM +import qualified Math.HiddenMarkovModel.Distribution as Distr + +import qualified Data.Packed.Matrix as Matrix +import qualified Data.Packed.Vector as Vector + +import Text.Read.HT (maybeRead) + +import Control.Monad (liftM2) + +import qualified Data.Map as Map +import qualified Data.NonEmpty as NonEmpty +import qualified Data.List.HT as ListHT +import Data.NonEmpty ((!:)) + + + +data Color = Red | Yellow | Green + deriving (Eq, Ord, Enum, Show, Read) + +{- | +Using 'show' and 'read' is not always a good choice +since they must format and parse Haskell expressions +which is not of much use to the outside world. +-} +instance Distr.CSVSymbol Color where + cellFromSymbol = show + symbolFromCell = maybeRead + + +hmm :: HMM.Discrete Double Color +hmm = + HMM.Cons { + HMM.initial = Vector.fromList [1/3, 1/6, 1/3, 1/6], + HMM.transition = + Matrix.fromLists $ + [0.8, 0.0, 0.0, 0.2] : + [0.2, 0.8, 0.0, 0.0] : + [0.0, 0.2, 0.8, 0.0] : + [0.0, 0.0, 0.2, 0.8] : + [], + HMM.distribution = + Distr.Discrete $ Map.fromList $ + (Red, Vector.fromList [1,0,0,0]) : + (Yellow, Vector.fromList [0,1,0,1]) : + (Green, Vector.fromList [0,0,1,0]) : + [] + } + +hmmDisturbed :: HMM.Discrete Double Color +hmmDisturbed = + HMM.Cons { + HMM.initial = Vector.fromList [1/4, 1/4, 1/4, 1/4], + HMM.transition = + Matrix.fromLists $ + [0.3, 0.2, 0.2, 0.3] : + [0.3, 0.3, 0.2, 0.2] : + [0.2, 0.3, 0.3, 0.2] : + [0.2, 0.2, 0.3, 0.3] : + [], + HMM.distribution = + Distr.Discrete $ Map.fromList $ + (Red, Vector.fromList [0.6, 0.2, 0.2, 0.2]) : + (Yellow, Vector.fromList [0.2, 0.6, 0.2, 0.6]) : + (Green, Vector.fromList [0.2, 0.2, 0.6, 0.2]) : + [] + } + + +red, yellowRG, green, yellowGR :: (HMM.State, Color) +red = (HMM.state 0, Red) +yellowRG = (HMM.state 1, Yellow) +green = (HMM.state 2, Green) +yellowGR = (HMM.state 3, Yellow) + +labeledSequences :: NonEmpty.T [] (NonEmpty.T [] (HMM.State, Color)) +labeledSequences = + (red !: red : red : red : + yellowRG : yellowRG : + green : green : green : green : green : + yellowGR : + red : red : red : + []) !: + (green !: green : green : + yellowGR : + red : red : red : red : + yellowRG : + green : green : green : green : green : + yellowGR : yellowGR : + []) : + [] + +{- | +Construct a Hidden Markov model by watching a set +of manually created sequences of emissions and according states. +-} +hmmTrainedSupervised :: HMM.Discrete Double Color +hmmTrainedSupervised = + HMM.trainMany (HMM.trainSupervised 4) labeledSequences + + +stateSequences :: NonEmpty.T [] (NonEmpty.T [] Color) +stateSequences = fmap (fmap snd) labeledSequences + +{- | +Construct a Hidden Markov model starting from a known model +and a set of sequences that contain only the emissions, but no states. +-} +hmmTrainedUnsupervised :: HMM.Discrete Double Color +hmmTrainedUnsupervised = + HMM.trainMany (HMM.trainUnsupervised hmm) stateSequences + +{- | +Repeat unsupervised training until convergence. +-} +hmmIterativelyTrained :: HMM.Discrete Double Color +hmmIterativelyTrained = + snd $ head $ dropWhile fst $ + ListHT.mapAdjacent (\hmm0 hmm1 -> (HMM.deviation hmm0 hmm1 > 1e-5, hmm1)) $ + iterate + (flip HMM.trainMany stateSequences . HMM.trainUnsupervised) + hmmDisturbed + + +verifyRevelation :: + HMM.Discrete Double Color -> NonEmpty.T [] (HMM.State, Color) -> Bool +verifyRevelation model xs = + fmap fst xs == HMM.reveal model (fmap snd xs) + +verifyRevelations :: [Bool] +verifyRevelations = + liftM2 verifyRevelation + [hmm, hmmDisturbed, hmmTrainedSupervised, hmmTrainedUnsupervised] + (NonEmpty.flatten labeledSequences) diff --git a/src/Math/HiddenMarkovModel/Named.hs b/src/Math/HiddenMarkovModel/Named.hs new file mode 100644 --- /dev/null +++ b/src/Math/HiddenMarkovModel/Named.hs @@ -0,0 +1,92 @@ +module Math.HiddenMarkovModel.Named ( + T(..), + Discrete, + Gaussian, + fromModelAndNames, + toCSV, + fromCSV, + ) where + +import qualified Math.HiddenMarkovModel.Distribution as Distr +import qualified Math.HiddenMarkovModel.Private as HMM +import qualified Math.HiddenMarkovModel.CSV as HMMCSV +import Math.HiddenMarkovModel.Distribution (State(..)) +import Math.HiddenMarkovModel.Utility (attachOnes) + +import qualified Numeric.LinearAlgebra.Algorithms as Algo +import qualified Data.Packed.Vector as Vector + +import qualified Text.CSV.Lazy.String as CSV +import Text.Printf (printf) + +import qualified Control.Monad.Exception.Synchronous as ME +import qualified Control.Monad.Trans.State as MS + +import qualified Data.Map as Map +import qualified Data.List as List +import Data.Tuple.HT (swap) +import Data.Map (Map) + + +{- | +A Hidden Markov Model with names for each state. + +Although 'nameFromStateMap' and 'stateFromNameMap' are exported +you must be careful to keep them consistent when you alter them. +-} +data T distr prob = + Cons { + model :: HMM.T distr prob, + nameFromStateMap :: Map State String, + stateFromNameMap :: Map String State + } + deriving (Show, Read) + +type Discrete prob symbol = T (Distr.Discrete prob symbol) prob +type Gaussian a = T (Distr.Gaussian a) a + + +fromModelAndNames :: HMM.T distr prob -> [String] -> T distr prob +fromModelAndNames md names = + let m = Map.fromList $ zip [State 0 ..] names + in Cons { + model = md, + nameFromStateMap = m, + stateFromNameMap = inverseMap m + } + +inverseMap :: Map State String -> Map String State +inverseMap = + Map.fromListWith (error "duplicate label") . + map swap . Map.toList + + +toCSV :: + (Distr.CSV distr, Algo.Field prob, Show prob) => + T distr prob -> String +toCSV hmm = + CSV.ppCSVTable $ snd $ CSV.toCSVTable $ HMMCSV.padTable "" $ + Map.elems (nameFromStateMap hmm) : HMM.toCells (model hmm) + +fromCSV :: + (Distr.CSV distr, Algo.Field prob, Read prob) => + String -> ME.Exceptional String (T distr prob) +fromCSV = + MS.evalStateT parseCSV . map HMMCSV.fixShortRow . CSV.parseCSV + +parseCSV :: + (Distr.CSV distr, Algo.Field prob, Read prob) => + HMMCSV.CSVParser (T distr prob) +parseCSV = do + names <- HMMCSV.parseStringList =<< HMMCSV.getRow + let duplicateNames = + Map.keys $ Map.filter (> (1::Int)) $ + Map.fromListWith (+) $ attachOnes names + in HMMCSV.assert (null duplicateNames) $ + "duplicate names: " ++ List.intercalate ", " duplicateNames + md <- HMM.parseCSV + let n = length names + m = Vector.dim (HMM.initial md) + in HMMCSV.assert (n == m) $ + printf "got %d state names for %d state" n m + return $ fromModelAndNames md names diff --git a/src/Math/HiddenMarkovModel/Normalized.hs b/src/Math/HiddenMarkovModel/Normalized.hs new file mode 100644 --- /dev/null +++ b/src/Math/HiddenMarkovModel/Normalized.hs @@ -0,0 +1,168 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{- | +Counterparts to functions in "Math.HiddenMarkovModel.Private" +that normalize interim results. +We need to do this in order to prevent +to round very small probabilities to zero. +-} +module Math.HiddenMarkovModel.Normalized where + +import qualified Math.HiddenMarkovModel.Distribution as Distr +import Math.HiddenMarkovModel.Private + (T(..), Trained(..), emission, matrixMaxMul, sumTransitions) +import Math.HiddenMarkovModel.Distribution (State(State)) +import Math.HiddenMarkovModel.Utility (normalizeFactor, normalizeProb) + +import qualified Numeric.Container as NC +import qualified Data.Packed.Development as Dev +import qualified Data.Packed.Vector as Vector +import Numeric.Container ((<>)) +import Data.Packed.Matrix (Matrix) +import Data.Packed.Vector (Vector) + +import qualified Control.Functor.HT as Functor + +import qualified Data.NonEmpty.Class as NonEmptyC +import qualified Data.NonEmpty as NonEmpty +import qualified Data.Foldable as Fold +import qualified Data.List as List +import Data.Traversable (Traversable, mapAccumL) +import Data.Tuple.HT (mapFst, mapSnd, swap) + + +{- | +Logarithm of the likelihood to observe the given sequence. +We return the logarithm because the likelihood can be so small +that it may be rounded to zero in the choosen number type. +-} +logLikelihood :: + (Distr.EmissionProb distr, Floating prob, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission, + Traversable f) => + T distr prob -> NonEmpty.T f emission -> prob +logLikelihood hmm = Fold.sum . fmap (log . fst) . alpha hmm + +alpha :: + (Distr.EmissionProb distr, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission, + Traversable f) => + T distr prob -> + NonEmpty.T f emission -> NonEmpty.T f (prob, Vector prob) +alpha hmm (NonEmpty.Cons x xs) = + let normMulEmiss y = normalizeFactor . NC.mul (emission hmm y) + in NonEmpty.scanl + (\(_,alphai) xi -> normMulEmiss xi (transition hmm <> alphai)) + (normMulEmiss x (initial hmm)) + xs + +beta :: + (Distr.EmissionProb distr, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission, + Traversable f, NonEmptyC.Zip f, NonEmptyC.Reverse f) => + T distr prob -> + f (prob, emission) -> NonEmpty.T f (Vector prob) +beta hmm = + nonEmptyScanr + (\(ci,xi) betai -> + NC.scale (recip ci) $ NC.mul (emission hmm xi) betai <> transition hmm) + (NC.constant 1 (NC.dim $ initial hmm)) + +alphaBeta :: + (Distr.EmissionProb distr, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission, + Traversable f, NonEmptyC.Zip f, NonEmptyC.Reverse f) => + T distr prob -> + NonEmpty.T f emission -> + (NonEmpty.T f (prob, Vector prob), NonEmpty.T f (Vector prob)) +alphaBeta hmm xs = + let calphas = alpha hmm xs + in (calphas, + beta hmm $ NonEmpty.tail $ NonEmptyC.zip (fmap fst calphas) xs) + + +xiFromAlphaBeta :: + (Distr.EmissionProb distr, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission, + Traversable f, NonEmptyC.Zip f) => + T distr prob -> + NonEmpty.T f emission -> + NonEmpty.T f (prob, Vector prob) -> + NonEmpty.T f (Vector prob) -> + f (Matrix prob) +xiFromAlphaBeta hmm xs calphas betas = + let (cs,alphas) = Functor.unzip calphas + in NonEmptyC.zipWith4 + (\x alpha0 c1 beta1 -> + NC.scale (recip c1) $ + NC.mul + (NC.outer (NC.mul (emission hmm x) beta1) alpha0) + (transition hmm)) + (NonEmpty.tail xs) + (NonEmpty.init alphas) + (NonEmpty.tail cs) + (NonEmpty.tail betas) + +zetaFromAlphaBeta :: + (NC.Container Vector prob, NonEmptyC.Zip f) => + NonEmpty.T f (prob, Vector prob) -> + NonEmpty.T f (Vector prob) -> + NonEmpty.T f (Vector prob) +zetaFromAlphaBeta calphas betas = + NonEmptyC.zipWith (NC.mul . snd) calphas betas + + +{- | +Reveal the state sequence +that led most likely to the observed sequence of emissions. +It is found using the Viterbi algorithm. +-} +reveal :: + (Distr.EmissionProb distr, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission, + Traversable f, NonEmptyC.Reverse f) => + T distr prob -> NonEmpty.T f emission -> NonEmpty.T f State +reveal hmm (NonEmpty.Cons x xs) = + fmap State $ + uncurry (nonEmptyScanr Dev.at') $ + mapFst NC.maxIndex $ + mapAccumL + (\alphai xi -> + swap $ mapSnd (NC.mul (emission hmm xi)) $ + matrixMaxMul (transition hmm) $ normalizeProb alphai) + (NC.mul (emission hmm x) (initial hmm)) xs + + +{- | +Variant of NonEmpty.scanr with less stack consumption. +-} +nonEmptyScanr :: + (Traversable f, NonEmptyC.Reverse f) => + (a -> b -> b) -> b -> f a -> NonEmpty.T f b +nonEmptyScanr f x = + NonEmptyC.reverse . NonEmpty.scanl (flip f) x . NonEmptyC.reverse + + +{- | +Consider a superposition of all possible state sequences +weighted by the likelihood to produce the observed emission sequence. +Now train the model with respect to all of these sequences +with respect to the weights. +This is done by the Baum-Welch algorithm. +-} +trainUnsupervised :: + (Distr.Estimate tdistr, Distr.Distribution tdistr ~ distr, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) => + T distr prob -> NonEmpty.T [] emission -> Trained tdistr prob +trainUnsupervised hmm xs = + let (alphas, betas) = alphaBeta hmm xs + zetas = zetaFromAlphaBeta alphas betas + + in Trained { + trainedInitial = NonEmpty.head zetas, + trainedTransition = + sumTransitions hmm $ xiFromAlphaBeta hmm xs alphas betas, + trainedDistribution = + Distr.accumulateEmissions $ map (zip (NonEmpty.flatten xs)) $ + List.transpose $ map Vector.toList $ NonEmpty.flatten zetas + } diff --git a/src/Math/HiddenMarkovModel/Pattern.hs b/src/Math/HiddenMarkovModel/Pattern.hs new file mode 100644 --- /dev/null +++ b/src/Math/HiddenMarkovModel/Pattern.hs @@ -0,0 +1,108 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{- | +This module provides a simple way to train +the transition matrix and initial probability vector +using simple patterns of state sequences. + +You may create a trained model using semigroup combinators like this: + +> let a = atom $ HMM.state 0 +> b = atom $ HMM.state 1 +> distr = +> Distr.DiscreteTrained $ Map.fromList $ +> ('a', Vector.fromList [1,2]) : +> ('b', Vector.fromList [4,3]) : +> ('c', Vector.fromList [0,1]) : +> [] +> in finish 2 distr $ replicate 5 $ replicate 10 a <> replicate 20 b +-} +module Math.HiddenMarkovModel.Pattern ( + T, + atom, + append, + replicate, + finish, + ) where + +import qualified Math.HiddenMarkovModel.Distribution as Distr +import qualified Math.HiddenMarkovModel as HMM +import Math.HiddenMarkovModel.Private (Trained(..)) +import Math.HiddenMarkovModel.Distribution (State(State)) + +import qualified Numeric.LinearAlgebra.Algorithms as Algo +import qualified Numeric.Container as NC +import qualified Data.Packed.Vector as Vector +import Data.Packed.Matrix (Matrix) +import Data.Packed.Vector (Vector) + +import qualified Data.Map as Map +import Data.Semigroup (Semigroup, (<>), times1p) + +import Prelude hiding (replicate) + + +newtype T prob = Cons (Int -> (State, Matrix prob, State)) + +atom :: + (NC.Container Vector prob) => + State -> T prob +atom s = Cons $ \n -> (s, NC.konst 0 (n,n), s) + + +instance (Algo.Field prob) => Semigroup (T prob) where + (<>) = append + times1p k = replicate $ fromIntegral (k-1) + + +infixl 5 `append` + +append :: + (NC.Container Vector prob) => + T prob -> T prob -> T prob +append (Cons f) (Cons g) = + Cons $ \n -> + case (f n, g n) of + ((sai, ma, sao), (sbi, mb, sbo)) -> + (sai, increment (sbi,sao) 1 $ NC.add ma mb, sbo) + +replicate :: + (NC.Container Vector prob) => + Int -> T prob -> T prob +replicate ki (Cons f) = + Cons $ \n -> + case f n of + (si, m, so) -> + let k = fromIntegral ki + in (si, increment (si,so) (k-1) $ NC.scale k m, so) + +increment :: + (NC.Container Vector a) => + (State, State) -> a -> Matrix a -> Matrix a +increment (State i, State j) x m = NC.accum m (+) [((i,j), x)] + + +finish :: + (NC.Container Vector prob) => + Int -> tdistr -> T prob -> Trained tdistr prob +finish n tdistr (Cons f) = + case f n of + (State si, m, _so) -> + Trained { + trainedInitial = NC.assoc n 0 [(si,1)], + trainedTransition = m, + trainedDistribution = tdistr + } + + +_example :: HMM.DiscreteTrained Double Char +_example = + let a = atom $ HMM.state 0 + b = atom $ HMM.state 1 + distr = + Distr.DiscreteTrained $ Map.fromList $ + ('a', Vector.fromList [1,2]) : + ('b', Vector.fromList [4,3]) : + ('c', Vector.fromList [0,1]) : + [] + in finish 2 distr $ replicate 5 $ replicate 10 a <> replicate 20 b diff --git a/src/Math/HiddenMarkovModel/Private.hs b/src/Math/HiddenMarkovModel/Private.hs new file mode 100644 --- /dev/null +++ b/src/Math/HiddenMarkovModel/Private.hs @@ -0,0 +1,286 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE UndecidableInstances #-} +module Math.HiddenMarkovModel.Private where + +import qualified Math.HiddenMarkovModel.Distribution as Distr +import qualified Math.HiddenMarkovModel.CSV as HMMCSV +import Math.HiddenMarkovModel.Distribution (State(State)) + +import qualified Numeric.LinearAlgebra.Algorithms as Algo +import qualified Numeric.LinearAlgebra.Util as LinAlg +import qualified Numeric.Container as NC +import qualified Data.Packed.Development as Dev +import qualified Data.Packed.Matrix as Matrix +import qualified Data.Packed.Vector as Vector +import Numeric.Container ((<>)) +import Data.Packed.Matrix (Matrix) +import Data.Packed.Vector (Vector) + +import qualified Data.NonEmpty.Class as NonEmptyC +import qualified Data.NonEmpty as NonEmpty +import qualified Data.Semigroup as Sg +import qualified Data.List as List +import Data.Traversable (Traversable, mapAccumL) +import Data.Tuple.HT (mapPair, mapFst, mapSnd, swap) + + +{- | +A Hidden Markov model consists of a number of (hidden) states +and a set of emissions. +There is a vector for the initial probability of each state +and a matrix containing the probability for switching +from one state to another one. +The 'distribution' field points to probability distributions +that associate every state with emissions of different probability. +Famous distribution instances are discrete and Gaussian distributions. +See "Math.HiddenMarkovModel.Distribution" for details. + +The transition matrix is transposed +with respect to popular HMM descriptions. +But I think this is the natural orientation, because this way +you can write \"transition matrix times probability column vector\". + +The type has two type parameters, +although the one for the distribution would be enough. +However, replacing @prob@ by @Distr.Probability distr@ +would prohibit the derived Show and Read instances. +-} +data T distr prob = + Cons { + initial :: Vector prob, + transition :: Matrix prob, + distribution :: distr + } + deriving (Show, Read) + + +emission :: + (Distr.EmissionProb distr, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) => + T distr prob -> emission -> Vector prob +emission = Distr.emissionProb . distribution + + +forward :: + (Distr.EmissionProb distr, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission, + Traversable f) => + T distr prob -> NonEmpty.T f emission -> prob +forward hmm = NC.sumElements . NonEmpty.last . alpha hmm + +alpha :: + (Distr.EmissionProb distr, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission, + Traversable f) => + T distr prob -> + NonEmpty.T f emission -> NonEmpty.T f (Vector prob) +alpha hmm (NonEmpty.Cons x xs) = + NonEmpty.scanl + (\alphai xi -> NC.mul (emission hmm xi) (transition hmm <> alphai)) + (NC.mul (emission hmm x) (initial hmm)) + xs + + +backward :: + (Distr.EmissionProb distr, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission, + Traversable f) => + T distr prob -> NonEmpty.T f emission -> prob +backward hmm (NonEmpty.Cons x xs) = + NC.sumElements $ + NC.mul (initial hmm) $ + NC.mul (emission hmm x) $ + NonEmpty.head $ beta hmm xs + +beta :: + (Distr.EmissionProb distr, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission, + Traversable f) => + T distr prob -> + f emission -> NonEmpty.T f (Vector prob) +beta hmm = + NonEmpty.scanr + (\xi betai -> NC.mul (emission hmm xi) betai <> transition hmm) + (NC.constant 1 (NC.dim $ initial hmm)) + + +alphaBeta :: + (Distr.EmissionProb distr, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission, + Traversable f) => + T distr prob -> + NonEmpty.T f emission -> + (prob, NonEmpty.T f (Vector prob), NonEmpty.T f (Vector prob)) +alphaBeta hmm xs = + let alphas = alpha hmm xs + betas = beta hmm $ NonEmpty.tail xs + recipLikelihood = recip $ NC.sumElements $ NonEmpty.last alphas + in (recipLikelihood, alphas, betas) + + + +xiFromAlphaBeta :: + (Distr.EmissionProb distr, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) => + T distr prob -> prob -> + NonEmpty.T [] emission -> + NonEmpty.T [] (Vector prob) -> + NonEmpty.T [] (Vector prob) -> + [Matrix prob] +xiFromAlphaBeta hmm recipLikelihood xs alphas betas = + zipWith3 + (\x alpha0 beta1 -> + NC.scale recipLikelihood $ + NC.mul + (NC.outer (NC.mul (emission hmm x) beta1) alpha0) + (transition hmm)) + (NonEmpty.tail xs) + (NonEmpty.init alphas) + (NonEmpty.tail betas) + +zetaFromXi :: + (Distr.Probability distr ~ prob, Num prob, NC.Product prob) => + T distr prob -> [Matrix prob] -> [Vector prob] +zetaFromXi hmm xis = + map (NC.constant 1 (Matrix.rows $ transition hmm) <>) xis + +zetaFromAlphaBeta :: + (NC.Container Vector prob) => + prob -> + NonEmpty.T [] (Vector prob) -> + NonEmpty.T [] (Vector prob) -> + NonEmpty.T [] (Vector prob) +zetaFromAlphaBeta recipLikelihood alphas betas = + fmap (NC.scale recipLikelihood) $ + NonEmptyC.zipWith NC.mul alphas betas + + +{- | +In constrast to Math.HiddenMarkovModel.reveal +this does not normalize the vector. +This is slightly simpler but for long sequences +the product of probabilities might be smaller +than the smallest representable number. +-} +reveal :: + (Distr.EmissionProb distr, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission, + Traversable f) => + T distr prob -> NonEmpty.T f emission -> NonEmpty.T f State +reveal hmm (NonEmpty.Cons x xs) = + fmap State $ + uncurry (NonEmpty.scanr Dev.at') $ + mapFst NC.maxIndex $ + mapAccumL + (\alphai xi -> + swap $ mapSnd (NC.mul (emission hmm xi)) $ + matrixMaxMul (transition hmm) alphai) + (NC.mul (emission hmm x) (initial hmm)) xs + +matrixMaxMul :: + (NC.Container Vector a) => + Matrix a -> Vector a -> (Vector Int, Vector a) +matrixMaxMul m v = + mapPair (Vector.fromList, Vector.fromList) $ unzip $ + map ((\x -> (NC.maxIndex x, NC.maxElement x)) . NC.mul v) $ + Matrix.toRows m + + + +{- | +A trained model is a temporary form of a Hidden Markov model +that we need during the training on multiple training sequences. +It allows to collect knowledge over many sequences with 'mergeTrained', +even with mixed supervised and unsupervised training. +You finish the training by converting the trained model +back to a plain modul using 'finishTraining'. + +You can create a trained model in three ways: + +* supervised training using an emission sequence with associated states, + +* unsupervised training using an emission sequence and an existing Hidden Markov Model, + +* derive it from state sequence patterns, cf. "Math.HiddenMarkovModel.Pattern". +-} +data Trained distr prob = + Trained { + trainedInitial :: Vector prob, + trainedTransition :: Matrix prob, + trainedDistribution :: distr + } + deriving (Show, Read) + + +sumTransitions :: + (NC.Container Vector e, Num e) => + T distr t -> [Matrix e] -> Matrix e +sumTransitions hmm = + foldl NC.add (NC.konst 0 $ LinAlg.size $ transition hmm) + +{- | +Baum-Welch algorithm +-} +trainUnsupervised :: + (Distr.Estimate tdistr, Distr.Distribution tdistr ~ distr, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) => + T distr prob -> NonEmpty.T [] emission -> Trained tdistr prob +trainUnsupervised hmm xs = + let (recipLikelihood, alphas, betas) = alphaBeta hmm xs + zetas = zetaFromAlphaBeta recipLikelihood alphas betas + + in Trained { + trainedInitial = NonEmpty.head zetas, + trainedTransition = + sumTransitions hmm $ + xiFromAlphaBeta hmm recipLikelihood xs alphas betas, + trainedDistribution = + Distr.accumulateEmissions $ map (zip (NonEmpty.flatten xs)) $ + List.transpose $ map Vector.toList $ NonEmpty.flatten zetas + } + + +mergeTrained :: + (Distr.Estimate tdistr, Distr.Distribution tdistr ~ distr, + Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) => + Trained tdistr prob -> Trained tdistr prob -> Trained tdistr prob +mergeTrained hmm0 hmm1 = + Trained { + trainedInitial = NC.add (trainedInitial hmm0) (trainedInitial hmm1), + trainedTransition = + NC.add (trainedTransition hmm0) (trainedTransition hmm1), + trainedDistribution = + Distr.combine + (trainedDistribution hmm0) (trainedDistribution hmm1) + } + +instance + (Distr.Estimate tdistr, Distr.Distribution tdistr ~ distr, + Distr.Probability distr ~ prob) => + Sg.Semigroup (Trained tdistr prob) where + (<>) = mergeTrained + + +toCells :: + (Distr.CSV distr, Algo.Field prob, Show prob) => + T distr prob -> [[String]] +toCells hmm = + (HMMCSV.cellsFromVector $ initial hmm) : + (HMMCSV.cellsFromMatrix $ transition hmm) ++ + [] : + (Distr.toCells $ distribution hmm) + +parseCSV :: + (Distr.CSV distr, Algo.Field prob, Read prob) => + HMMCSV.CSVParser (T distr prob) +parseCSV = do + v <- HMMCSV.parseNonEmptyVectorCells + m <- HMMCSV.parseSquareMatrixCells $ Vector.dim v + HMMCSV.skipEmptyRow + distr <- Distr.parseCells $ Vector.dim v + return $ Cons { + initial = v, + transition = m, + distribution = distr + } diff --git a/src/Math/HiddenMarkovModel/Test.hs b/src/Math/HiddenMarkovModel/Test.hs new file mode 100644 --- /dev/null +++ b/src/Math/HiddenMarkovModel/Test.hs @@ -0,0 +1,143 @@ +module Math.HiddenMarkovModel.Test where + +import qualified Math.HiddenMarkovModel as HMM +import qualified Math.HiddenMarkovModel.Normalized as Normalized +import qualified Math.HiddenMarkovModel.Private as Priv +import qualified Math.HiddenMarkovModel.Distribution as Distr + +import qualified Numeric.Container as NC +import qualified Data.Packed.Matrix as Matrix +import qualified Data.Packed.Vector as Vector +import Data.Packed.Matrix (Matrix) +import Data.Packed.Vector (Vector) + +import qualified System.Random as Rnd + +import qualified Data.NonEmpty.Class as NonEmptyC +import qualified Data.NonEmpty as NonEmpty +import qualified Data.Foldable as Fold +import qualified Data.Map as Map +import Data.NonEmpty ((!:)) + + +hmm :: HMM.Discrete Double Char +hmm = + HMM.Cons { + HMM.initial = Vector.fromList [0.1, 0.2, 0.3, 0.4], + HMM.transition = + Matrix.fromLists $ + [0.7, 0.1, 0.0, 0.2] : + [0.1, 0.6, 0.1, 0.0] : + [0.1, 0.2, 0.7, 0.0] : + [0.1, 0.1, 0.2, 0.8] : + [], + HMM.distribution = + Distr.Discrete $ Map.fromList $ + ('a', Vector.fromList [1,0,0,0]) : + ('b', Vector.fromList [0,1,0,1]) : + ('c', Vector.fromList [0,0,1,0]) : + [] + } + + +sequ :: NonEmpty.T [] Char +sequ = 'a' !: take 20 (HMM.generate hmm (Rnd.mkStdGen 42)) + +{- | +Should all be equal. +-} +sequLikelihood :: ((Double, Double), Double, NonEmpty.T [] Double) +sequLikelihood = + ((Priv.forward hmm sequ, Priv.backward hmm sequ), + exp $ Normalized.logLikelihood hmm sequ, + NonEmptyC.zipWith NC.dot + (Priv.alpha hmm sequ) + (Priv.beta hmm $ NonEmpty.tail sequ)) + +{- | +Should all be one. +-} +sequLikelihoodNormalized :: NonEmpty.T [] Double +sequLikelihoodNormalized = + let (calphas,betas) = Normalized.alphaBeta hmm sequ + in NonEmptyC.zipWith NC.dot (fmap snd calphas) betas + + +{- | +Lists should be equal, but the first list contains one less element. +-} +zetas :: + ([Vector Double], + NonEmpty.T [] (Vector Double), + NonEmpty.T [] (Vector Double)) +zetas = + let (recipLikelihood, alphas, betas) = Priv.alphaBeta hmm sequ + in (Priv.zetaFromXi hmm $ + Priv.xiFromAlphaBeta hmm recipLikelihood sequ alphas betas, + Priv.zetaFromAlphaBeta recipLikelihood alphas betas, + uncurry Normalized.zetaFromAlphaBeta $ + Normalized.alphaBeta hmm sequ) + +{- | +Quick test of zetas - result should be @(True, very small, very small)@. +-} +zetasDiff :: (Bool, Double, Double) +zetasDiff = + case zetas of + (z0,z1,z2) -> + (length z0 == length (NonEmpty.tail z1) && + length z0 == length (NonEmpty.tail z2), + maximum $ map NC.normInf $ zipWith NC.sub z0 $ NonEmpty.init z1, + NonEmpty.maximum $ fmap NC.normInf $ NonEmptyC.zipWith NC.sub z1 z2) + +{- | +Lists should be equal +-} +xis :: ([Matrix Double], [Matrix Double]) +xis = + let (recipLikelihood, alphas, betas) = Priv.alphaBeta hmm sequ + in (Priv.xiFromAlphaBeta hmm recipLikelihood sequ alphas betas, + uncurry (Normalized.xiFromAlphaBeta hmm sequ) $ + Normalized.alphaBeta hmm sequ) + +{- | +Quick test of xis - result should be @(True, very small)@. +-} +xisDiff :: (Bool, Double) +xisDiff = + case xis of + (x0,x1) -> + (length x0 == length x1, + maximum $ map (NC.normInf . Matrix.flatten) $ zipWith NC.sub x0 x1) + + +reveal :: Bool +reveal = + Normalized.reveal hmm sequ == Priv.reveal hmm sequ + + +trainUnsupervised :: + (HMM.DiscreteTrained Double Char, + HMM.DiscreteTrained Double Char) +trainUnsupervised = + (Priv.trainUnsupervised hmm sequ, + Normalized.trainUnsupervised hmm sequ) + +trainUnsupervisedDiff :: (Double, Double, (Bool, Double)) +trainUnsupervisedDiff = + case trainUnsupervised of + (hmm0,hmm1) -> + (NC.normInf $ Matrix.flatten $ NC.sub + (Priv.trainedTransition hmm0) (Priv.trainedTransition hmm1), + NC.normInf $ NC.sub + (Priv.trainedInitial hmm0) (Priv.trainedInitial hmm1), + case (Priv.trainedDistribution hmm0, Priv.trainedDistribution hmm1) of + (Distr.DiscreteTrained m0, Distr.DiscreteTrained m1) -> + (Map.size m0 == Map.size m1, + Fold.maximum $ fmap NC.normInf $ + Map.intersectionWith NC.sub m0 m1)) + + +nonEmptyScanr :: Int -> [Int] -> Bool +nonEmptyScanr x xs = + Normalized.nonEmptyScanr (-) x xs == NonEmpty.scanr (-) x xs diff --git a/src/Math/HiddenMarkovModel/Utility.hs b/src/Math/HiddenMarkovModel/Utility.hs new file mode 100644 --- /dev/null +++ b/src/Math/HiddenMarkovModel/Utility.hs @@ -0,0 +1,35 @@ +{-# LANGUAGE FlexibleContexts #-} +module Math.HiddenMarkovModel.Utility where + +import qualified Numeric.Container as NC +import Data.Packed.Vector (Vector) + +import qualified System.Random as Rnd + +import qualified Control.Monad.Trans.State as MS + + +normalizeProb :: + (NC.Container Vector a, Fractional a) => Vector a -> Vector a +normalizeProb = snd . normalizeFactor + +normalizeFactor :: + (NC.Container Vector a, Fractional a) => + Vector a -> (a, Vector a) +normalizeFactor xs = + let c = NC.sumElements xs + in (c, NC.scale (recip c) xs) + +-- see htam:Stochastic +randomItemProp :: + (Rnd.RandomGen g, Rnd.Random b, Num b, Ord b) => + [(a,b)] -> MS.State g a +randomItemProp props = + let (keys,ps) = unzip props + in do p <- MS.state (Rnd.randomR (0, sum ps)) + return $ + fst $ head $ dropWhile ((0<=) . snd) $ + zip keys $ tail $ scanl (-) p ps + +attachOnes :: (Num b) => [a] -> [(a,b)] +attachOnes = map (flip (,) 1)