Haskell implementations of Hidden Markov Models & related algorithms from AIMA book
1{-# LANGUAGE DataKinds #-}
2{-# LANGUAGE TypeOperators #-}
3{-# LANGUAGE RankNTypes #-}
4{-# LANGUAGE ScopedTypeVariables #-}
5---------------------------------
6-- |
7-- Module : CHMM
8-- Copyright : (C) Kaushik Chakraborty, 2019
9-- License : Apache v2 (see the file LICENSE)
10-- Maintainer : Kaushik Chakraborty <git@kaushikc.org>
11-- Stability : experimental
12--
13-- Cost based HMM
14-- https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0202337
15---------------------------------
16
17module CHMM where
18
19import Papa
20import GHC.TypeNats
21import Data.Proxy (Proxy(..))
22import Linear.V
23import Linear
24import Combinatorics
25import qualified Data.Vector as Vec
26
27import HMM
28-----------------------
29
30-- * Documentation
31
32-- | For a specific resiliency property that need evaluation, following algorithm infer the probability of the same property in linear time taking into account corresponding observations and state trajectories
33cHMMInfer :: forall s u t b i. (KnownNat s, KnownNat u, KnownNat t, Eq b, Ord b, Integral i, Show i, Show b)
34 => HMM s u R b -- ^ HMM model
35 -> V t b -- ^ evidence/observation trajectory for time @1 .. t@
36 -> V t i -- ^ state trajectory for time @1 .. t@; represented as @0@-based index of state
37 -> R -- ^ property probability
38cHMMInfer hmm evs xs = (\x y z -> (x * y) / z) <$> y0 <*> y1 <*> y2 $ zs
39 where
40 y0 :: Vec.Vector (i, b) -> R
41 y0 = foldl' (\m (s,o') ->
42 m *
43 sensorDiagonal hmm o' ^?! _Just . ix (fromIntegral s) . ix (fromIntegral s))
44 1
45
46 y1 :: Vec.Vector (i, b) -> R
47 y1 vs = ifoldl (\idx m (s,_) ->
48 m *
49 bool
50 (hmm ^?! tModel . ix (vs ^?! ix (idx - 1) . _1 . to fromIntegral) . ix (fromIntegral s))
51 ((Linear.transpose (hmm ^. prior) !*! hmm ^. tModel) ^?! ix 0 . ix (fromIntegral s))
52 (idx == 0))
53 1 vs
54
55 y2 :: KnownNat s => Vec.Vector (i, b) -> R
56 y2 vs = sum . extract1 $
57 ifoldl (\idx m (_,o') ->
58 bool
59 ((sensorDiagonal hmm o' ^?! _Just) !*! Linear.transpose (hmm ^. tModel) !*! m)
60 ((sensorDiagonal hmm o' ^?! _Just) !*! Linear.transpose (hmm ^. tModel) !*! hmm ^. prior)
61 (idx == 0)
62 )
63 (unitColumn :: Message s R)
64 vs
65
66 zs :: Vec.Vector (i, b)
67 zs = Vec.zip (toVector xs) (toVector evs)
68
69-- | a trajectory is a sequence of variables which can be random variables (r.v) as part of either sensor or transition model of a HMM or values from the support of some r.v.
70trajectory :: forall n t a . (KnownNat n, KnownNat t)
71 => Proxy t -- ^ choice parameter
72 -> V n a -- ^ superset of choice variables which is divided into sets each of length provided by the choice parameter
73 -> V (n ^ t) (V t a) -- ^ final set with @n^t@ ordered subsets each of length @t@
74trajectory p = V . Vec.fromList . fmap (V . Vec.fromList) . variateRep t' . toList
75 where
76 t' = fromIntegral $ natVal p
77
78-- | HMM model representing the swarm robotics use case of the paper
79--
80-- Transition Model
81--
82-- @
83-- [0.7393859 , 0.2233470 , 0.0330670 , 0.0042001]
84-- [0.0928202 , 0.7172515 , 0.1719556 , 0.0179727]
85-- [0.0139341 , 0.1755173 , 0.6936519 , 0.1168966]
86-- [0.0047083 , 0.0403520 , 0.2561893 , 0.6987504]
87-- @
88--
89-- Sensor Model
90--
91-- @
92-- [1.00000, 0.00000, 0.00000, 0.00000]
93-- [0.10000, 0.90000, 0.00000, 0.00000]
94-- [0.01000, 0.18000, 0.81000, 0.00000]
95-- [0.00100, 0.02700, 0.24300, 0.72900]
96-- @
97--
98swarmRobotsHMM :: HMM 4 4 R Int
99swarmRobotsHMM = mkHMM_
100 (Just $ V $ Vec.replicate 4 (toV $ V1 0.25))
101 (V $ Vec.fromList
102 [
103 V $ Vec.fromList [0.7393859 , 0.2233470 , 0.0330670 , 0.0042001],
104 V $ Vec.fromList [0.0928202 , 0.7172515 , 0.1719556 , 0.0179727],
105 V $ Vec.fromList [0.0139341 , 0.1755173 , 0.6936519 , 0.1168966],
106 V $ Vec.fromList [0.0047083 , 0.0403520 , 0.2561893 , 0.6987504]
107 ])
108 (V $ Vec.fromList [1..4])
109 (Linear.transpose (V $ Vec.fromList
110 [
111 V $ Vec.fromList [1.00000, 0.00000, 0.00000, 0.00000],
112 V $ Vec.fromList [0.10000, 0.90000, 0.00000, 0.00000],
113 V $ Vec.fromList [0.01000, 0.18000, 0.81000, 0.00000],
114 V $ Vec.fromList [0.00100, 0.02700, 0.24300, 0.72900]
115 ]))
116
117-- | inferring the /l-resistance/ property of the swarm robots as represented in the 'swarmRobotsHMM' model given a series of observation and a paticular /l/ value to check against
118--
119-- E.g.
120--
121-- >>> runCHMMInfer [3,2,3] 2
122-- 0.7756162284376005
123--
124runCHMMInfer :: [Int] -> Int -> R
125runCHMMInfer evs l = let v = Vec.fromList [0 .. dim (swarmRobotsHMM ^. tModel) - 1]
126 n = length evs
127 in
128 reifyDimNat n $ \t' ->
129 reifyVectorNat v $ \n' ->
130 let ts = trajectory t' n'
131 valid_ts = V $ Vec.filter (Vec.all (>=l) . toVector) $ toVector ts
132 in
133 sum $ cHMMInfer swarmRobotsHMM (V $ Vec.fromList evs) <$> valid_ts
134
135
136
137
138
139