diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -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ₙ = j ∩ d = k], +P\[Xₙ = j and 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. @@ -52,7 +52,7 @@ where q = 0.1; r = 0.3; s = 0.6 ``` -We can now easily see a list of states, deaths, and the probabilities. +We can now see a list of states, extinctions, and the probabilities. __`> chain [pure 0 :: Sum Int :* Product Rational :* Extinction] !! 3`__ ``` @@ -66,4 +66,5 @@ This means that starting from a state of zero, after three time steps there is a 51/500 chance -that the state is zero and there has been one death. +that the state is zero and there has been one extinction. + diff --git a/markov-realization.cabal b/markov-realization.cabal --- a/markov-realization.cabal +++ b/markov-realization.cabal @@ -1,19 +1,13 @@ cabal-version: 1.12 --- 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.3.3 +version: 0.4 description: Please see the README on GitHub at 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 +copyright: 2020 Alex Loomis license: BSD3 license-file: LICENSE build-type: Simple @@ -55,6 +49,7 @@ -Wall -Wcompat -Wincomplete-uni-patterns -Wincomplete-record-updates -Wredundant-constraints build-depends: base >=4.7 && <5 - , HTF , markov-realization + , HTF + , MonadRandom default-language: Haskell2010 diff --git a/src/Markov.hs b/src/Markov.hs --- a/src/Markov.hs +++ b/src/Markov.hs @@ -140,8 +140,10 @@ -- and combined additively for equal states. -- E.g., probabilities. newtype Product a = Product a - deriving newtype (Num, Fractional, Enum, Show) + deriving newtype (Enum, Fractional, Integral, Num, Real, Show) +-- Don't sort by probability. +-- |WARNING! Defined @compare _ _ = EQ@! instance Ord (Product a) where compare _ _ = EQ diff --git a/src/Markov/Example.hs b/src/Markov/Example.hs --- a/src/Markov/Example.hs +++ b/src/Markov/Example.hs @@ -115,6 +115,9 @@ -- |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. +-- +-- >>> randomPath (mkStdGen 70) (Urn (2,5)) !! 8 :: (Product Double, Urn) +-- (0.1648351648351649, Urn (2,13)) newtype Urn = Urn (Int,Int) deriving newtype (Eq, Ord, Show) diff --git a/src/Markov/Extra.hs b/src/Markov/Extra.hs --- a/src/Markov/Extra.hs +++ b/src/Markov/Extra.hs @@ -1,9 +1,10 @@ -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE TypeOperators #-} {-| Module : Markov.Extra Maintainer : atloomis@math.arizona.edu +Description : Utilities for "Markov". Stability : Experimental -} @@ -14,9 +15,9 @@ , (>*<) ) where -import Markov import qualified Control.Monad.Random as MR -import qualified Data.List as DL +import qualified Data.List as DL +import Markov --------------------------------------------------------------------------------------- -- Misc @@ -27,8 +28,9 @@ 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 = fmap (`MR.evalRand` g) . iterate (>>= (randomProduct . step)) $ pure x +randomPath :: (Markov ((,) a) b, Real a, MR.RandomGen g) => g -> b -> [(a,b)] +randomPath g x = fmap (`MR.evalRand` g) . iterate (>>= (randomProduct . step)) + $ pure (1,x) -- |Create a transition function from a transition matrix. -- Inputs should obey: @@ -52,3 +54,4 @@ (>*<) :: a -> b -> a :* b a >*< b = (a,b) infixl 5 >*< + diff --git a/test/Test.hs b/test/Test.hs --- a/test/Test.hs +++ b/test/Test.hs @@ -1,6 +1,7 @@ {-# OPTIONS_GHC -F -pgmF htfpp #-} -{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE TypeOperators #-} +import Control.Monad.Random (mkStdGen) import Data.Ratio import Test.Framework @@ -10,47 +11,46 @@ main = htfMain htf_thisModulesTests --- Examples in the documentation. test_fromLists = assertEqual - (chain [pure (FromLists 't') :: (Product Double, FromLists)] !! 100) [ (0.5060975609756099, FromLists 'a') , (0.29268292682926833, FromLists 'l') , (0.201219512195122, FromLists 't') ] + (chain [pure (FromLists 't') :: (Product Double, FromLists)] !! 100) test_m0Simple = assertEqual - (take 3 $ chain0 [Simple 0]) [ [0] , [-1,1] , [-2,0,2] ] + (take 3 $ chain0 [Simple 0]) 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)] ] + (take 3 $ chain [pure 0 :: (Product Double, Simple)]) test_piSimple = assertEqual - (take 3 $ chain [pure 0 :: (Product Int, Simple)]) [ [(1,0)] , [(1,-1),(1,1)] , [(1,-2),(2,0),(1,2)] ] + (take 3 $ chain [pure 0 :: (Product Int, Simple)]) test_siSimple = assertEqual - (chain [pure 0 :: (Sum Int, Simple)] !! 2) [ (0,0), (0,1), (0,2), (1,-1), (1,0), (2,-2) ] + (chain [pure 0 :: (Sum Int, Simple)] !! 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) ] + (filter (\((_,Merge xs),_) -> xs == "aaa") $ chain + [1 >*< Merge "" >*< 1 :: Product Rational :* Merge String :* Room] !! 3) test_expLoss = assertEqual @@ -59,5 +59,11 @@ test_expLossMore = assertEqual - (expectedLoss [pure $ initial [1,0,3,2,2,5] :: (Product Double, FillBin)]) 9.764425505050506 + (expectedLoss [pure $ initial [1,0,3,2,2,5] :: (Product Double, FillBin)]) + +test_randomPath = + assertEqual + (0.1648351648351649, Urn (2,13)) + (randomPath (mkStdGen 70) (Urn (2,5)) !! 8 :: (Product Double, Urn)) +