diff --git a/ToDo b/ToDo --- a/ToDo +++ b/ToDo @@ -1,6 +1,3 @@ -the library could benefit from Monad transformers like WriterT and Product, see - http://www.randomhacks.net/articles/2007/02/21/refactoring-probability-distributions - Examples: Election and prognoses diff --git a/probability.cabal b/probability.cabal --- a/probability.cabal +++ b/probability.cabal @@ -1,5 +1,5 @@ Name: probability -Version: 0.2.2.1 +Version: 0.2.3 License: BSD3 License-File: COPYRIGHT Author: Martin Erwig , Steve Kollmansberger @@ -13,7 +13,8 @@ The monad is similar to the List monad for non-deterministic computations, but extends the List monad by a measure of probability. Small interface to R plotting. -Tested-With: GHC==6.4.1, GHC==6.8.2 +Tested-With: GHC==6.4.1, GHC==6.8.2, GHC==6.12.3 +Tested-With: GHC==7.0.2, GHC==7.2.2 Cabal-Version: >=1.6 Build-Type: Simple Data-Files: @@ -27,14 +28,16 @@ Source-Repository this type: darcs location: http://code.haskell.org/~thielema/probability/ - tag: 0.2.2.1 + tag: 0.2.3 Flag splitBase description: Choose the new smaller, split-up base package. Library - Build-Depends: transformers >=0.0.1 && <0.3 + Build-Depends: + utility-ht >=0.0.6 && <0.1, + transformers >=0.0.1 && <0.3 If flag(splitBase) Build-Depends: containers >=0.1 && <1, @@ -42,6 +45,7 @@ base >=2 && <5 Else Build-Depends: + special-functors >=1.0 && <1.1, base >=1.0 && <2 Hs-Source-Dirs: src @@ -65,6 +69,7 @@ Numeric.Probability.Example.Diagnosis Numeric.Probability.Example.Dice Numeric.Probability.Example.DiceAccum + Numeric.Probability.Example.Kruskal Numeric.Probability.Example.MontyHall Numeric.Probability.Example.NBoys Numeric.Probability.Example.Predator @@ -75,3 +80,4 @@ Numeric.Probability.Monad Numeric.Probability.PrintList Numeric.Probability.Show + Numeric.Probability.Either diff --git a/src/Numeric/Probability/Distribution.hs b/src/Numeric/Probability/Distribution.hs --- a/src/Numeric/Probability/Distribution.hs +++ b/src/Numeric/Probability/Distribution.hs @@ -5,10 +5,15 @@ import Numeric.Probability.Show (showR) import qualified Numeric.Probability.Shape as Shape +import Control.Applicative (Applicative(..)) import Control.Monad (liftM, liftM2, join, ) +import qualified Data.List.HT as ListHT import qualified Data.Map as Map import qualified Data.List as List +import qualified Data.Maybe as Maybe +import Data.Ord.HT (comparing, ) +import Data.Eq.HT (equating, ) import Prelude hiding (map, filter) @@ -44,10 +49,14 @@ If you need efficiency, you should remove redundant elements by 'norm'. 'norm' converts to 'Data.Map' and back internally and you can hope that the compiler fuses the lists with the intermediate Map structure. + +The defined monad is equivalent to +@WriterT (Product prob) [] a@. +See . -} newtype T prob a = Cons {decons :: [(a,prob)]} -certainly :: Num prob => a -> T prob a +certainly :: Num prob => a -> T prob a certainly x = Cons [(x,1)] instance Num prob => Monad (T prob) where @@ -55,6 +64,10 @@ d >>= f = Cons [(y,q*p) | (x,p) <- decons d, (y,q) <- decons (f x)] fail _ = Cons [] +instance Num prob => Applicative (T prob) where + pure = certainly + fm <*> m = Cons [(f x,q*p) | (f,p) <- decons fm, (x,q) <- decons m] + {- Dist cannot be an instance of MonadPlus, because there is no mzero @@ -115,24 +128,24 @@ size :: T prob a -> Int size = length . decons -check :: RealFloat prob => T prob a -> T prob a +check :: (RealFloat prob, Show prob) => T prob a -> T prob a check (Cons d) = if abs (1-sumP d) < errorMargin then Cons d else error ("Illegal distribution: total probability = "++show (sumP d)) -- | can fail because of rounding errors, better use 'fromFreqs' -cons :: RealFloat prob => [(a,prob)] -> T prob a +cons :: (RealFloat prob, Show prob) => [(a,prob)] -> T prob a cons = check . Cons sumP :: Num prob => [(a,prob)] -> prob sumP = sum . List.map snd sortP :: Ord prob => [(a,prob)] -> [(a,prob)] -sortP = List.sortBy (\x y -> compare (snd y) (snd x)) +sortP = List.sortBy (comparing snd) sortElem :: Ord a => [(a,prob)] -> [(a,prob)] -sortElem = List.sortBy (\x y -> compare (fst y) (fst x)) +sortElem = List.sortBy (comparing fst) -- ** Normalization = grouping @@ -145,8 +158,11 @@ norm'' :: (Num prob, Ord a) => [(a,prob)] -> [(a,prob)] norm'' = - List.map (\ xs@((x,_):_) -> (x, sum (List.map snd xs))) . - List.groupBy (\x y -> fst x == fst y) . sortElem + List.map (\xs -> + case xs of + ((x,_):_) -> (x, sum (List.map snd xs)) + _ -> error "Probability.Distribution.norm'': every sub-list in groupBy must be non-empty") . + ListHT.groupBy (equating fst) . sortElem -- | pretty printing @@ -164,16 +180,48 @@ (//%) :: (Ord a, Show a) => T Rational a -> () -> IO () (//%) x () = putStr (pretty show x) -instance (Num prob, Ord prob, Ord a, Show a) => +instance (Num prob, Ord prob, Show prob, Ord a, Show a) => Show (T prob a) where showsPrec p (Cons xs) = showParen (p>10) (showString "fromFreqs " . showsPrec 11 (sortP (norm' xs))) + +{- | +We would like to have an equality test of type + +> (==) :: T prob a -> T prob a -> T prob Bool + +that is consistent with the 'Num' instance. +We would certainly define + +> x==y = norm (liftM2 (==) x y) . + +However the 'Eq' class enforces the type + +> T prob a -> T prob a -> Bool . + +We could implement this as check for equal distributions. +This would be inconsistent with the 'Num' instance +because it compares entire distributions, +not only individual outcomes. +Thus we provide this function as 'equal'. + +I would prefer to omit the 'Eq' instance completely, +but unfortunately the 'Num' instance requires 'Eq' as superclass. +-} instance Eq (T prob a) where - (==) = error "Probability.Dist.== cannot be implemented sensibly. It only exists for Num instance. Haskell98's numeric type class hierarchy sucks." + (==) = error "Probability.Distribution.== cannot be implemented sensibly." {- +instance (Num prob, Ord a) => Eq (T prob a) where + (==) = equal +-} + +equal :: (Num prob, Eq prob, Ord a) => T prob a -> T prob a -> Bool +equal = equating (decons . norm) + +{- The Num operations consider their operands as independent distributions (like all operations on distributions do). All functions normalize their results if normalization is lost by the plain operation. @@ -271,8 +319,10 @@ cond b d d' = b >>= \c -> if c then d else d' truth :: (Num prob) => T prob Bool -> prob -truth (Cons ((b,p):_:[])) = if b then p else 1-p -truth (Cons _) = error "Probability.truth: corrupt boolean random variable" +truth dist = + case decons (norm dist) of + (b,p):_ -> if b then p else 1-p + [] -> error "Probability.truth: corrupt boolean random variable" infixl 1 >>=? @@ -310,6 +360,11 @@ filter :: (Fractional prob) => (a -> Bool) -> T prob a -> T prob a filter p = fromFreqs . List.filter (p . fst) . decons + +mapMaybe :: (Fractional prob) => + (a -> Maybe b) -> T prob a -> T prob b +mapMaybe f = + fromFreqs . Maybe.mapMaybe (\(x,p) -> fmap (flip (,) p) $ f x) . decons -- | selecting from distributions diff --git a/src/Numeric/Probability/Either.hs b/src/Numeric/Probability/Either.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Probability/Either.hs @@ -0,0 +1,37 @@ +module Numeric.Probability.Either where + +import Control.Monad.Instances () + +import Control.Applicative (Applicative, pure, (<*>), liftA2, ) + + +newtype EitherT a m b = EitherT (m (Either a b)) + +instance Functor m => Functor (EitherT a m) where + fmap f (EitherT m) = EitherT $ fmap (fmap f) m + +instance Applicative m => Applicative (EitherT a m) where + pure a = EitherT $ pure $ Right a + EitherT af <*> EitherT am = + EitherT $ + liftA2 (\ef em -> + case ef of + Left b -> Left b + Right f -> + case em of + Left b -> Left b + Right m -> Right $ f m) af am + +instance Monad m => Monad (EitherT a m) where + return a = EitherT $ return $ Right a + EitherT m >>= f = EitherT $ do + e <- m + case e of + Left b -> return $ Left b + Right a -> + case f a of + EitherT n -> n + fail s = EitherT $ fail s + +throw :: Applicative m => a -> EitherT a m b +throw a = EitherT $ pure $ Left a diff --git a/src/Numeric/Probability/Example/Collection.hs b/src/Numeric/Probability/Example/Collection.hs --- a/src/Numeric/Probability/Example/Collection.hs +++ b/src/Numeric/Probability/Example/Collection.hs @@ -11,7 +11,7 @@ import Control.Monad.Trans.State (StateT(StateT, runStateT), evalStateT, ) import Control.Monad (liftM2, replicateM, ) -import qualified Data.List as List +import qualified Data.List.HT as ListHT import System.Random (Random) @@ -21,16 +21,10 @@ type Probability = Rational -{- | -see also the proposal - --} selectOne :: (Fractional prob) => StateT (Collection a) (Dist.T prob) a selectOne = - StateT $ \c -> - Dist.uniform $ init $ - zipWith (\xs (y:ys) -> (y, xs++ys)) (List.inits c) (List.tails c) + StateT $ Dist.uniform . ListHT.removeEach select1 :: (Fractional prob) => Collection a -> Dist.T prob a select1 = evalStateT selectOne diff --git a/src/Numeric/Probability/Example/DiceAccum.hs b/src/Numeric/Probability/Example/DiceAccum.hs --- a/src/Numeric/Probability/Example/DiceAccum.hs +++ b/src/Numeric/Probability/Example/DiceAccum.hs @@ -8,7 +8,6 @@ -} module Numeric.Probability.Example.DiceAccum where -import qualified Numeric.Probability.Example.Dice as Dice import qualified Numeric.Probability.Random as Rnd import qualified Numeric.Probability.Distribution as Dist import qualified Numeric.Probability.Transition as Trans diff --git a/src/Numeric/Probability/Example/Kruskal.hs b/src/Numeric/Probability/Example/Kruskal.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Probability/Example/Kruskal.hs @@ -0,0 +1,302 @@ +{- | +Given a row of n (~50) dice and +two players starting with a random dice within the first m (~5) dice. +Every players moves along the row, according the pips on the dice. +They stop if a move would exceed the row. +What is the probability that they stop at the same die? +(It is close to one.) + +Wuerfelschlange (german) +http://faculty.uml.edu/rmontenegro/research/kruskal_count/kruskal.html + +Kruskal's trick +http://www.math.de/exponate/wuerfelschlange.html/ +-} +module Numeric.Probability.Example.Kruskal where + +import qualified Numeric.Probability.Distribution as Dist +import qualified Numeric.Probability.Transition as Trans +import qualified Numeric.Probability.Random as Random +import qualified Numeric.Probability.Object as Obj +import qualified System.Random as Rnd +import qualified Text.Printf as P +import qualified Data.List as List +import Control.Monad (replicateM, ) +import Data.Function.HT (nest, compose2, ) +import Data.Tuple.HT (mapSnd, ) +import Data.Bool.HT (if', ) + + +type Die = Int + +type Probability = Rational +type Dist = Dist.T Probability + +die :: + (Obj.C prob experiment, Fractional prob) => + Score -> experiment Die +die maxPips = Obj.uniform [1..maxPips] + + + +type Score = Int + +{- | +We reformulate the problem to the following game: +There are two players, both of them collect a number of points. +In every round the player with the smaller score throws a die +and adds the pips to his score. +If the two players somewhen get the same score, then the game ends +and the score is the result of the game (@Just score@). +If one of the players exceeds the maximum score n, +then the game stops and players lose (@Nothing@). +-} +game :: + (Obj.C prob experiment, Fractional prob) => + Score -> Score -> (Score,Score) -> experiment (Maybe Score) +game maxPips maxScore = + let go (x,y) = + if maxScore < max x y + then return Nothing + else case compare x y of + EQ -> return (Just x) + LT -> do + d <- die maxPips + go (x+d, y) + GT -> do + d <- die maxPips + go (x, y+d) + in go + +gameRound :: + Score -> Score -> + Dist (Either (Maybe Score) (Score,Score)) -> + Dist (Either (Maybe Score) (Score,Score)) +gameRound maxPips maxScore current = Dist.norm $ do + e <- current + case e of + Left end -> return $ Left end + Right (x,y) -> + if maxScore < max x y + then return $ Left Nothing + else case compare x y of + EQ -> return $ Left (Just x) + LT -> do + d <- die maxPips + return $ Right (x+d, y) + GT -> do + d <- die maxPips + return $ Right (x, y+d) + +gameFast :: Score -> Score -> Dist (Score,Score) -> Dist (Maybe Score) +gameFast maxPips maxScore start = + Dist.mapMaybe (either Just (error "the game must be finished after maxScore moves")) $ + nest (maxScore+1) (gameRound maxPips maxScore) (fmap Right start) + +gameFastEither :: Score -> Score -> Dist (Score,Score) -> Dist (Maybe Score) +gameFastEither maxPips maxScore = Trans.untilLeft $ \(x,y) -> + if maxScore < max x y + then return $ Left Nothing + else case compare x y of + EQ -> return $ Left (Just x) + LT -> do + d <- die maxPips + return $ Right (x+d, y) + GT -> do + d <- die maxPips + return $ Right (x, y+d) + +{- | +This version could be generalized +to both Random and Distribution monad +while remaining efficient. +-} +gameFastFix :: Score -> Score -> Dist (Score,Score) -> Dist (Maybe Score) +gameFastFix maxPips maxScore = + Trans.fix $ \go (x,y) -> + if maxScore < max x y + then return Nothing + else case compare x y of + EQ -> return (Just x) + LT -> do + d <- die maxPips + go (x+d, y) + GT -> do + d <- die maxPips + go (x, y+d) + +{- | +In 'gameFastFix' we group the scores by rounds. +This leads to a growing probability distribution, +but we do not need the round number. +We could process the game in a different way: +We only consider the game states +where the lower score matches the round number. +-} +gameLeastScore :: Score -> Score -> Dist (Score,Score) -> Dist (Maybe Score) +gameLeastScore maxPips maxScore = + (Trans.fix $ \go (n,(x,y)) -> + if n > maxScore + then return Nothing + else + let next = go . (,) (succ n) + in case (x==n, y==n) of + (False, False) -> next (x,y) + (True, False) -> do + d <- die maxPips + next (x+d, y) + (False, True) -> do + d <- die maxPips + next (x, y+d) + (True, True) -> return (Just x)) + . + fmap ((,) 0) + +{- | +'gameLeastScore' can be written in terms of a matrix power. +For n pips we need a n² × n² matrix. +Using symmetries, we reduce it to a square matrix with size n·(n+1)/2. + +/ p[n+1,(n+1,n+1)] \ / p[n,(n+0,n+0)] \ +| p[n+1,(n+1,n+2)] | | p[n,(n+0,n+1)] | +| p[n+1,(n+1,n+3)] | | p[n,(n+0,n+2)] | +| ... | | ... | +| p[n+1,(n+1,n+6)] | = M/6 · | p[n,(n+0,n+5)] | +| p[n+1,(n+2,n+2)] | | p[n,(n+1,n+1)] | +| ... | | ... | +| p[n+1,(n+2,n+6)] | | p[n,(n+1,n+5)] | +| ... | | ... | +\ p[n+1,(n+6,n+6)] / \ p[n,(n+5,n+5)] / + +M[(n+1,(x,y)),(n,(x,y))] = 6 + +M[(n+1,(min y (n+d), max y (n+d))), (n,(n,y))] = 1 + +M[(n+1,(x1,y1)),(n,(x0,y0))] = 0 +-} +flattenedMatrix :: Score -> [Int] +flattenedMatrix maxPips = do + x1 <- [1..maxPips] + y1 <- [x1..maxPips] + x0 <- [0..maxPips-1] + y0 <- [x0..maxPips-1] + return $ + if' ((x0,y0) == (x1,y1)) maxPips $ + if' (x0==0 && (y0==x1 || y0==y1)) 1 0 + +{- +let e0 = [1,0,0,...,0] + +The cumulated probability is +e0 * (I + M + M^2 + ... + M^(n-1)) * startVector +and with M = V*D*V^-1 we get +e0 * V * (I + D + D^2 + ... + D^(n-1)) * V^-1 * startVector + +e0 * (I - M^n) * (I-M)^(-1) * startVector +-} +startVector :: Score -> [Int] +startVector maxPips = do + x <- [0..maxPips-1] + y <- [x..maxPips-1] + return $ if' (x==y) 1 2 + + +compareMaybe :: (Ord a) => Maybe a -> Maybe a -> Ordering +compareMaybe Nothing _ = GT +compareMaybe _ Nothing = LT +compareMaybe (Just a) (Just b) = compare a b + +cumulate :: (Ord a) => Dist (Maybe a) -> [(Maybe a, Probability)] +cumulate = + uncurry zip . mapSnd (scanl1 (+)) . unzip . + List.sortBy (compose2 compareMaybe fst) . Dist.decons + +runExact :: Score -> IO () +runExact maxPips = + mapM_ (\(m,p) -> + case m of + Just n -> + P.printf "%4d %7.2f %s\n" + n (fromRational (100*p) :: Double) (show p) + Nothing -> putStrLn $ "total: " ++ show p) $ + cumulate $ + gameFastFix maxPips 120 $ fmap ((,) 0) $ die maxPips + + +trace :: Score -> [Score] -> [Score] +trace s xs@(x:_) = s : trace (s+x) (drop x xs) +trace s [] = [s] + +chop :: [Score] -> [[Score]] +chop [] = [] +chop xs@(x:_) = uncurry (:) $ mapSnd chop $ splitAt x xs + +meeting :: [Score] -> [Score] -> Maybe Score +meeting xt@(x:xs) yt@(y:ys) = + case compare x y of + LT -> meeting xs yt + GT -> meeting xt ys + EQ -> Just x +meeting _ _ = Nothing + +{- | +This is a bruteforce implementation of the original game: +We just roll the die @maxScore@ times +and then jump from die to die according to the number of pips. +-} +bruteforce :: Score -> Score -> (Score,Score) -> Random.T (Maybe Score) +bruteforce maxPips maxScore (x,y) = do + points <- replicateM maxScore $ die maxPips + let run s = trace s (drop s points) + return (meeting (run x) (run y)) + +runSimulation :: Score -> IO () +runSimulation maxPips = + mapM_ (\(m,p) -> + case m of + Just n -> + P.printf "%4d %7.2f\n" + n (fromRational (100*p) :: Double) + Nothing -> putStrLn $ "total: " ++ show p) $ + cumulate $ + Random.runSeed (Rnd.mkStdGen 42) $ Random.dist $ + replicate 100000 $ bruteforce maxPips 120 . (,) 0 =<< die maxPips + + +latexDie :: Score -> String +latexDie pips = + "\\epsdice{" ++ show pips ++ "}" + +latexMarkedDie :: Score -> String +latexMarkedDie pips = + "\\epsdice[black]{" ++ show pips ++ "}" + +latexFromChain :: [Score] -> String +latexFromChain = + unlines . map latexDie + +latexChoppedFromChain :: [Score] -> String +latexChoppedFromChain = + unlines . + concatMap (\(p:ps) -> latexMarkedDie p : map latexDie ps) . + chop + +makeChains :: IO () +makeChains = do + let chains = + map + (\seed -> + Random.runSeed (Rnd.mkStdGen seed) $ + replicateM 42 $ die 6) + [30..42] + writeFile "KruskalDice.tex" $ + "\\noindent\n" + ++ + (List.intercalate "\\\\[4ex]\n" $ + map (("$\\rightarrow$" ++) . latexFromChain) chains) + ++ + "\\newpage\n" ++ + "\\noindent\n" + ++ + (List.intercalate "\\\\[4ex]\n" $ + map (("$\\rightarrow$" ++) . latexChoppedFromChain) chains) diff --git a/src/Numeric/Probability/Monad.hs b/src/Numeric/Probability/Monad.hs --- a/src/Numeric/Probability/Monad.hs +++ b/src/Numeric/Probability/Monad.hs @@ -1,16 +1,13 @@ -- | Monad helper functions module Numeric.Probability.Monad where +import Control.Monad.HT ((<=<), ) import Control.Monad (liftM, ) import Prelude hiding (iterate, ) --- | binary composition, available in GHC-6.8 Control.Monad -(>=>) :: Monad m => (a -> m b) -> (b -> m c) -> a -> m c -f >=> g = (>>= g) . f - -- | composition of a list of monadic functions compose :: Monad m => [a -> m a] -> a -> m a -compose = foldl (>=>) return +compose = foldl (flip (<=<)) return iterate :: Monad m => Int -> (a -> m a) -> (a -> m a) @@ -18,6 +15,8 @@ -- | like 'iterate' but returns all intermediate values +-- like Control.Monad.HT.iterateLimit, but counts differently. +-- Is it actually useful this way? walk :: (Monad m) => Int -> (a -> m a) -> (a -> m [a]) walk n f = let recourse 0 _ = return [] diff --git a/src/Numeric/Probability/Object.hs b/src/Numeric/Probability/Object.hs --- a/src/Numeric/Probability/Object.hs +++ b/src/Numeric/Probability/Object.hs @@ -1,4 +1,7 @@ -{-# OPTIONS_GHC -fglasgow-exts #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE UndecidableInstances #-} {- | Portability: Multi-parameter type class with functional dependency @@ -15,10 +18,13 @@ import qualified Numeric.Probability.Distribution as Dist import qualified Numeric.Probability.Random as Rnd import qualified Numeric.Probability.Shape as Shape +import qualified Numeric.Probability.Either as PE import qualified Data.List as List +import Control.Monad (liftM, ) + class Monad obj => C prob obj | obj -> prob where fromFrequencies :: [(a,prob)] -> obj a @@ -28,6 +34,10 @@ instance Fractional prob => C prob (Dist.T prob) where fromFrequencies = Dist.fromFreqs + +instance C prob obj => C prob (PE.EitherT b obj) where + fromFrequencies = + PE.EitherT . liftM Right . fromFrequencies diff --git a/src/Numeric/Probability/Percentage.hs b/src/Numeric/Probability/Percentage.hs --- a/src/Numeric/Probability/Percentage.hs +++ b/src/Numeric/Probability/Percentage.hs @@ -21,7 +21,7 @@ percent x = Cons (x/100) -showPfix :: (RealFrac prob) => Int -> prob -> String +showPfix :: (RealFrac prob, Show prob) => Int -> prob -> String showPfix precision f = if precision==0 then showR 3 (round (f*100) :: Integer)++"%" diff --git a/src/Numeric/Probability/Random.hs b/src/Numeric/Probability/Random.hs --- a/src/Numeric/Probability/Random.hs +++ b/src/Numeric/Probability/Random.hs @@ -7,7 +7,8 @@ import qualified System.Random as Random import System.Random (Random, ) -import Control.Monad.Trans.State (State, state, evalState, ) +import Control.Applicative (Applicative(..)) +import Control.Monad.Trans.State (State, state, evalState, runState, ) import qualified System.IO as IO import Prelude hiding (print) @@ -26,6 +27,9 @@ instance Functor T where fmap f = Cons . fmap f . decons +instance Applicative T where + pure x = Cons (pure x) + fm <*> m = Cons (decons fm <*> decons m) randomR :: Random.Random a => (a, a) -> T a randomR rng = @@ -35,7 +39,7 @@ Run random action in 'IO' monad. -} run :: T a -> IO a -run r = fmap (evalState (decons r)) Random.getStdGen +run = Random.getStdRandom . runState . decons {- | Run random action without 'IO' using a seed. @@ -66,7 +70,9 @@ {- | 'dist' converts a list of randomly generated values into -a distribution by taking equal weights for all values +a distribution by taking equal weights for all values. +Thus @dist (replicate n rnd)@ simulates @rnd@ @n@ times +and returns an estimation of the distribution represented by @rnd@. -} dist :: (Fractional prob, Ord a) => [T a] -> Distribution prob a dist = fmap (Dist.norm . Dist.uniform) . sequence diff --git a/src/Numeric/Probability/Show.hs b/src/Numeric/Probability/Show.hs --- a/src/Numeric/Probability/Show.hs +++ b/src/Numeric/Probability/Show.hs @@ -1,12 +1,19 @@ module Numeric.Probability.Show where +import qualified Data.List.HT as ListHT + +{- showL :: Show a => Int -> a -> String showL n x = s ++ replicate (n-length s) ' ' where s=show x +showL n x = + ListHT.padRight1 ' ' n (show x) +-} + showR :: Show a => Int -> a -> String -showR n x = replicate (n-length s) ' ' ++ s - where s=show x +showR n x = + ListHT.padLeft ' ' n (show x) --showP :: Float -> String --showP f = showR 3 (round (f*100))++"%" diff --git a/src/Numeric/Probability/Transition.hs b/src/Numeric/Probability/Transition.hs --- a/src/Numeric/Probability/Transition.hs +++ b/src/Numeric/Probability/Transition.hs @@ -3,6 +3,11 @@ import qualified Numeric.Probability.Distribution as Dist +import qualified Numeric.Probability.Either as PE + +import qualified Data.Map as Map + +import qualified Data.List.HT as ListHT import qualified Data.List as List import Prelude hiding (map, maybe, id, ) @@ -49,6 +54,34 @@ [T prob a] -> T prob a compose = foldl (\acc x v -> Dist.norm (acc v >>= x)) return + +untilLeft :: (Num prob, Ord a, Ord b) => + (a -> Dist.T prob (Either b a)) -> Dist.T prob a -> Dist.T prob b +untilLeft f = + let go final dist = + if null (Dist.decons dist) + then Dist.Cons $ Map.toList final + else + case ListHT.unzipEithers $ + List.map (\(e,p) -> either (\l -> Left (l,p)) (\r -> Right (r,p)) e) $ + Dist.decons $ Dist.norm $ dist >>= f of + (newFinal, stillActive) -> + go (Map.unionWith (+) (Map.fromListWith (+) newFinal) final) $ + Dist.Cons stillActive + in go Map.empty + +{- | +In @fix $ \go a -> do ...; go xy@ +any action after a 'go' is ignored. +-} +fix :: (Num prob, Ord a, Ord b) => + ((a -> PE.EitherT a (Dist.T prob) b) -> + (a -> PE.EitherT a (Dist.T prob) b)) -> + Dist.T prob a -> Dist.T prob b +fix f = + untilLeft $ \a -> + case f PE.throw a of + PE.EitherT m -> fmap (either Right Left) m -- * Spreading changes into transitions diff --git a/src/Numeric/Probability/Visualize.hs b/src/Numeric/Probability/Visualize.hs --- a/src/Numeric/Probability/Visualize.hs +++ b/src/Numeric/Probability/Visualize.hs @@ -9,6 +9,7 @@ import qualified Numeric.Probability.Distribution as Dist +import Control.Monad (when, ) import Data.List (nub, sort, ) @@ -204,11 +205,11 @@ title fe++"\",xlab=\""++ xLabel fe++"\",ylab=\""++ yLabel fe++"\")") - mapM out1' (zipWith3 drawy [1..length ys'] ps ys') - if null (concatMap label ps) - then return () - else out1' $ legend (incr minx) maxy ps - out1' ("dev2bitmap(\""++(fileName fe)++".pdf\", type=\"pdfwrite\")") + mapM_ out1' (zipWith3 drawy [1 ..] ps ys') + when (not $ null $ concatMap label ps) $ + out1' $ legend (incr minx) maxy ps + out1' ("dev2bitmap(" ++ show (fileName fe ++ ".pdf") ++ + ", type=\"pdfwrite\")") {-