diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +The MIT License (MIT) + +Copyright (c) Leopold Tal G, Vikram Verma + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/markov-processes.cabal b/markov-processes.cabal new file mode 100644 --- /dev/null +++ b/markov-processes.cabal @@ -0,0 +1,50 @@ +name: markov-processes +version: 0.0.2 +license-file: LICENSE +license: MIT +synopsis: Hidden Markov processes. +author: Leopold Tal G , Yorick van Pelt , Vikram Verma +maintainer: Vikram Verma +category: AI +build-type: Simple +cabal-version: >=1.10 + +library + exposed-modules: + AI.Markov.HMM, + Data.Distribution + other-modules: + Data.Functor.Extras, + Data.List.Extras, + System.Random.Extras + other-extensions: + RecordWildCards, + TemplateHaskell + hs-source-dirs: + src + build-depends: + base == 4.*, + bifunctors, + memoize, + MonadRandom, + random + default-language: + Haskell2010 + +test-suite markov-tests + type: + exitcode-stdio-1.0 + main-is: + Main.lhs + hs-source-dirs: + examples, + test + build-depends: + base == 4.*, + bifunctors, + assertions, + markov-processes, + memoize, + random + default-language: + Haskell2010 diff --git a/src/AI/Markov/HMM.lhs b/src/AI/Markov/HMM.lhs new file mode 100644 --- /dev/null +++ b/src/AI/Markov/HMM.lhs @@ -0,0 +1,296 @@ +> {-# LANGUAGE RecordWildCards, TemplateHaskell #-} +> +> module AI.Markov.HMM (HMM(..), observe, evaluate, inspect, sequenceP) where +> +> import Control.Applicative ((<$>), pure) +> import Control.Monad (forM) +> import Data.Bifunctor (Bifunctor(first)) +> import Data.Distribution (Distribution(..), Probability, (), (<~~)) +> import Data.Function (on) +> import Data.Function.Memoize (Memoizable(..), deriveMemoize, deriveMemoizable, memoize4) +> import Data.List (maximumBy) +> import Data.List.Extras (pairs, argmax) +> import Data.Ratio (Ratio) +> import System.Random (RandomGen(..)) +> import System.Random.Extras (split3) + +_Hidden Markov models_ (HMMs) are used to model generative sequences +characterisable by a doubly-embedded stochastic process in which the +underlying process is hidden, and can only be observed through an +upper-level process, which produces output. + +Structurally: + +You know what a [Mealy machine][mealy] is, right? A discrete hidden Markov +model is structurally similar to a Mealy machine, except that its transitions +and output values are governed by probability distributions, rather than tokens +from an input alphabet. Rather than having a fixed initial state, this is +chosen by an experiment on the prior distribution for the state set; a set of +conditional probability distributions determine the likelihood of transition to +any state from a given state, and another set the conditional probability of +emitting each symbol in the output alphabet. + + [mealy]: https://en.wikipedia.org/wiki/Mealy_machine + +More formally, a HMM is a five-tuple consisting of: + +> data HMM state symbol = HMM + + 1. The distinct states of the process, S. (|S| = N) + +> { states :: [state] + + 2. A finite dictionary of possible observations, E. + +> , symbols :: [symbol] + + 3. An initial state distribution, the likelihood of the process + starting in each state s ∈ S. + +> , start :: Distribution state + + 4. A transition distribution; given some current state i ∈ S, this + provides the likelihood of the process next transiting to any + j ∈ S. By the Markov assumption, this probability is dependent only + on i. + +> , transition :: state -> Distribution state + + 5. The observation symbol distributions; given some s ∈ S, this + provides the likelihood of every observation o ∈ E being observed + at that i. By the independence assumption, the output observation + at any time is dependent only on the current state. + +> , emission :: state -> Distribution symbol +> } + +Having characterised some sequence generator as a HMM (and so given +a well-parametrised configuration of such a model), we can simulate +its output without collecting any further data; inferring future +behaviour from demonstrated statistical properties. This is particularly +useful in cases where gathering raw data is expensive. + +> observe :: RandomGen seed => seed -> HMM state symbol -> [symbol] +> observe seed hmm@HMM{..} = observe' ns hmm (start <~~ ts) +> where (ts, ns) = split seed +> +> observe' :: RandomGen seed => seed -> HMM state symbol -> state -> [symbol] +> observe' seed hmm@HMM{..} state = obs : observe' s2 hmm nxt +> where (s0,s1,s2) = split3 seed +> (obs, nxt) = (emission state <~~ s0, transition state <~~ s1) + +Rabiner [^rabiner1989] outlined three fundamental inference problems for HMMs: + + 1. Evaluation: given a model and a sequence $(O_n)_{n=1}^T$ of observations, compute the + likelihood those observations were produced by that HMM. This can also be + interpreted as a scoring problem -- given a sequence produced by the real + signal source, we can compare the accuracy of models. + +We first consider a straightforward (albeit intractable) approach; +computing the likelihood our observations were produced by each possible +state sequence $(I_n)_{n=1}^T$ (of appropriate length), and summing the result +probabilities. I.e: $$P(O|HMM)=\sum_{j=0}^{n}P(O|I_n,HMM)P(I_n|HMM)$$ + +First, some precursors; given a set of states, `sequencesOfN` finds all +$n$-length state sequences: + +> sequencesOfN :: Int -> HMM state symbol -> [[state]] +> sequencesOfN n = sequence . replicate n . states + +`sequenceP` determines the likelihood of a state sequence, given a model; +$P(I|HMM)=P(I_1)P(I_2|I_1)\ldots P(I_r|I_{r-1})$: + +> sequenceP :: Eq state => HMM state symbol -> [state] -> Probability +> sequenceP HMM{..} sequence = product +> $ head sequence : map (uncurry (?>) . first transition) (pairs sequence) + +`sequenceObservationsP` computes the likelihood of an observation sequence +given a state sequence; $P(O|I, HMM)=P(O_1|I_1)P(O_2|I_2) +\ldots P(O_n|I_n)$: + +> sequenceObservationsP :: Eq symbol => HMM state symbol -> [(state, symbol)] -> Probability +> sequenceObservationsP HMM{..} = product . map (uncurry (?>) . first emission) + +Using the above primitives, we can now express the procedure: + +> inefficientEvaluate :: (Eq state, Eq symbol) => HMM state symbol -> [symbol] -> Probability +> inefficientEvaluate hmm observations = sum $ zipWith (*) statesP statesObsP +> where states = sequencesOfN (length observations) hmm +> statesP = sequenceP hmm <$> states +> statesObsP = sequenceObservationsP hmm <$> map (`zip` observations) states + +But this has abysmal runtime-performance! Let $N$ be the number of +states, and $T$ equal to the number of observations (and thus the +length of each state sequence); there are thus $N^T$ possible state +sequences, and for each sequence we require about $2T$ calculations +(Rabiner's [^rabiner1989] figure, I haven't validated this); meaning +$2TN^T$ calculations in total: + + $N$ $T$ $2TN^T$ + ---- ---- ------------------ + 5 100 $\approx 10^{72}$ + 10 100 $\approx 10^{102}$ + 15 100 $\approx 10^{120}$ + +Given this profile, this function is only included here for didactic +purposes, and is not exported by this module. + +A more efficient solution exists in the _forward algorithm_, which +arranges the computation so that redundant calculations may be cached: + +Given a partial observation sequence $\bold{O} = {O_1,O_2,\ldots,O_n}$ and +a terminal state $i$, the _forward variable_ provides the likelihood of +having observed $\bold{O}$ and being in state $i$ after time $n$ -- +$\alpha_n(i) = P(O_1,O_2,\ldots,O_n,I_n=i|HMM)$: + +> forwardVariable' :: (Memoizable state, Memoizable symbol, Eq state, Eq symbol, Enum state, Bounded state) => Int -> HMM state symbol -> [symbol] -> state -> Probability +> forwardVariable' 0 HMM{..} observations state = (state forwardVariable' t hmm@HMM{..} observations state = (*) (observations !! t predecessor <- states +> let a = forwardVariable (t-1) hmm observations predecessor +> b = state return $ a * b + +Below, we compute the terminal forward variable for each state in the HMM, +using a constant set of observations. As the computation of `a` is independent +of the state under consideration, it's time-saving to memoise this value: + +> forwardVariable :: (Enum state, Bounded state, Eq state, Eq symbol, Memoizable state, Memoizable symbol) => Int -> HMM state symbol -> [symbol] -> state -> Probability +> forwardVariable = memoize4 forwardVariable' +> +> instance (Memoizable state, Memoizable symbol, Eq state, Eq symbol, Enum state, Bounded state) => Memoizable (HMM state symbol) where +> memoize = $(deriveMemoize ''HMM) +> +> deriveMemoizable ''Ratio + +As the final observation must be emitted in some (unknown) state, taking the +sum of terminal forward variables for each state yields our desired probability +-- $P(O_1,O_2,\ldots,O_n|HMM) = \sum_{n=1}^{|I|} \alpha_{|\bold{O}|}(I_n)$: + +> forward :: (Memoizable state, Memoizable symbol, Eq state, Eq symbol, Enum state, Bounded state) => HMM state symbol -> [symbol] -> Probability +> forward hmm@HMM{..} observations = sum [forwardVariable t hmm observations state | state <- states] +> where t = pred $ length observations +> +> evaluate :: (Memoizable state, Memoizable symbol, Eq state, Eq symbol, Enum state, Bounded state) => HMM state symbol -> [symbol] -> Probability +> evaluate = forward + +This evaluation algorithm requires about $TN^2$ calculations, making it +many orders of magnitudes more efficient than the naïve approach: + + $N$ $T$ $TN^2$ + ---- ---- ------- + 5 100 2,500 + 10 100 10,000 + 15 100 22,500 + + 2. Inspection: uncovering the hidden part of the model; from an + observation sequence and a model, uncover the state sequence best + explaining those observations, as provided by some optimality + criterion. + +There are several possible ways of solving this problem: there are +several possible optimality criteria, as the specification is ambiguous +as to the definition of an optimal state sequence. + +One possible criterion involves maximising the expected number of +correct individual states: computing at each point in time, the most +likely state given the observation sequence and model. + +To implement this solution, we need implement a variant of the forward +variable (namely, the backward variable), that describes the likelihood +of observing a sequence of succeeding observations from some known state +-- $\beta_n(i) = P(O_{n+1},O_{n+2},\ldots,O_{T},I_n=i,HMM)$: + +> backwardVariable' :: (Memoizable state, Memoizable symbol, Eq state, Eq symbol, Enum state, Bounded state) => Int -> HMM state symbol -> [symbol] -> state -> Probability +> backwardVariable' n HMM{..} observations state = if succ n == length observations then 1 else sum $ do +> successor <- states +> let a = transition state ?> successor +> b = emission successor ?> (observations !! n) +> c = backwardVariable (n+1) HMM{..} observations successor +> return $ a * b * c +> +> backwardVariable :: (Memoizable state, Memoizable symbol, Eq state, Eq symbol, Enum state, Bounded state) => Int -> HMM state symbol -> [symbol] -> state -> Probability +> backwardVariable = memoize4 backwardVariable' + +To recap: the forward variable determines the likelihood of +reaching some state $i$ being reached at time $n$, and a sequence +$O_1,O_2,\ldots,O_n$ being observed preceding it. The backward +variable accounts for the likelihood of the same state $i$ being +reached at time $n$, and the succeeding observation sequence being +$O_{n+1},O_{n+2},\ldots,O_{T}$. + +Between them we can compute, provided a model and observation sequence, +the likelihood of some state $i$ given the observation sequence--the +smoothed probability value: $\gamma_n(i)=\frac{\alpha_n(i)\beta_n(i)} +{\sum^{N}_{s=1}\alpha_n{s}\beta_n{s}} = \frac{P(I_n=i,O_1,O_2,\ldots,O_T|HMM)} +{\sum^{N}_{s=1} P(I_n=s,O_1,O_2,\ldots,O_T|HMM)} = +\frac{P(I_n=i,O_1,O_2,\ldots,O_T|HMM)} {P(O_1,O_2,\ldots,O_T|HMM)} = +P(I_n=i|O_1,O_2,\ldots,O_T,HMM)$ + +> forwardBackwardVariable :: (Memoizable state, Memoizable symbol, Eq state, Eq symbol, Enum state, Bounded state) => Int -> HMM state symbol -> [symbol] -> state -> Probability +> forwardBackwardVariable n hmm observations state = forwardVariable n hmm observations state +> * backwardVariable n hmm observations state +> +> smooth :: (Memoizable state, Memoizable symbol, Eq state, Eq symbol, Enum state, Bounded state) => Int -> HMM state symbol -> [symbol] -> state -> Probability +> smooth n hmm@HMM{..} observations state = numerator / denominator +> where numerator = forwardBackward n state +> denominator = sum $ zipWith forwardBackward [0..] states +> forwardBackward n = forwardBackwardVariable n hmm observations + +By maximising the smoothing value at each position in the sequence, we +can find the most likely state at each position: + +> forwardBackward :: (Memoizable state, Memoizable symbol, Eq state, Eq symbol, Enum state, Bounded state) => HMM state symbol -> [symbol] -> [state] +> forwardBackward hmm@HMM{..} observations = do +> position <- [0..pred $ length observations] +> return $ argmax (smooth position hmm observations) states + +This is a bad criterion, though; in considering only individual states, +we neglect information about the probability that state sequences will +occur. Consider the case in which a HMM has state transitions with zero +probability: the optimal state sequence may not even be valid! For this +reason, as before, this function is not exported by this module. + +A more reasonable optimality criterion is to find the most probable +contiguous sequence of states; i.e. determining the state sequence that +maximises $P(O,I|HMM)$. This can be found by application of _the Viterbi +algorithm_, which involves maximising over likelihood estimates for each +possible state sequence. + +Given an observation sequence $\bold{O} = {O_1,O_2,\ldots,O_T}$ a state +$i$, the _Viterbi step_ finds the state sequence most likely to account +for the first $n$ observations and terminating at state $i$; that is: +$$\delta_n(i) = \argmax_{I_1,I_2,\ldots,I_{n-1}} P(I_1,I_2,\ldots,I_n=1, +\bold{O}|HMM)$$: + +> viterbiStep' :: (Memoizable state, Memoizable symbol, Eq state, Eq symbol, Enum state, Bounded state) => Int -> HMM state symbol -> [symbol] -> state -> ([state], Probability) +> viterbiStep' 0 hmm@HMM{..} observations state = (pure state, start ?> state * emission state ?> head observations) +> viterbiStep' n hmm@HMM{..} observations state = argmax snd $ do +> predecessor <- states +> let (path, prob) = viterbiStep (n-1) hmm observations predecessor +> likelihood = prob * transition predecessor ?> state * emission state ?> (observations !! n) +> return $ (state:path, likelihood) +> +> viterbiStep :: (Memoizable state, Memoizable symbol, Eq state, Eq symbol, Enum state, Bounded state) => Int -> HMM state symbol -> [symbol] -> state -> ([state], Probability) +> viterbiStep = memoize4 viterbiStep' + +As the final observation must be emitted in some state, maximising the +Viterbi step over $S$ yields the desired state sequence: + +> viterbi :: (Memoizable state, Memoizable symbol, Eq state, Eq symbol, Enum state, Bounded state) => HMM state symbol -> [symbol] -> [state] +> viterbi hmm@HMM{..} observations = reverse . fst . argmax snd $ do +> let t = length observations - 1 +> terminal <- states +> return $ viterbiStep t hmm observations terminal +> +> inspect :: (Memoizable state, Memoizable symbol, Eq state, Eq symbol, Enum state, Bounded state) => HMM state symbol -> [symbol] -> [state] +> inspect = viterbi + + 3. Training: given some observation sequence, determine the parameters of + some HMM that best model the data. If we so adapt model parameters to + observed training data, we can accurately simulate real signal sources. + + [^rabiner1989]: Lawrence R. Rabiner, _A Tutorial on Hidden Markov Models and + Selected Applications in Speech Recognition_, Proceedings of + the IEEE 77 (1989): 257-286 diff --git a/src/Data/Distribution.lhs b/src/Data/Distribution.lhs new file mode 100644 --- /dev/null +++ b/src/Data/Distribution.lhs @@ -0,0 +1,47 @@ +> module Data.Distribution (Distribution(..), Probability, probability, (), choose, (<~~), chooseMany) where +> +> import Data.Maybe (fromMaybe) +> import Data.Functor.Extras ((<$$>)) +> import Data.Ratio ((%)) +> import Control.Monad.Random (evalRand, fromList) +> import System.Random (RandomGen(..)) + +We represent probabilities as ratios of `Integer`s , with no strict +bounds-checking: + +> type Probability = Rational + +Distribution represents a discrete probability distribution. More +specifically: given some type parameter, a Distribution relates each +value of that type with the probability of its occurrence. + +> type Distribution a = [(a, Probability)] + +`choose` takes a random seed and a Distribution, and produces a value +according to that Distribution. + +> choose :: RandomGen g => g -> Distribution a -> a +> choose seed = (`evalRand` seed) . fromList +> +> (<~~) :: RandomGen g => Distribution a -> g -> a +> (<~~) = flip choose + +`chooseMany` takes a random seed and a Distribution, and produces an +infinite list of values according to that Distribution. Exercise: why +shouldn't this be implemented as `map choose . repeat . fromList`? + +> chooseMany :: RandomGen g => g -> Distribution a -> [a] +> chooseMany seed = (`evalRand` seed) . sequence . repeat . fromList + +`probability` takes a Distribution and a value, and returns the +likelihood of that value occurring per that distribution. Should a +probability not be defined for that value, this likelihood is zero: + +> probability :: Eq a => a -> Distribution a -> Probability +> probability = fromMaybe (0%1) <$$> lookup +> +> ( a -> Distribution a -> Probability +> ( +> (?>) :: Eq a => Distribution a -> a -> Probability +> (?>) = flip probability diff --git a/src/Data/Functor/Extras.hs b/src/Data/Functor/Extras.hs new file mode 100644 --- /dev/null +++ b/src/Data/Functor/Extras.hs @@ -0,0 +1,8 @@ +module Data.Functor.Extras ((<$$>), for) where + +(<$$>) :: (Functor f, Functor g) => (a -> b) -> f(g(a)) -> f(g(b)) +(<$$>) = fmap fmap fmap +infixl 4 <$$> + +for :: Functor f => f a -> (a -> b) -> f b +for = flip fmap diff --git a/src/Data/List/Extras.lhs b/src/Data/List/Extras.lhs new file mode 100644 --- /dev/null +++ b/src/Data/List/Extras.lhs @@ -0,0 +1,16 @@ +> module Data.List.Extras (pairs, argmax) where +> +> import Control.Applicative ((<*>)) +> import Data.Function (on) +> import Data.List (maximumBy) + +`pairs` takes a list, and returns a list of consecutive elements pairs: + +> pairs :: [a] -> [(a,a)] +> pairs = zip <*> tail + +`argmax` takes a function and a list, and returns the element of the +list for which the given function attains its maximum value: + +> argmax :: Ord b => (a -> b) -> [a] -> a +> argmax = maximumBy . on compare diff --git a/src/System/Random/Extras.lhs b/src/System/Random/Extras.lhs new file mode 100644 --- /dev/null +++ b/src/System/Random/Extras.lhs @@ -0,0 +1,11 @@ +> module System.Random.Extras (split3) where +> +> import System.Random (RandomGen(..)) + +Obtain a vector of three distinct random generators from one. TODO: +learn implementation's randomness properties. + +> split3 :: RandomGen g => g -> (g, g, g) +> split3 seed = (b, c, d) +> where (a, b) = split seed +> (c, d) = split a diff --git a/test/Main.lhs b/test/Main.lhs new file mode 100644 --- /dev/null +++ b/test/Main.lhs @@ -0,0 +1,44 @@ +> module Main (main) where +> +> import AI.Markov.HMM (HMM(..), sequenceP, evaluate, inspect) +> import Control.Applicative ((<$>)) +> import Data.Bifunctor (Bifunctor(first)) +> import Example.Clinic (clinic, Health(..), Symptom(..)) +> import Test.Assert (runAssertions) + +`sequenceP` should determine the likelihood of a state sequence, given a model; +$P(I|HMM)=P(i_1)P(i_2|i_1)\ldots P(i_r|i_{r-1})$: + +> sequencePTests :: [(String, Bool)] +> sequencePTests = let test = sequenceP clinic in first ("sequenceP: " ++) <$> +> [ ("Properly evaluates start probabilities." , test [Healthy] == 0.6) +> , ("Properly evaluates transition probabilities." , test [Healthy, Fever] == 0.6*0.3) +> ] + +The evaluate function should predict the likelihood some sequence of +observations were produced by that HMM; $P(O_1,O_2,\ldots,O_n|HMM)$. + +> evaluateTests :: [(String, Bool)] +> evaluateTests = let test = evaluate clinic in first ("evaluate: " ++) <$> +> [ ("Likelihood of [N] for clinic." , test [Normal] == 0.6*0.5+0.4*0.1) +> , ("Likelihood of [N,N] for clinic." , test [Normal,Normal] == 0.6*0.5*0.7*0.5+0.6*0.5*0.3*0.1+0.4*0.1*0.4*0.5+0.4*0.1*0.6*0.1) +> ] + +The inspect function should find the state sequence most likely +responsible for some sequence of observations, given a HMM; +$argmax_{q\inQ}P(q,O_1,O_2,\ldots,O_n|HMM)$: + +> inspectTests :: [(String, Bool)] +> inspectTests = let test = inspect clinic in first ("inspect: " ++) <$> +> [ ("Sequence behind [N] for clinic.", test [Normal] == [Healthy]) +> , ("Sequence behind [N,C,D] for clinic.", test [Normal, Cold, Dizzy] == [Healthy, Healthy, Fever]) +> ] + +To run our assertions: + +> main :: IO () +> main = runAssertions $ concat +> [ sequencePTests +> , evaluateTests +> , inspectTests +> ]