markov-realization 0.1.0 → 0.2.1

raw patch · 8 files changed

+635/−565 lines, 8 filesdep +HTFdep +comonaddep +configuration-toolsdep ~MonadRandomdep ~contravariantdep ~discriminationPVP ok

version bump matches the API change (PVP)

Dependencies added: HTF, comonad, configuration-tools, markov-realization

Dependency ranges changed: MonadRandom, contravariant, discrimination, generic-deriving

API changes (from Hackage documentation)

Files

• README.md +11/−13
• markov-realization.cabal +53/−30
• src/Markov.hs +89/−91
• src/Markov/Example.hs +417/−0
• src/Markov/Examples.hs +0/−422
• src/Markov/Instance.hs +8/−0
• src/Markov/Instances.hs +0/−9
• test/Test.hs +57/−0

@@ -1,6 +1,6 @@ # Markov Tutorial -Let X<sub>n</sub> denote the nth state of a Markov chain with state space ℕ.+Let Xₙ denote the nth state of a Markov chain with state space ℕ. For x ≠ 0 define transition probabilities  p(x,0) = q,@@ -14,7 +14,7 @@ p(x,x+1) = s. Let p(x,y) = 0 in all other cases. Suppose we wanted to find-P\[X<sub>n</sub> = j ∩ d = k],+P\[Xₙ = j ∩ d = k], where d denotes the number of transitions from a positive integer to zero. There are three values we need to track — extinctions, probability, and state.@@ -49,17 +49,15 @@  We can now easily see a list of states, deaths, and the probabilities. -<pre>-<b>-> chain [pure 0 :: Sum Int :* Product Rational :* Extinction] !! 3-</b>-((0,8 % 125),0)-((0,111 % 500),1)-((1,51 % 500),0)-((0,9 % 25),2)-((1,9 % 250),1)-((0,27 % 125),3-</pre>+__`> chain [pure 0 :: Sum Int :* Product Rational :* Extinction] !! 3`__+```+[ ((0,8 % 125),0)+, ((0,111 % 500),1)+, ((1,51 % 500),0)+, ((0,9 % 25),2)+, ((1,9 % 250),1)+, ((0,27 % 125),3) ]+```  This means that starting from a state of zero, after three time steps there is a 51/500 chance

@@ -1,38 +1,61 @@ cabal-version: 1.12-name: markov-realization-version: 0.1.0-license: BSD3-license-file: LICENSE-copyright: 2019 Alex Loomis-maintainer: atloomis@math.arizona.edu-author: Alex Loomis-homepage: https://github.com/alexloomis/markov-bug-reports: https://github.com/alexloomis/markov/issues-synopsis: Realizations of Markov chains.-description:-    Please see the README on GitHub at <https://github.com/alexloomis/markov#markov-tutorial>-category: Statistics-build-type: Simple++-- This file has been generated from package.yaml by hpack version 0.31.1.+--+-- see: https://github.com/sol/hpack+--+-- hash: b164fe328e9edef0858d48e61b56997280bd4bb2ea6ffe923afb24084a14efe6++name:           markov-realization+version:        0.2.1+description:    Please see the README on GitHub at <https://github.com/alexloomis/markov#markov-tutorial>+homepage:       https://github.com/alexloomis/markov+bug-reports:    https://github.com/alexloomis/markov/issues+author:         Alex Loomis+maintainer:     atloomis@math.arizona.edu+copyright:      2019 Alex Loomis+license:        BSD3+license-file:   LICENSE+build-type:     Simple extra-source-files:     README.md     ChangeLog.md+category:       Statistics+synopsis:       Realizations of Markov chains.  source-repository head-    type: git-    location: git://github.com/alexloomis/markov.git+  type: git+  location: git://github.com/alexloomis/markov.git  library-    exposed-modules:-        Markov-        Markov.Examples-        Markov.Instances-    hs-source-dirs: src-    other-modules:-        Paths_markov_realization-    default-language: Haskell2010-    build-depends:-        base >=4.7 && <5,-        contravariant >=1.5.1 && <1.6,-        discrimination ==0.4.*,-        generic-deriving >=1.12.4 && <1.13,-        MonadRandom >=0.5.1.1 && <0.6+  exposed-modules:+      Markov+      Markov.Example+      Markov.Instance+  other-modules:+      Paths_markov_realization+  hs-source-dirs:+      src+  build-depends:+      base >=4.7 && <5+    , comonad+    , configuration-tools+    , contravariant+    , discrimination+    , generic-deriving+    , MonadRandom+  default-language: Haskell2010++test-suite markov-test+  type: exitcode-stdio-1.0+  main-is: Test.hs+  other-modules:+      Paths_markov_realization+  hs-source-dirs:+      test+  ghc-options: -threaded -rtsopts -with-rtsopts=-N+  build-depends:+      base >=4.7 && <5+    , HTF+    , markov-realization+  default-language: Haskell2010

@@ -1,6 +1,11 @@-{-# LANGUAGE MultiParamTypeClasses, GeneralizedNewtypeDeriving, DeriveGeneric,-DeriveAnyClass, DerivingStrategies, TypeOperators #-}-{-|+{-# LANGUAGE DeriveGeneric              #-}+{-# LANGUAGE DeriveAnyClass             #-}+{-# LANGUAGE DerivingStrategies         #-}+{-# LANGUAGE FlexibleContexts           #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE MultiParamTypeClasses      #-}+{-# LANGUAGE TypeOperators              #-}+{- | Module      : Markov Description : Realization of Markov processes with known parameters. Maintainer  : atloomis@math.arizona.edu@@ -9,42 +14,45 @@ Three type classes for deterministically analyzing Markov chains with known parameters. 'Markov0' is intended to list possible outcomes,-'Markov' should allow for more sophisticated analysis,-and 'MultiMarkov' is intended to make implementing-hidden Markov models easier.+'Markov' should allow for more sophisticated analysis. See "Examples" for examples. See README for a detailed description. -} module Markov (-              -- *Markov0-                Markov0 (..)-              -- *Markov-              , Markov (..)-              -- *MultiMarkov-              , randomProduct-              , randomPath-              , MultiMarkov (..)-              -- *Combine-              , Combine (..)-              , Merge (..)-              , Sum (..)-              , Product (..)-              -- *Misc-              , (:*)-              , (>*<)-              , fromLists-              -- *Testing-              ) where+     -- *Markov0+       Markov0 (..)+     , chain0 -import Markov.Instances ()-import Control.Applicative ((<**>))-import Generics.Deriving (Generic)+     -- *Markov+     , Markov (..)+     , chain++     -- *Combine+     , Combine (..)+     , Merge (..)+     , Sum (..)+     , Product (..)++     -- *Misc+     , (:*)+     , (>*<)+     , fromLists+     , randomProduct+     , randomPath+     ) where++-- import Configuration.Utils.Operators ((<*<))+import Control.Comonad import Data.Discrimination (Grouping, grouping)+import Generics.Deriving (Generic)++import Markov.Instance ()++import qualified Control.Monad.Random as MR import qualified Data.Discrimination as DD+import qualified Data.Functor.Contravariant as FC import qualified Data.List as DL import qualified Data.List.NonEmpty as NE-import qualified Data.Functor.Contravariant as FC-import qualified Control.Monad.Random as MR  --------------------------------------------------------------- -- Markov0@@ -52,48 +60,43 @@  -- |A basic implementation of Markov chains. class (Eq m) => Markov0 m where-    -- |The transition functions from a state.     transition0 :: m -> [m -> m]     step0       :: m -> [m]     -- |Iterated steps.-    chain0      :: [m] -> [[m]]-    step0 x = fmap ($ x) (transition0 x)-    chain0  = DL.iterate' $ DL.nub . concatMap step0+    transition0 x = const <$> step0 x+    step0 x = ($ x) <$> transition0 x+    {-# MINIMAL transition0 | step0 #-} +chain0 :: Markov0 m => [m] -> [[m]]+chain0 = DL.iterate' $ DL.nub . concatMap step0+ --------------------------------------------------------------------------------------- -- Markov ---------------------------------------------------------------------------------------  -- |An implementation of Markov chains.--- To speed up @chain@, try instead:------ > chain = DL.iterate' $ map summarize' . NE.group . DL.sort . concatMap step--- >     where summarize' xs@((_,b)NE.:|_) = (summarize . fmap fst $ xs, b)-class (Combine t, Grouping t, Grouping m, Monoid t) => Markov t m where-    transition :: m -> [(t, m -> m)]-    step       :: (t,m) -> [(t,m)]-    chain      :: [(t,m)] -> [[(t,m)]]-    step x = fmap (x <**>) (transition $ snd x)-    -- |Iterated steps, with equal states combined using 'summarize' operation.-    chain  = DL.iterate' $ map (summarize' . NE.fromList)-             . DD.group . concatMap step-             where summarize' xs@((_,b)NE.:|_) = (summarize . fmap fst $ xs, b)-             -- WARNING: DD.group does not currently respect equivalence classes.-------------------------------------------------------------------------------------------- Multi-Transition Markov----------------------------------------------------------------------------------------+class (Applicative t, Comonad t) => Markov t m where+    transition :: m -> [t (m -> m)]+    step       :: t m -> [t m]+    sequential :: [m -> [t (m -> m)]]+    transition = fmap (fmap const) . step . pure+    step x = foldr (concatMap . step') [x] sequential+      where step' f y = (<*> y) <$> f (extract y)+    sequential = [transition]+    {-# MINIMAL transition | step | sequential #-}+    -- Could also be defined as follows:+    --+    -- transition = foldr compose stayPut sequential+      -- where stayPut = const [pure id]+            -- compose g f a = composeWith g a =<< f a+            -- composeWith g a x = (<*< x) <$> g (extract $ fmap ($ a) x)+    -- step x = (<*> x) <$> transition (extract x)+    -- sequential = [fmap (fmap const) . step . pure] --- |An implementation of Markov chains that allows multi-transition steps.-class (Combine m, Grouping m, Semigroup m) => MultiMarkov m where-    multiTransition :: m -> [m -> [m]]-    multiStep       :: m -> [m]-    multiChain      :: [m] -> [[m]]-    multiStep x = foldr phi [x] (multiTransition x)-        where phi f = concatMap (delta f)-              delta f y = map (y <>) (f y)-    multiChain  = DL.iterate' $ map (summarize . NE.fromList)-                  . DD.group . concatMap multiStep+-- WARNING: DD.group does not currently respect equivalence classes.+-- |Iterated steps, with equal states combined using 'summarize' operation.+chain :: (Combine (t m), Grouping (t m), Markov t m) => [t m] -> [[t m]]+chain = DL.iterate' $ fmap (summarize . NE.fromList) .  DD.group . concatMap step  --------------------------------------------------------------------------------------- -- Combine@@ -109,7 +112,9 @@ class Combine a where     combine  :: a -> a -> a     summarize :: NE.NonEmpty a -> a+    combine a b = summarize . NE.fromList $ [a,b]     summarize (a NE.:| b) = foldr combine a b+    {-# MINIMAL combine | summarize #-}  instance (Combine a, Combine b) => Combine (a,b) where     combine (w,x) (y,z) = (combine w y, combine x z)@@ -123,10 +128,7 @@  -- |Easier way to write nested 2-tuples. type a :* b = (a,b)--- |Easier way to write nested 2-tuples,--- since @a >*\< b >*\< c >*< d@--- is much easier to read than--- @(((a,b),c),d)@.+-- |Easier way to write nested 2-tuples. -- Left associative, binds weaker than @+@ -- but stronger than @==@. (>*<) :: a -> b -> a :* b@@ -139,35 +141,33 @@  -- Does not group to combine unless equal. -- |Values from a 'Monoid' which have the respective--- binary operation applied each step.+-- binary operation applied each step,+-- where different values mean states should not be combined. -- E.g., strings with concatenation. newtype Merge a = Merge a     deriving (Eq, Generic)     deriving newtype (Semigroup, Monoid, Enum, Num, Fractional, Show)     deriving anyclass Grouping -instance Combine (Merge a) where-    combine = const+instance Combine (Merge a) where combine = const  --------------------------------------------------------------------------------------- -- Sum --------------------------------------------------------------------------------------- --- |Values which are added each step.+-- |Values which are added each step+-- where different values mean states should not be combined. -- E.g., number of times a red ball is picked from an urn. newtype Sum a = Sum a     deriving Generic     deriving newtype (Eq, Enum, Num, Fractional, Show)     deriving anyclass Grouping -instance Combine (Sum a) where-    combine = const+instance Combine (Sum a) where combine = const -instance Num a => Semigroup (Sum a) where-    x <> y = x + y+instance Num a => Semigroup (Sum a) where x <> y = x + y -instance Num a => Monoid (Sum a) where-    mempty = 0+instance Num a => Monoid (Sum a) where mempty = 0  --------------------------------------------------------------------------------------- -- Product@@ -186,17 +186,14 @@     grouping = FC.contramap (const ()) grouping  -- This causes Data.List.group to act more like Data.Discrimination.group-instance Eq (Product a) where-    _ == _ = True+-- |WARNING! Defined @_ == _ = True@!+instance Eq (Product a) where _ == _ = True -instance Num a => Combine (Product a) where-    combine = (+)+instance Num a => Combine (Product a) where combine = (+) -instance Num a => Semigroup (Product a) where-    x <> y = x * y+instance Num a => Semigroup (Product a) where x <> y = x * y -instance Num a => Monoid (Product a) where-    mempty = 1+instance Num a => Monoid (Product a) where mempty = 1  --------------------------------------------------------------------------------------- -- Misc@@ -204,16 +201,17 @@  -- |Randomly choose from a list by probability. randomProduct :: (Real a, MR.MonadRandom m) => [(a, b)] -> m (a, b)-randomProduct xs = MR.fromList . map (\x -> (x, toRational $ fst x)) $ xs+randomProduct = MR.fromList . fmap (\x -> (x, toRational $ fst x))  -- |Returns a single realization of a Markov chain.-randomPath :: (Markov a b, Real a, MR.RandomGen g) => (a,b) -> g -> [(a,b)]-randomPath x g = map (flip MR.evalRand g) . iterate (>>= (randomProduct . step)) $ pure x+randomPath :: (Markov ((,) a) b, Real a, MR.RandomGen g) => (a,b) -> g -> [(a,b)]+randomPath x g = fmap (`MR.evalRand` g) . iterate (>>= (randomProduct . step)) $ pure x  -- |Create a transition function from a transition matrix.--- If [[a]] is an n x n matrix, length [b] should be n.-fromLists :: Eq  b => [[a]] -> [b] -> b -> [(a, b -> b)]+--+-- prop> all (== length matrix) (map length matrix)+-- prop> length matrix == length states+fromLists :: Eq  b => [[a]] -> [b] -> b -> [(a, c -> b)] fromLists matrix states b = case DL.elemIndex b states of     Nothing -> []-    Just n  -> zip (matrix!!n) toState-    where toState = map const states+    Just n  -> zip (matrix!!n) $ fmap const states

@@ -0,0 +1,417 @@+{-# LANGUAGE DeriveGeneric              #-}+{-# LANGUAGE DeriveAnyClass             #-}+{-# LANGUAGE DerivingStrategies         #-}+{-# LANGUAGE FlexibleContexts           #-}+{-# LANGUAGE FlexibleInstances          #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE LambdaCase                 #-}+{-# LANGUAGE MultiParamTypeClasses      #-}+{-# LANGUAGE TypeOperators              #-}++{-|+Module      : Markov.Example+Description : Examples of Markov chains implemented using "Markov".+Maintainer  : atloomis@math.arizona.edu+Stability   : experimental++Several examples of Markov chains.+It is probably more helpful to read the source code than the Haddock documentation.+-}+module Markov.Example+     ( FromMatrix (..)+     , Simple (..)+     , Urn (..)+     , Extinction (..)+     , Tidal (..)+     , Room (..)+     , FillBin+     , initial+     , expectedLoss+     ) where++import Markov+import Generics.Deriving (Generic)+import Data.Discrimination (Grouping)++---------------------------------------------------------------+-- From a matrix+---------------------------------------------------------------++-- |An example defined from a matrix.+--+-- >>> chain [pure (FromMatrix 't') :: (Product Double, FromMatrix)] !! 100+-- [ (0.5060975609756099,'a')+-- , (0.201219512195122,'t')+-- , (0.29268292682926833,'l') ]+newtype FromMatrix = FromMatrix Char+    deriving Generic+    deriving newtype (Eq, Show)+    deriving anyclass Grouping++instance Combine FromMatrix where combine = const++instance Markov ((,) (Product Double)) FromMatrix where+    transition = let mat = [ [0.4, 0.3, 0.3]+                           , [0.2, 0.1, 0.7]+                           , [0.9, 0.1, 0.0] ]+                     chars = map FromMatrix ['a','t','l']+                 in fromLists mat chars++---------------------------------------------------------------+-- Simple random walk+---------------------------------------------------------------++-- |A simple random walk.+-- Possible outcomes of the first three steps:+--+-- >>> take 3 $ chain0 [Simple 0]+-- [ [0]+-- , [-1,1]+-- , [-2,0,2] ]+--+-- Probability of each outcome:+--+-- >>> take 3 $ chain [pure 0 :: (Product Double, Simple)]+-- [ [(1.0,0)]+-- , [(0.5,-1),(0.5,1)]+-- , [(0.25,-2),(0.5,0),(0.25,2)] ]+--+-- Number of ways to achieve each outcome:+--+-- >>> take 3 $ chain [pure 0 :: (Product Int, Simple)]+-- [ [(1,0)]+-- , [(1,-1),(1,1)]+-- , [(1,-2),(2,0),(1,2)] ]+--+-- Number of times @pred@ was applied,+-- allowing steps in place (@id@)+-- for more interesting output:+--+-- >>> chain [pure 0 :: (Sum Int, Simple)] !! 2+-- [ (2,-2), (1,-1), (1,0), (0,0), (0,1), (0,2) ]++newtype Simple = Simple Int+    deriving Generic+    deriving newtype (Num, Enum, Eq, Ord, Show)+    deriving anyclass Grouping++instance Combine Simple where combine = const++instance Markov0 Simple where+    transition0 _ = [pred, succ]++instance Markov ((,) (Product Double)) Simple where+    transition _ = [ 0.5 >*< pred+                   , 0.5 >*< succ ]++instance Markov ((,) (Product Int)) Simple where+    transition _ = [ 1 >*< pred+                   , 1 >*< succ ]++instance Markov ((,) (Sum Int)) Simple where+    transition _ = [ 1 >*< pred+                   , 0 >*< id+                   , 0 >*< succ ]++---------------------------------------------------------------+-- Urn model+---------------------------------------------------------------++-- |An urn contains balls of two colors.+-- At each step, a ball is chosen uniformly at random from the urn+-- and a ball of the same color is added.+newtype Urn = Urn (Int,Int)+    deriving Generic+    deriving newtype (Eq, Ord, Show)+    deriving anyclass Grouping++instance Combine Urn where combine = const++instance Markov ((,) (Product Double)) Urn where+    transition x = [ probLeft x >*< addLeft+                   , 1 - probLeft x >*< addRight ]++addLeft :: Urn -> Urn+addLeft  (Urn (a,b)) = Urn (a+1,b)++addRight :: Urn -> Urn+addRight (Urn (a,b)) = Urn (a,b+1)++probLeft :: Fractional a => Urn -> a+probLeft (Urn (a,b)) =  fromIntegral a / fromIntegral (a + b)++---------------------------------------------------------------+-- Tutorial+---------------------------------------------------------------++-- |This is the chain from the README.+newtype Extinction = Extinction Int+    deriving Generic+    deriving newtype (Eq, Num, Show)+    deriving anyclass Grouping++instance Combine Extinction where combine = const++instance Markov ((,) (Sum Int, Product Rational)) Extinction where+    transition = \case+        0 -> [ 0 >*< (q+r) >*< id+             , 0 >*< s >*< (+) 1 ]+        _ -> [ 1 >*< q >*< const 0+             , 0 >*< r >*< id+             , 0 >*< s >*< (+) 1 ]+      where q = 0.1; r = 0.3; s = 0.6++---------------------------------------------------------------+-- More complex random walk+---------------------------------------------------------------++-- |A time inhomogenous random walk that vaguely models tides+-- by periodically switching directions+-- and falling back from a shore at the origin.+data Tidal = Tidal { time     :: Double+                   , position :: Int }+                   deriving (Eq, Ord, Show, Generic)+                   deriving anyclass Grouping++instance Combine Tidal where combine = const++instance Markov ((,) (Product Double)) Tidal where+    transition tw = [ probRight tw >*< stepPos (+1)+                    , 1 - probRight tw >*< stepPos (flip (-) 1) ]++stepPos :: (Int -> Int) -> Tidal -> Tidal+stepPos f tw = Tidal (time tw + 1) (f $ position tw)++probRight :: Tidal -> Product Double+probRight tw = Product $ timeBias * positionBias+  where timeBias = (1 + sin (2 * pi * time tw / stepsPerCycle))/2+        positionBias+            | position tw >= 0 = 1 / steepness+            | otherwise       = 1+        stepsPerCycle = 10+        steepness     = 1.3 -- Double from 1 (flat) to +infty++---------------------------------------------------------------+-- Hidden Markov Model+---------------------------------------------------------------++-- |A hidden Markov model.+--+-- >>> :{ filter (\((_,Merge xs),_) -> xs == "aaa") $ chain+--  [1 >*< Merge "" >*< 1 :: Product Rational :* Merge String :* Room] !! 3+-- :}+-- [ ((3243 % 200000,"aaa"),Room 1)+-- , ((9729 % 500000,"aaa"),Room 2)+-- , ((4501 % 250000,"aaa"),Room 3) ]+--+-- Given that all three tokens recieved were @"a"@,+-- there is a probability of approximately @0.34@+-- that the current room is @Room 3@.+newtype Room = Room Int+    deriving (Generic, Show)+    deriving newtype (Eq, Num)+    deriving anyclass Grouping++instance Combine Room where combine = const++-- Note that changeState is applied before giveToken.+-- In spirit, we have @transition = giveToken . changeState@+instance Markov ((,) (Product Rational, Merge String)) Room where+    sequential = [giveToken, changeState]+      where changeState = \case+                1 -> [ 0.3 >*< mempty >*< const 1+                     , 0.6 >*< mempty >*< const 2+                     , 0.1 >*< mempty >*< const 3 ]+                2 -> [ 1.0 >*< mempty >*< const 3 ]+                3 -> [ 0.3 >*< mempty >*< const 1+                     , 0.6 >*< mempty >*< const 2+                     , 0.1 >*< mempty >*< const 3 ]+                _ -> error "State out of bounds in transition"+            giveToken = \case+                1 -> [ 0.5 >*< Merge "a" >*< const 1+                     , 0.5 >*< Merge "b" >*< const 1 ]+                2 -> [ 0.3 >*< Merge "a" >*< const 2+                     , 0.7 >*< Merge "b" >*< const 2 ]+                3 -> [ 0.4 >*< Merge "a" >*< const 3+                     , 0.4 >*< Merge "b" >*< const 3+                     , 0.2 >*< Merge "c" >*< const 3 ]+                _ -> error "State out of bounds in transition"++---------------------------------------------------------------+-- Yet more complex example+---------------------------------------------------------------++-- |Represents bins with free slots and items.+type Bin   = (Open,Full)+type Index = Int+-- |Represents space between bins where they can expand.+type Gap   = Int+type Full  = Int+type Open  = Int+type Trans = FillBin -> FillBin++-- |A collection of bins with gaps between them.+-- At each step an empty space is chosen+-- form a bin or from a gap.+-- If it is in a bin, the space is filled.+-- If it is in a gap, it is assigned to an adjacent bin,+-- which expands to contain it and any intervening spaces,+-- and then the space filled.+data FillBin = End Gap | Ext Gap Bin FillBin deriving (Eq, Ord, Generic, Grouping)++instance Show FillBin where+    show (Ext g b s) = show g ++ " " ++ show b ++ " " ++ show s+    show (End g) = show g++instance Combine FillBin where combine = const++instance Markov ((,) (Product Double)) FillBin where+    transition x = case probId x of+        0 -> filter (\(Product y,_) -> y /= 0) -- Careful, Product _ == Product _ = True+            $  [probAdd i x >*< addItem i | i <- indices]+            ++ [probGrowL i x >*< addItem i . growLeft  j i+                | i <- indices, j <- [1..gapN (i-1) x]]+            ++ [probGrowR i x >*< addItem i . growRight j i+                | i <- indices, j <- [1..gapN i x]]+        1 -> [pure id]+        _ -> error "Pattern not matched in transition"+        where indices = [1..size x]++-- |>>> fBFromLists [1,3,5,10] [(3,5),(9,9),(8,3)]+-- 1 (3,5) 3 (9,9) 5 (8,3) 10+fBFromLists :: [Gap] -> [Bin] -> FillBin+fBFromLists gaps bins = case (gaps,bins) of+    (g:_  , []  ) -> End g+    ([g]  , _   ) -> End g+    (g:gs , b:bs) -> Ext g b $ fBFromLists gs bs+    ([]   , _   ) -> End 0++-- |Create state where all bins start as (0,0).+--+-- >>> initial [5,7,0]+-- 5 (0,0) 7 (0,0) 0+initial :: [Int] -> FillBin+initial gs = fBFromLists gs $ repeat (0,0)++-- |The number of bins.+size :: FillBin -> Int+size = \case+    End _ -> 0+    Ext _ _ s -> 1 + size s++-- |The bins of a state.+getBins :: FillBin -> [Bin]+getBins = \case+    End _ -> []+    Ext _ b s -> b:getBins s++-- |The open values of a state.+getOpen :: FillBin -> [Open]+getOpen x = map fst $ getBins x++-- |The open value of the Nth bin.+openN :: Index -> FillBin -> Open+openN i x = getOpen x !!(i-1)++-- |The full values of a state.+getFull :: FillBin -> [Full]+getFull x = map snd $ getBins x++-- |The full value of the Nth bin.+fullN :: Index -> FillBin -> Full+fullN i x = getFull x !!(i-1)++-- |The gap values of a state.+getGap :: FillBin -> [Gap]+getGap = \case+    End g -> [g]+    Ext g _ s -> g:getGap s++-- |Warning! Indexed from zero!+gapN :: Index -> FillBin -> Gap+gapN i x = getGap x !! i++-- |The command @iApply i f s@ is analagous to+-- @take i s ++ f (drop i s)@.+iApply :: Trans -> Index -> Trans+iApply f idx x = case (idx,x) of+    (1, y) -> f y+    (i, Ext g b s) -> Ext g b $ iApply f (i-1) s+    _ -> error "Pattern not matched in iApply"++-- |Add an item to the ith bin.+addItem :: Index -> Trans+addItem = iApply h+    where h (Ext g (o,f) s) = Ext g (o-1,f+1) s+          h _ = error "pattern not matched in h in addItem"++-- |Expand the ith bin to the left by j.+-- The Markov chain will use @addItem i . growLeft j i@.+growLeft :: Int -> Index -> Trans+growLeft j = iApply h+    where h (Ext g (o,f) s) = Ext (g-j) (o+j,f) s+          h _ = error "pattern not matched in h in growLeft"++growRight :: Int -> Index -> Trans+growRight j = iApply h+    where h (Ext g (o,f) s) = Ext g (o+j,f) (shrink s)+          h _ = error "pattern not matched in h in growRight"+          shrink = \case+              End g -> End (g-j)+              Ext g b t -> Ext (g-j) b t++-- |The sum of all open slots in bins and gaps.+slots :: FillBin -> Int+slots x = sum $ getGap x ++ getOpen x++-- |The probability that a state returns to itself.+probId :: Num a => FillBin -> a+probId x+    | slots x == 0 = 1+    | otherwise    = 0++divInt :: (Integral a, Integral b, Fractional c) => a -> b -> c+divInt x y = fromIntegral x / fromIntegral y++-- |The probability that the ith bin gains an item.+probAdd :: Fractional a => Index -> FillBin -> a+probAdd i x = openN i x `divInt` slots x++-- |The probability that the ith bin expands to the left.+probGrowL :: Fractional a => Index -> FillBin -> a+probGrowL i x+    | test      = 1 `divInt` slots x+    | otherwise = 0+    where test = i == 1 || fullN i x < fullN (i-1) x++-- |The probability that the ith bin expands to the right.+probGrowR :: Fractional a => Index -> FillBin -> a+probGrowR i x+    | test      = 1 `divInt` slots x+    | otherwise = 0+    where test = i == size x || fullN i x <= fullN (i+1) x++---------------------------------------------------------------+-- Several functions to help study the previous process+---------------------------------------------------------------++-- |The \(l^2\) distance between a finished state+-- and a state with perfectly balanced bins.+individualLoss :: Fractional a => FillBin -> a+individualLoss x = sum . map f . getFull $ x+    where f y = (fromIntegral y - ideal)^2+          ideal = sum (getFull x) `divInt` size x++probLoss :: Fractional a => (Product a, FillBin) -> a+probLoss (Product x, y) = x * individualLoss y++-- |Expected loss of a set of states of @['FillBin']@.+-- Loss is the \(l^2\) distance between a finished state+-- and a state with perfectly balanced bins.+--+-- >>> expectedLoss [pure $ initial [1,0,3] :: (Product Double, FillBin)]+-- 2.0+expectedLoss :: (Fractional a, Markov ((,) (Product a)) FillBin) +    => [Product a :* FillBin] -> a+expectedLoss xs = sum . map probLoss $ chain xs !! idx+    where idx = slots . snd . head $ xs

@@ -1,422 +0,0 @@-{-# LANGUAGE MultiParamTypeClasses, FlexibleContexts, DeriveGeneric, DeriveAnyClass,-DerivingStrategies, GeneralizedNewtypeDeriving, FlexibleInstances, TypeOperators #-}-{-|-Module      : Examples-Description : Examples of Markov chains implemented using "Markov".-Maintainer  : atloomis@math.arizona.edu-Stability   : experimental--Several examples of Markov chains.-It is probably more helpful to read the source code than the Haddock documentation.--}-module Markov.Examples ( FromMatrix (..)-                       , Simple (..)-                       , Urn (..)-                       , Extinction (..)-                       , Tidal (..)-                       , Room (..)-                       , FillBin-                       , initial-                       , expectedLoss-                       ) where--import Markov-import Generics.Deriving (Generic)-import Data.Discrimination (Grouping)-------------------------------------------------------------------- From a matrix-------------------------------------------------------------------- |An example defined from a matrix.------ >>> chain [pure 't' :: Product Double :* Char] !! 100--- [ (0.5060975609756099,'a')--- , (0.201219512195122,'t')--- , (0.29268292682926833,'l') ]-newtype FromMatrix = FromMatrix Char-    deriving Generic-    deriving newtype (Eq, Show)-    deriving anyclass Grouping--instance Markov (Product Double) FromMatrix where-    transition = let mat = [ [0.4, 0.3, 0.3]-                           , [0.2, 0.1, 0.7]-                           , [0.9, 0.1, 0.0] ]-                     chars = map FromMatrix ['a','t','l']-                 in fromLists mat chars-    ------------------------------------------------------------------- Simple random walk-------------------------------------------------------------------- |A simple random walk.--- Possible outcomes of the first three steps:------ >>> take 3 $ chain0 [Simple 0]--- [ [0]--- , [-1,1]--- , [-2,0,2]]------ Probability of each outcome:------ >>> take 3 $ chain [pure 0 :: Product Double :* Simple]--- [ [(1.0,0)]--- , [(0.5,-1),(0.5,1)]--- , [(0.25,-2),(0.5,0),(0.25,2)] ]------ Number of ways to achieve each outcome:------ >>> take 3 $ chain [pure 0 :: Product Int :* Simple]--- [ [(1,0)]--- , [(1,-1),(1,1)]--- , [(1,-2),(2,0),(1,2)] ]------ Number of times @pred@ was applied,--- allowing steps in place (@id@)--- for more interesting output:------ >>> chain [pure 0 :: Sum Int :* Simple] !! 2--- [ (2,-2)--- , (1,-1)--- , (1,0)--- , (0,0)--- , (0,1)--- , (0,2) ]--newtype Simple = Simple Int-    deriving Generic-    deriving newtype (Num, Enum, Eq, Ord, Show)-    deriving anyclass Grouping--instance Markov0 Simple where-    transition0 _ = [pred, succ]--instance Markov (Product Double) Simple where-    transition _ = [ 0.5 >*< pred-                   , 0.5 >*< succ ]--instance Markov (Product Int) Simple where-    transition _ = [ 1 >*< pred-                   , 1 >*< succ ]--instance Markov (Sum Int) Simple where-    transition _ = [ 1 >*< pred-                   , 0 >*< id-                   , 0 >*< succ ]-              -- = [ 1 >*< pred-              --   , pure id-              --   , pure succ ]-------------------------------------------------------------------- Urn model-------------------------------------------------------------------- |An urn contains balls of two colors.--- At each step, a ball is chosen uniformly at random from the urn--- and a ball of the same color is added.-newtype Urn = Urn (Int,Int)-    deriving Generic-    deriving newtype (Eq, Ord, Show)-    deriving anyclass Grouping--instance Markov (Product Double) Urn where-    transition x = [ probLeft x >*< addLeft-                   , 1 - probLeft x >*< addRight ]--addLeft :: Urn -> Urn-addLeft  (Urn (a,b)) = Urn (a+1,b)--addRight :: Urn -> Urn-addRight (Urn (a,b)) = Urn (a,b+1)--probLeft :: Fractional a => Urn -> a-probLeft (Urn (a,b)) = (fromIntegral a)/(fromIntegral $ a + b)-------------------------------------------------------------------- Tutorial-------------------------------------------------------------------- |This is the chain from the README.-newtype Extinction = Extinction Int-    deriving Generic-    deriving newtype (Eq, Num, Show)-    deriving anyclass Grouping--instance Markov (Sum Int, Product Rational) Extinction where-    transition x = case x of-        0 -> [ 0 >*< (q+r) >*< id-             , 0 >*< s >*< (+1) ]-        _ -> [ 1 >*< q >*< const 0-             , 0 >*< r >*< id-             , 0 >*< s >*< (+1) ]-        where q = 0.1; r = 0.3; s = 0.6---- This is equivalent to the definition above.-instance Combine Extinction where-    combine = const--instance Semigroup Extinction where-    (<>) = flip const--instance MultiMarkov (Sum Int :* Product Rational :* Extinction) where-    multiTransition _ = [trans]-        where trans ((_,_),z) = case z of-                  0 -> [ 0 >*< (q+r) >*< 0-                       , 0 >*< s >*< 1 ]-                  x -> [ 1 >*< q >*< 0-                       , 0 >*< r >*< x-                       , 0 >*< s >*< x+1 ]-                  where q = 0.1; r = 0.3; s = 0.6-------------------------------------------------------------------- More complex random walk-------------------------------------------------------------------- |A time inhomogenous random walk that vaguely models tides--- by periodically switching directions--- and falling back from a shore at the origin.-data Tidal = Tidal { time     :: Double-                   , position :: Int }-                   deriving (Eq, Ord, Show, Generic)-                   deriving anyclass Grouping--instance Markov (Product Double) Tidal where-    transition tw = [ probRight tw >*< stepPos (+1)-                    , 1 - (probRight tw) >*< stepPos (flip (-) 1) ]--stepPos :: (Int -> Int) -> Tidal -> Tidal-stepPos f tw = Tidal (time tw + 1) (f $ position tw)--probRight :: Tidal -> Product Double-probRight tw = Product $ timeBias * positionBias-    where timeBias = (1 + sin (2 * pi * (time tw) / stepsPerCycle))/2-          positionBias-              | position tw >= 0 = 1 / steepness-              | otherwise       = 1-          stepsPerCycle = 10-          steepness     = 1.3 -- Double from 1 (flat) to +infty-------------------------------------------------------------------- Hidden Markov Model-------------------------------------------------------------------- |A hidden Markov model.------ >>> filter (\((_,Merge xs),_) -> xs == "aaa") $ multiChain [1 >*< Merge "" >*< 1 :: Product Rational :* Merge String :* Room] !! 3--- [ ((3243 % 200000,"aaa"),Room 1)--- , ((9729 % 500000,"aaa"),Room 2)--- , ((4501 % 250000,"aaa"),Room 3) ]------ Given that all three tokens recieved were @"a"@,--- there is a probability of approximately @0.34@--- that the current room is @Room 3@.-newtype Room = Room Int-    deriving (Generic, Show)-    deriving newtype (Eq, Num)-    deriving anyclass Grouping--instance Semigroup Room where-    (<>) = flip const--instance Combine Room where-    combine = const---- Note that changeState is applied before giveToken.--- In spirit, we have stepj = giveToken . changeState-instance MultiMarkov (Product Rational :* Merge String :* Room) where-    multiTransition _ = [giveToken, changeState]-        where changeState ((_,_),z) = case z of-                  1 -> [ 0.3 >*< mempty >*< 1-                       , 0.6 >*< mempty >*< 2-                       , 0.1 >*< mempty >*< 3 ]-                  2 -> [ 1.0 >*< mempty >*< 3 ]-                  3 -> [ 0.3 >*< mempty >*< 1-                       , 0.6 >*< mempty >*< 2-                       , 0.1 >*< mempty >*< 3 ]-                  _ -> error "State out of bounds in transitionk"-              giveToken ((_,_),z) = case z of-                  1 -> [ 0.5 >*< Merge "a" >*< 1-                       , 0.5 >*< Merge "b" >*< 1 ]-                  2 -> [ 0.3 >*< Merge "a" >*< 2-                       , 0.7 >*< Merge "b" >*< 2 ]-                  3 -> [ 0.4 >*< Merge "a" >*< 3-                       , 0.4 >*< Merge "b" >*< 3-                       , 0.2 >*< Merge "c" >*< 3 ]-                  _ -> error "State out of bounds in transitionk"-------------------------------------------------------------------- Yet more complex example-------------------------------------------------------------------- |Represents bins with free slots and items.-type Bin   = (Open,Full)-type Index = Int--- |Represents space between bins where they can expand.-type Gap   = Int-type Full  = Int-type Open  = Int-type Trans = FillBin -> FillBin---- |A collection of bins with gaps between them.--- At each step an empty space is chosen--- form a bin or from a gap.--- If it is in a bin, the space is filled.--- If it is in a gap, it is assigned to an adjacent bin,--- which expands to contain it and any intervening spaces,--- and then the space filled.-data FillBin = End Gap | Ext Gap Bin FillBin deriving (Eq, Ord, Generic, Grouping)--instance Show FillBin where-    show (Ext g b s) = show g ++ " " ++ show b ++ " " ++ show s-    show (End g) = show g--instance Markov (Product Double) FillBin where-    transition x = case probId x of-        0 -> filter (\(Product y,_) -> y /= 0) -- Careful, Product _ == Product _ = True-            $  [probAdd i x >*< addItem i | i <- indices]-            ++ [probGrowL i x >*< addItem i . growLeft  j i-                | i <- indices, j <- [1..gapN (i-1) x]]-            ++ [probGrowR i x >*< addItem i . growRight j i-                | i <- indices, j <- [1..gapN i x]]-        1 -> [pure id]-        _ -> error "Pattern not matched in transition"-        where indices = [1..size x]---- |>>> fBFromLists [1,3,5,10] [(3,5),(9,9),(8,3)]--- 1 (3,5) 3 (9,9) 5 (8,3) 10-fBFromLists :: [Gap] -> [Bin] -> FillBin-fBFromLists gaps bins = case (gaps,bins) of-    (g:_  , []  ) -> End g-    ([g]  , _   ) -> End g-    (g:gs , b:bs) -> Ext g b $ fBFromLists gs bs-    ([]   , _   ) -> End 0---- |Create state where all bins start as (0,0).------ >>> initial [5,7,0]--- 5 (0,0) 7 (0,0) 0-initial :: [Int] -> FillBin-initial gs = fBFromLists gs $ repeat (0,0)---- |The number of bins.-size :: FillBin -> Int-size x = case x of-    End _ -> 0-    Ext _ _ s -> 1 + size s---- |The bins of a state.-getBins :: FillBin -> [Bin]-getBins x = case x of-    End _ -> []-    Ext _ b s -> b:getBins s---- |The open values of a state.-getOpen :: FillBin -> [Open]-getOpen x = map fst $ getBins x---- |The open value of the Nth bin.-openN :: Index -> FillBin -> Open-openN i x = (getOpen x)!!(i-1)---- |The full values of a state.-getFull :: FillBin -> [Full]-getFull x = map snd $ getBins x---- |The full value of the Nth bin.-fullN :: Index -> FillBin -> Full-fullN i x = (getFull x)!!(i-1)---- |The gap values of a state.-getGap :: FillBin -> [Gap]-getGap x = case x of-    End g -> [g]-    Ext g _ s -> g:getGap s---- |Warning! Indexed from zero!-gapN :: Index -> FillBin -> Gap-gapN i x = (getGap x)!!i---- |The command @iApply i f s@ is analagous to--- @take i s ++ f (drop i s)@.-iApply :: Trans -> Index -> Trans-iApply f idx x = case (idx,x) of-    (1, y) -> f y-    (i, Ext g b s) -> Ext g b $ iApply f (i-1) s-    _ -> error "Pattern not matched in iApply"---- |Add an item to the ith bin.-addItem :: Index -> Trans-addItem = iApply h-    where h (Ext g (o,f) s) = Ext g (o-1,f+1) s-          h _ = error "pattern not matched in h in addItem"---- |Expand the ith bin to the left by j.--- The Markov chain will use @addItem i . growLeft j i@.-growLeft :: Int -> Index -> Trans-growLeft j = iApply h-    where h (Ext g (o,f) s) = Ext (g-j) (o+j,f) s-          h _ = error "pattern not matched in h in growLeft"--growRight :: Int -> Index -> Trans-growRight j = iApply h-    where h (Ext g (o,f) s) = Ext g (o+j,f) (shrink s)-          h _ = error "pattern not matched in h in growRight"-          shrink s = case s of-              End g -> End (g-j)-              Ext g b t -> Ext (g-j) b t---- |The sum of all open slots in bins and gaps.-slots :: FillBin -> Int-slots x = sum $ getGap x ++ getOpen x---- |The probability that a state returns to itself.-probId :: Num a => FillBin -> a-probId x = case slots x == 0 of-    True  -> 1-    False -> 0--divInt :: (Integral a, Integral b, Fractional c) => a -> b -> c-divInt x y = (fromIntegral x)/(fromIntegral y)---- |The probability that the ith bin gains an item.-probAdd :: Fractional a => Index -> FillBin -> a-probAdd i x = openN i x `divInt` slots x---- |The probability that the ith bin expands to the left.-probGrowL :: Fractional a => Index -> FillBin -> a-probGrowL i x = case test of-    True  -> 1 `divInt` slots x-    False -> 0-    where test = i == 1 || fullN i x < fullN (i-1) x---- |The probability that the ith bin expands to the right.-probGrowR :: Fractional a => Index -> FillBin -> a-probGrowR i x = case test of-    True  -> 1 `divInt` slots x-    False -> 0-    where test = i == size x || fullN i x <= fullN (i+1) x-------------------------------------------------------------------- Several functions to help study the previous process-------------------------------------------------------------------- |The \(l^2\) distance between a finished state--- and a state with perfectly balanced bins.-individualLoss :: Fractional a => FillBin -> a-individualLoss x = sum . map f . getFull $ x-    where f y = (fromIntegral y - ideal)^2-          ideal = sum (getFull x) `divInt` size x--probLoss :: Fractional a => (Product a, FillBin) -> a-probLoss (Product x, y) = x * individualLoss y---- |Expected loss of a set of pstates of @['FillBin']@.--- Loss is the \(l^2\) distance between a finished state--- and a state with perfectly balanced bins.------ >>> expectedLoss [pure $ initial [1,0,3] :: Product Double :* FillBin]--- 2.0-expectedLoss :: (Fractional a, Markov (Product a) FillBin) => [Product a :* FillBin] -> a-expectedLoss xs = sum . map probLoss $ (chain xs) !! idx-    where idx = slots . snd . head $ xs

@@ -0,0 +1,8 @@+module Markov.Instance where++import Data.Discrimination (Grouping, grouping)+import Data.Functor.Contravariant (contramap)+import GHC.Float (castFloatToWord32, castDoubleToWord64)++instance Grouping Float where grouping = contramap castFloatToWord32 grouping+instance Grouping Double where grouping = contramap castDoubleToWord64 grouping

@@ -1,9 +0,0 @@-module Markov.Instances where--import Data.Discrimination (Grouping, grouping)-import Data.Functor.Contravariant (contramap)-import Data.Functor.Contravariant.Divisible (conquer)-import GHC.Float (castFloatToWord32, castDoubleToWord64)--instance Grouping Float where grouping = contramap castFloatToWord32 grouping-instance Grouping Double where grouping = contramap castDoubleToWord64 grouping

@@ -0,0 +1,57 @@+{-# OPTIONS_GHC -F -pgmF htfpp #-}+{-# LANGUAGE TypeOperators     #-}++import Data.Ratio+import Test.Framework++import Markov+import Markov.Example++main = htfMain htf_thisModulesTests++-- Examples in the documentation.+test_fromMatrix =+    assertEqual+    (chain [pure (FromMatrix 't') :: (Product Double, FromMatrix)] !! 100)+    [ (0.5060975609756099, FromMatrix 'a')+    , (0.201219512195122, FromMatrix 't')+    , (0.29268292682926833, FromMatrix 'l') ]++test_m0Simple =+    assertEqual+    (take 3 $ chain0 [Simple 0])+    [ [0]+    , [-1,1]+    , [-2,0,2] ]++test_pdSimple =+    assertEqual+    (take 3 $ chain [pure 0 :: (Product Double, Simple)])+    [ [(1.0,0)]+    , [(0.5,-1),(0.5,1)]+    , [(0.25,-2),(0.5,0),(0.25,2)] ]++test_piSimple =+    assertEqual+    (take 3 $ chain [pure 0 :: (Product Int, Simple)])+    [ [(1,0)]+    , [(1,-1),(1,1)]+    , [(1,-2),(2,0),(1,2)] ]++test_siSimple =+    assertEqual+    (chain [pure 0 :: (Sum Int, Simple)] !! 2)+    [ (2,-2), (1,-1), (1,0), (0,0), (0,1), (0,2) ]++test_HMM =+    assertEqual+    (filter (\((_,Merge xs),_) -> xs == "aaa") $ chain+    [1 >*< Merge "" >*< 1 :: Product Rational :* Merge String :* Room] !! 3)+    [ ((Product $ 3243 % 200000, Merge "aaa"),Room 1)+    , ((Product $ 9729 % 500000, Merge "aaa"),Room 2)+    , ((Product $ 4501 % 250000, Merge "aaa"),Room 3) ]++test_expLoss =+    assertEqual+    (expectedLoss [pure $ initial [1,0,3] :: (Product Double, FillBin)])+    2