markov-realization (empty) → 0.1.0
raw patch · 8 files changed
+789/−0 lines, 8 filesdep +MonadRandomdep +basedep +contravariantsetup-changed
Dependencies added: MonadRandom, base, contravariant, discrimination, generic-deriving
Files
- ChangeLog.md +3/−0
- LICENSE +30/−0
- README.md +66/−0
- Setup.hs +2/−0
- markov-realization.cabal +38/−0
- src/Markov.hs +219/−0
- src/Markov/Examples.hs +422/−0
- src/Markov/Instances.hs +9/−0
+ ChangeLog.md view
@@ -0,0 +1,3 @@+# Changelog for markov++## Unreleased changes
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright Author name here (c) 2019++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 Author name here 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.
+ README.md view
@@ -0,0 +1,66 @@+# Markov Tutorial++Let X<sub>n</sub> denote the nth state of a Markov chain with state space ℕ.+For x ≠ 0 define transition probabilities++p(x,0) = q,++p(x,x) = r, and++p(x,x+1) = s.++When x = 0, let+p(x,0) = q+r,+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],+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.+Extinctions add a value to a counter each time they happen+and the counter takes integral values,+so they can be represented by `Sum Int`.+Probabilities are multiplied each step,+and added when duplicate steps are combined.+We want decimal probabilities, so+we can represent this with `Product Rational`.+We will make a new type for the state.++```haskell+data Extinction = Extinction Int+ deriving Generic+ deriving newtype (Eq, Num, Show)+ deriving anyclass Grouping+```++All that remains is to make an instance of `Markov`.++```haskell+instance Markov (Sum Int, Product Rational) Extinction where+ transition x = case state 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+```++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>++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.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ markov-realization.cabal view
@@ -0,0 +1,38 @@+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+extra-source-files:+ README.md+ ChangeLog.md++source-repository head+ 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
+ src/Markov.hs view
@@ -0,0 +1,219 @@+{-# LANGUAGE MultiParamTypeClasses, GeneralizedNewtypeDeriving, DeriveGeneric,+DeriveAnyClass, DerivingStrategies, TypeOperators #-}+{-|+Module : Markov+Description : Realization of Markov processes with known parameters.+Maintainer : atloomis@math.arizona.edu+Stability : experimental++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.+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++import Markov.Instances ()+import Control.Applicative ((<**>))+import Generics.Deriving (Generic)+import Data.Discrimination (Grouping, grouping)+import qualified Data.Discrimination as DD+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+---------------------------------------------------------------++-- |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++---------------------------------------------------------------------------------------+-- 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+---------------------------------------------------------------------------------------++-- |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++---------------------------------------------------------------------------------------+-- Combine+---------------------------------------------------------------------------------------++-- |Within equivalence classes, @combine@ should be associative,+-- commutative, and should be idempotent up to equivalence.+-- I.e. if @x == y == z@,+--+-- prop> (x `combine` y) `combine` z = x `combine` (y `combine` z)+-- prop> x `combine` y = y `combine` x+-- prop> x `combine` x == x+class Combine a where+ combine :: a -> a -> a+ summarize :: NE.NonEmpty a -> a+ summarize (a NE.:| b) = foldr combine a b++instance (Combine a, Combine b) => Combine (a,b) where+ combine (w,x) (y,z) = (combine w y, combine x z)++instance (Combine a, Combine b, Combine c) => Combine (a,b,c) where+ combine (a,w,x) (b,y,z) = (combine a b, combine w y, combine x z)++---------------------------------------------------------------------------------------+-- Easier way to write nested 2-tuples+---------------------------------------------------------------------------------------++-- |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)@.+-- Left associative, binds weaker than @+@+-- but stronger than @==@.+(>*<) :: a -> b -> a :* b+a >*< b = (a,b)+infixl 5 >*<++---------------------------------------------------------------------------------------+-- Merge+---------------------------------------------------------------------------------------++-- Does not group to combine unless equal.+-- |Values from a 'Monoid' which have the respective+-- binary operation applied each step.+-- 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++---------------------------------------------------------------------------------------+-- Sum+---------------------------------------------------------------------------------------++-- |Values which are added each step.+-- 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 Num a => Semigroup (Sum a) where+ x <> y = x + y++instance Num a => Monoid (Sum a) where+ mempty = 0++---------------------------------------------------------------------------------------+-- Product+---------------------------------------------------------------------------------------++-- Does not effect equality of tuple,+-- @combine x y = x + y@.+-- |Values which are multiplied each step,+-- and combined additively for equal states.+-- E.g., probabilities.+newtype Product a = Product a+ deriving Generic+ deriving newtype (Num, Fractional, Enum, Show)++instance Grouping (Product a) where+ grouping = FC.contramap (const ()) grouping++-- This causes Data.List.group to act more like Data.Discrimination.group+instance Eq (Product a) where+ _ == _ = True++instance Num a => Combine (Product a) where+ combine = (+)++instance Num a => Semigroup (Product a) where+ x <> y = x * y++instance Num a => Monoid (Product a) where+ mempty = 1++---------------------------------------------------------------------------------------+-- Misc+---------------------------------------------------------------------------------------++-- |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++-- |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++-- |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)]+fromLists matrix states b = case DL.elemIndex b states of+ Nothing -> []+ Just n -> zip (matrix!!n) toState+ where toState = map const states
+ src/Markov/Examples.hs view
@@ -0,0 +1,422 @@+{-# 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
+ src/Markov/Instances.hs view
@@ -0,0 +1,9 @@+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