probability 0.1 → 0.2
raw patch · 44 files changed
+2552/−1922 lines, 44 filesdep +mtldep −haskell98
Dependencies added: mtl
Dependencies removed: haskell98
Files
- Alarm.hs +0/−58
- Barber.hs +0/−56
- Bayesian.hs +0/−95
- Boys.hs +0/−47
- Collection.hs +0/−112
- Dice.hs +0/−39
- ListUtils.hs +0/−16
- MontyHall.hs +0/−79
- NBoys.hs +0/−46
- Predator.hs +0/−87
- PrintList.hs +0/−54
- Probability.hs +0/−655
- Queuing.hs +0/−144
- Show.hs +0/−16
- ToDo +26/−12
- TreeGrowth.hs +0/−143
- Visualize.hs +0/−241
- probability.cabal +32/−22
- src/Numeric/Probability/Distribution.hs +342/−0
- src/Numeric/Probability/Example/Alarm.hs +62/−0
- src/Numeric/Probability/Example/Barber.hs +69/−0
- src/Numeric/Probability/Example/Bayesian.hs +101/−0
- src/Numeric/Probability/Example/Boys.hs +55/−0
- src/Numeric/Probability/Example/Collection.hs +133/−0
- src/Numeric/Probability/Example/Diagnosis.hs +60/−0
- src/Numeric/Probability/Example/Dice.hs +43/−0
- src/Numeric/Probability/Example/DiceAccum.hs +77/−0
- src/Numeric/Probability/Example/MontyHall.hs +112/−0
- src/Numeric/Probability/Example/NBoys.hs +41/−0
- src/Numeric/Probability/Example/Predator.hs +87/−0
- src/Numeric/Probability/Example/Queuing.hs +159/−0
- src/Numeric/Probability/Example/TreeGrowth.hs +151/−0
- src/Numeric/Probability/Expectation.hs +57/−0
- src/Numeric/Probability/Monad.hs +53/−0
- src/Numeric/Probability/Object.hs +63/−0
- src/Numeric/Probability/Percentage.hs +127/−0
- src/Numeric/Probability/PrintList.hs +54/−0
- src/Numeric/Probability/Random.hs +102/−0
- src/Numeric/Probability/Shape.hs +32/−0
- src/Numeric/Probability/Show.hs +16/−0
- src/Numeric/Probability/Simulation.hs +96/−0
- src/Numeric/Probability/Trace.hs +50/−0
- src/Numeric/Probability/Transition.hs +108/−0
- src/Numeric/Probability/Visualize.hs +244/−0
− Alarm.hs
@@ -1,58 +0,0 @@-module Alarm where--import Probability (Dist, Probability, choose, (??), (|||))--type PBool = Dist Bool---flp :: Float -> PBool-flp p = choose p True False----- * Alarm network---- | prior burglary 1%-b :: PBool-b = flp 0.01---- | prior earthquake 0.1%-e :: PBool-e = flp 0.001---- | conditional probability of alarm given burglary and earthquake-a :: Bool -> Bool -> PBool-a b0 e0 =- case (b0,e0) of- (False,False) -> flp 0.01- (False,True) -> flp 0.1- (True,False) -> flp 0.7- (True,True) -> flp 0.8----- | conditional probability of john calling given alarm-j :: Bool -> PBool-j a0 = if a0 then flp 0.8 else flp 0.05---- | conditional probability of mary calling given alarm-m :: Bool -> PBool-m a0 = if a0 then flp 0.9 else flp 0.1---- | calculate the full joint distribution-data Burglary = B { burglary :: Bool,- earthquake :: Bool,- alarm :: Bool,- john :: Bool,- mary :: Bool }- deriving (Eq, Ord, Show)--bJoint :: Dist Burglary-bJoint = do b' <- b -- burglary- e' <- e -- earthquake- a' <- a b' e' -- alarm- j' <- j a' -- john- m' <- m a' -- mary- return (B b' e' a' j' m')---- | what is the probability that mary calls given that john calls?-pmj :: Probability-pmj = mary ?? bJoint ||| john
− Barber.hs
@@ -1,56 +0,0 @@-module Barber where--import Probability (Dist, RDist, Trans, normal)-import Queuing (Time, System, unit, evalSystem, idleAvgP, waiting)---- * barber shop--custServ :: Dist Time-custServ = normal [5..10]--nextCust :: Trans Time -- not dependant on serving time-nextCust _ = normal [3..6]--barbers :: Int-barbers = 1--customers :: Int-customers = 20--runs :: Int-runs = 50--barberEvent :: ((), (Dist Time, Time -> Dist Time))-barberEvent = unit (custServ, nextCust)--barberEvents :: [((), (Dist Time, Time -> Dist Time))]-barberEvents = replicate customers barberEvent--barberSystem :: (System () -> b) -> RDist b-barberSystem eval = evalSystem runs barbers barberEvents eval----- * category--data Category = ThreeOrLess | FourToTen | MoreThanTen- deriving (Eq,Ord,Show)--cat :: Time -> Category-cat n | n <= 3 = ThreeOrLess-cat n | n <= 10 = FourToTen-cat _ = MoreThanTen--perc :: Float -> String-perc n | n <= 0.25 = "0% to 25%"-perc n | n <= 0.5 = "25% to 50%"-perc n | n <= 0.75 = "50% to 75%"-perc _ = "75% to 100%"---- * evaluation---- | avg barber idle time-barberIdle :: RDist String-barberIdle = barberSystem (perc.(idleAvgP barbers))--- | avg customer waiting time (unserved customers)-customerWait :: RDist Category-customerWait = barberSystem ( cat.(`div` customers).(waiting barbers) )
− Bayesian.hs
@@ -1,95 +0,0 @@-module Bayesian where--import Probability (Dist, Probability, ProbRep, maybeT, sequ, (??), (|||))---{---Approach: model a node with k predecessors as a function with k- parameters---}------ * Abbreviations, smart constructors--type State a = [a]-type PState a = Dist (State a)-type STrans a = State a -> PState a-type SPred a = a -> State a -> Bool--event :: ProbRep -> a -> STrans a-event p e0 = maybeT p (e0:)--happens :: Eq a => SPred a-happens = elem--network :: [STrans a] -> PState a-network = flip sequ []---source :: ProbRep -> a -> STrans a-source = event--bin :: Eq a => a -> a -> ProbRep -> ProbRep -> ProbRep -> ProbRep -> a -> STrans a-bin x y a b c d z s | elem x s && elem y s = event a z s- | elem x s = event b z s- | elem y s = event c z s- | otherwise = event d z s----- | Two possible causes for one effect--data Nodes = A | B | E deriving (Eq,Ord,Show)--g :: PState Nodes-g = network [source 0.1 A,- source 0.2 B,- bin A B 1 1 0.5 0 E]---- * queries--e, aE, bE :: Probability-e = happens E ?? g-aE = happens A ?? g ||| happens E-bE = happens B ?? g ||| happens E---{--data State = State {causeA :: Bool, causeB :: Bool, effect :: Bool}- deriving (Eq,Ord,Show)--nCauseA s = s{causeA=True}--}------- Wet grass example------ cloudy = true 0.5------ sprinkler c = dep c 0.1 0.5------ rain c = dep c 0.8 0.2------ wetGrass s r = bin s r 0.99 0.9 0.9 0------ c = cloudy--- s = sprinkler cloudy--- r = rain cloudy--- w = wetGrass s r----- alarm :: Prob -> Prob -> Prob--- alarm b e = cond b (pTrue 0.8)--- (cond e (pTrue 0.1) (pTrue 0.01))------ john :: Prob -> Prob--- john a = cond a (pTrue 0.7) (pTrue 0.1)------ mary :: Prob -> Prob--- mary a = cond a (pTrue 0.6) (pTrue 0.2)--------- maryWhenJohn = mary a ?? john a--- where a = alarm (pTrue 0.5) (pTrue 0.1)
− Boys.hs
@@ -1,47 +0,0 @@-{- |-Consider a family of two children. Given that there is a boy in the family,-what is the probability that there are two boys in the family?--}--module Boys where--import Probability- (Dist, Probability, Trans, Event,- uniform, just, mapD, sequ, (??), (|||))---data Child = Boy | Girl- deriving (Eq,Ord,Show)--type Family = [Child]--birth :: Trans Family-birth f = uniform [Boy:f,Girl:f]--family :: Dist Family-family = sequ [birth,birth] []---- NOTE: could be fixed to 2--- could be renamed to allBoys----boys :: Int -> Event Family-boys n = just (replicate n Boy)--existsBoy :: Event Family-existsBoy = elem Boy---- NOTE: might not be needed, i.e., definition can be inlined instead----familyWithBoy :: Dist Family-familyWithBoy = family ||| existsBoy--twoBoys :: Probability-twoBoys = (boys 2) ?? familyWithBoy---countBoys :: Family -> Int-countBoys = length . filter (==Boy)--numBoys :: Dist Int-numBoys = mapD countBoys familyWithBoy-
− Collection.hs
@@ -1,112 +0,0 @@-module Collection where--import Probability- (Dist, Probability, uniform, mapD, (??), oneOf, while, random, (~.))--import qualified List (delete)----type Collection a = [a]---- this is a StateT-selectOne :: Eq a => Collection a -> Dist (a,Collection a)-selectOne c = uniform [(v,List.delete v c) | v <- c]--select1 :: Eq a => Collection a -> Dist a-select1 = mapD fst . selectOne--select2 :: Eq a => Collection a -> Dist (a,a)-select2 c = do (x,c') <- selectOne c- y <- select1 c'- return (x,y)---- this is a replicateM with respect to StateT-selectMany :: Eq a => Int -> Collection a -> Dist ([a],Collection a)-selectMany 0 c = return ([],c)-selectMany n c = do (x,c1) <- selectOne c- (xs,c2) <- selectMany (n-1) c1- return (x:xs,c2)--select :: Eq a => Int -> Collection a -> Dist [a]-select n = mapD (reverse . fst) . selectMany n----- * Example collections---- ** marbles--data Marble = R | G | B deriving (Eq,Ord,Show)--bucket :: Collection Marble-bucket = [R,R,R,R,R, G,G,G, B,B]--jar :: Collection Marble-jar = [R,R,G,G,B]---- pRGB = prob (just [R,G,B]) (select 3 bucket)-pRGB :: Probability-pRGB = (==[R,G,B]) ?? select 3 jar-pRG :: Probability-pRG = (oneOf [[R,G],[G,R]]) ?? select 2 jar---- ** cards--data Suit = Club | Spade | Heart | Diamond- deriving (Eq,Ord,Show,Enum)--data Rank = Plain Int | Jack | Queen | King | Ace- deriving (Eq,Ord,Show)--type Card = (Rank,Suit)--plains :: [Rank]-plains = map Plain [2..10]--faces :: [Rank]-faces = [Jack,Queen,King,Ace]--isFace :: Card -> Bool-isFace (r,_) = r `elem` faces--- isFace = (`elem` faces) . fst--isPlain :: Card -> Bool-isPlain (r,_) = r `elem` plains--ranks :: [Rank]-ranks = plains ++ faces--suits :: [Suit]-suits = [Club,Spade,Heart,Diamond]--deck :: Collection Card-deck = [(r,s) | r <- ranks, s <- suits]----- * Example--{- | mini-blackjack:-draw 2 cards, and if value is less than 14, continue drawing-until value equals or exceeds 14. if values exceeds 21,-you lose, otherwise you win.--}--value :: Card -> Int-value ((Plain n),_) = n-value (Ace,_) = 11-value _ = 10--draw :: ([Card], Collection Card) -> Dist ([Card], Collection Card)-draw (cards,cl) = fmap f (selectOne cl)- where- f (c,cl') = ((c:cards),cl')--drawTo16 :: t -> IO ([Card], Collection Card)-drawTo16 _ = while (\(cards,_)->(sum (map value cards) < 16))- (random draw) ([], deck)--win :: ([Card], b) -> Bool-win (cards,_) = sum (map value cards) <= 21--chanceWin :: IO (Dist Bool)-chanceWin = fmap (mapD win) ((100 ~. drawTo16) undefined)
− Dice.hs
@@ -1,39 +0,0 @@-module Dice where--import Probability (Dist, Probability, prod, uniform, (??))-import Monad (liftM2)---type Die = Int--die :: Dist Die-die = uniform [1..6]--twoDice :: Dist (Die,Die)-twoDice = prod die die--dice :: Int -> Dist [Die]-dice n = sequence $ replicate n die--- dice = replicateM---twoSixes :: Probability-twoSixes = (==(6,6)) ?? liftM2 (,) die die--{- |-@sixes p n@ computes the probability of getting-p sixes (@>1@, @==2@, ...) when rolling n dice--}-sixes :: (Int -> Bool) -> Int -> Probability-sixes p n = (p . length . filter (==6)) ?? dice n--droll :: Dist Die-droll =- liftM2 (+) (uniform [0,1]) die--g3 :: Probability-g3 = (>3) ?? die--addTwo :: Dist Die-addTwo =- liftM2 (+) die die
− ListUtils.hs
@@ -1,16 +0,0 @@-module ListUtils where----- | create a singleton list, you can also use 'return' for the list 'Monad'-singleton :: a -> [a]-singleton x = [x]----- | apply a function to the @n@th element of a list-onNth :: Int -> (a -> a) -> [a] -> [a]-onNth n f xs =- let (ys,zs) = splitAt n xs- in ys ++- case zs of- [] -> []- z:zs' -> f z : zs'
− MontyHall.hs
@@ -1,79 +0,0 @@-module MontyHall where---import Probability- (Dist, Trans, RDist, uniform, mapD, (~.), idT, sequ, certainly, )---import ListUtils (replicate)-import List ( (\\) )--- import Monad (liftM)--data Door = A | B | C- deriving (Eq,Ord,Show)--doors :: [Door]-doors = [A,B,C]--data State = Doors {prize :: Door, chosen :: Door, opened :: Door}- deriving (Eq,Ord,Show)----- | initial configuration of the game status-start :: State-start = Doors {prize=u,chosen=u,opened=u} where u=undefined---{- |-Steps of the game:-- (1) hide the prize-- (2) choose a door-- (3) open a non-open door, not revealing the prize-- (4) apply strategy: switch or stay--}-hide :: Trans State-hide s = uniform [s {prize = d} | d <- doors]--choose :: Trans State-choose s = uniform [s {chosen = d} | d <- doors]--open :: Trans State-open s = uniform [s {opened = d} | d <- doors \\ [prize s,chosen s]]--type Strategy = Trans State--switch :: Strategy-switch s = uniform [s {chosen = d} | d <- doors \\ [chosen s,opened s]]--stay :: Strategy-stay = idT--game :: Strategy -> Trans State-game s = sequ [hide,choose,open,s]----- * Playing the game--data Outcome = Win | Lose- deriving (Eq,Ord,Show)--result :: State -> Outcome-result s = if chosen s==prize s then Win else Lose--eval :: Strategy -> Dist Outcome-eval s = mapD result (game s start)--simEval :: Int -> Strategy -> RDist Outcome-simEval k s = mapD result `fmap` (k ~. game s) start----- * Alternative modeling--firstChoice :: Dist Outcome-firstChoice = uniform [Win,Lose,Lose]--switch' :: Trans Outcome-switch' Win = certainly Lose-switch' Lose = certainly Win
− NBoys.hs
@@ -1,46 +0,0 @@-{- |-Ceneralization of "Boys"--Consider a family of n children. Given that there are k boys in the family,-what is the probability that there are m boys in the family?--}--module NBoys where--import Probability- (Dist, Probability, Trans, Event,- uniform, mapD, sequ, (??), (|||))--data Child = Boy | Girl- deriving (Eq,Ord,Show)--type Family = [Child]--birth :: Trans Family-birth f = uniform [Boy:f,Girl:f]--family :: Int -> Dist Family-family n = sequ (replicate n birth) []--countBoys :: Family -> Int-countBoys = length . filter (==Boy)--boys :: Int -> Event Family-boys k f = countBoys f >= k--nBoys :: Int -> Int -> Int -> Probability-nBoys n k m = (boys m) ?? (family n ||| boys k)--numBoys :: Int -> Int -> Dist Int-numBoys n k = mapD countBoys (family n ||| boys k)-------- Special cases------- only boys in a family that has one boy----onlyBoys1 :: Int -> Probability-onlyBoys1 n = nBoys n 1 n-
− Predator.hs
@@ -1,87 +0,0 @@-{- |-Lotka-Volterra predator-prey model--parameters-- * @g@ : victims' growth factor-- * @d@ : predators' death factor-- * @s@ : search rate-- * @e@ : energetic efficiency--}--module Predator where--import Visualize (- Vis, Color(Green, Red),- figP, figure, title,- showParams, xLabel, yLabel, plotL, color, label,- )----- try: n>=500--- g = 1.05--- d = 0.95--- s = 0.01--- e = 0.01---g, d, s, e :: Float-g = 1.02-d = 0.98-s = 0.01-e = 0.01----- 'direct' function-over-time approach -- very inefficient due to recursion------ v :: Int -> Float--- v 0 = 20--- v t = ((1 + r - a*p(t-1)) * v (t-1)) `max` 0------ p :: Int -> Float--- p 0 = 15--- p t = ((1 - d + a*b*v(t-1)) * p (t-1)) `max` 0--------- fig1 = figP figure{title="Predator/Prey Simulation "++--- showParams [r,d,a,b] ["r","d","a","b"],--- xLabel="Time (generation)",--- yLabel="Population"}--- [(plotF (0,15,1) v){color=Green,label="Victim"},--- (plotF (0,15,1) p){color=Red,label="Prey"}]--v0 :: Float-v0 = 1--p0 :: Float-p0 = 1--dv :: (Float,Float) -> Float-dv (v,p) = (g*v - s*v*p) `max` 0--dp :: (Float,Float) -> Float-dp (v,p) = (d*p + e*v*p) `max` 0--dvp :: (Float, Float) -> (Float, Float)-dvp vp' = (dv vp', dp vp')--vp :: [(Float, Float)]-vp = (v0,p0):map dvp vp--vs :: [Float]-vs = map fst vp--ps :: [Float]-ps = map snd vp---fig1 :: Int -> Vis-fig1 n = figP figure{title="Predator/Prey Simulation "++- showParams [g,d,s,e] ["g","d","s","e"],- xLabel="Time (generation)",- yLabel="Population"}- [(plotL (take n vs)){color=Green,label="Victim"},- (plotL (take n ps)){color=Red,label="Prey"}]
− PrintList.hs
@@ -1,54 +0,0 @@--- | Utilities for printing lists-module PrintList where---import List (intersperse)---------------------------------------------------------------------------- PRINT UTILITIES-------------------------------------------------------------------------newtype Lines a = Lines [a]--instance Show a => Show (Lines a) where- show (Lines xs) = printList ("","\n","") show xs--asLines :: [a] -> Lines a-asLines = Lines---showNQ :: Show a => a -> String-showNQ = filter ('"'/=) . show--indent :: Int -> Int -> [Char]-indent i l = take (i*l) (repeat ' ')--printList :: ([a],[a],[a]) -> (b -> [a]) -> [b] -> [a]-printList (sep0,sep1,sep2) f xs =- sep0++concat (intersperse sep1 (map f xs))++sep2---asTuple, asSeq, asList, asSet, asLisp,- asString, asPlain, asPlain' :: (a -> [Char]) -> [a] -> [Char]--asTuple = printList ("(",",",")")-asSeq = printList ("",",","")-asList = printList ("[",",","]")-asSet = printList ("{",",","}")-asLisp = printList ("("," ",")")-asPlain f xs = if null xs then "" else printList (" "," ","") f xs-asPlain' f xs = if null xs then "" else printList (""," ","") f xs-asString = printList ("","","")--- asLines = printList ["","\n",""]--asCases :: Int -> (a -> [Char]) -> [a] -> [Char]-asCases l =- let ind = indent 4 l- in printList ("\n"++ind++" ","\n"++ind++" | ","")--asDefs :: [Char] -> (a -> [Char]) -> [a] -> [Char]-asDefs n = printList ("\n"++n,"\n"++n,"\n")--asParagraphs :: (a -> [Char]) -> [a] -> [Char]-asParagraphs = printList ("\n","\n\n","\n")
− Probability.hs
@@ -1,655 +0,0 @@-module Probability where--import qualified Random-import List (sort, sortBy, transpose)-import Monad (MonadPlus, mplus, mzero)--import ListUtils (singleton)-import Show (showR)---{- TO DO:--* create export list--* extend Dist by a constructor for continuous distributions:-- C (Float -> Float)--* prove correctness of |||--* Monad helpers into separate module---}----- * Auxiliary definitions---- ** Events-type Event a = a -> Bool--oneOf :: Eq a => [a] -> Event a-oneOf = flip elem--just :: Eq a => a -> Event a-just = oneOf . singleton----- ** Probabilities-newtype Probability = P ProbRep--type ProbRep = Float--precision :: Int-precision = 1--showPfix :: ProbRep -> String-showPfix f =- if precision==0- then showR 3 (round (f*100) :: Integer)++"%"- else showR (4+precision) (roundRel precision (f*100))++"%"--roundRel :: (RealFrac a) => Int -> a -> a-roundRel p x =- let d = 10^p- in fromIntegral (round (x*d) :: Integer)/d---- -- mixed precision--- ----- showP :: ProbRep -> String--- showP f | f>=0.1 = showR 3 (round (f*100))++"%"--- | otherwise = show (f*100)++"%"---- fixed precision----showP :: ProbRep -> String-showP = showPfix---instance Show Probability where- show (P p) = showP p--errorMargin :: ProbRep-errorMargin = 0.00001----- ** Monad composition---- | binary composition-(>@>) :: Monad m => (a -> m b) -> (b -> m c) -> a -> m c-f >@> g = (>>= g) . f---- | composition of a list of monadic functions-sequ :: Monad m => [a -> m a] -> a -> m a-sequ = foldl (>@>) return------ * Deterministic and probabilistic values---- ** Distributions---- | probability disribution-newtype Dist a = D {unD :: [(a,ProbRep)]}--instance Monad Dist where- return x = D [(x,1)]- d >>= f = D [(y,q*p) | (x,p) <- unD d, (y,q) <- unD (f x)]- fail _ = D []---- note: mzero is a zero for >>= and a unit for mplus----instance MonadPlus Dist where- mzero = D []- mplus d d' | isZero d || isZero d' = mzero- | otherwise = unfoldD $ choose 0.5 d d'--isZero :: Dist a -> Bool-isZero (D d) = null d---instance Functor Dist where- fmap f (D d) = D [(f x,p) | (x,p) <- d]--instance (Ord a,Eq a) => Eq (Dist a) where- D xs == D ys = map fst (norm' xs)==map fst (norm' ys) &&- all (\((_,p),(_,q))->abs (p-q)<errorMargin) (zip xs ys)----- *** Auxiliary functions for constructing and working with distributions-onD :: ([(a,ProbRep)] -> [(a,ProbRep)]) -> Dist a -> Dist a-onD f = D . f . unD--sizeD :: Dist a -> Int-sizeD = length . unD--checkD :: Dist a -> Dist a-checkD (D d) | abs (1-sumP d) < errorMargin = D d- | otherwise = error ("Illegal distribution: total probability = "++show (sumP d))--mkD :: [(a,ProbRep)] -> Dist a-mkD = checkD . D--sumP :: [(a,ProbRep)] -> ProbRep-sumP = sum . map snd--sortP :: [(a,ProbRep)] -> [(a,ProbRep)]-sortP = sortBy (\x y->compare (snd y) (snd x))----- *** Normalization = grouping-normBy :: Ord a => (a -> a -> Bool) -> Dist a -> Dist a-normBy f = onD $ accumBy f . sort--accumBy :: Num b => (a -> a -> Bool) -> [(a,b)] -> [(a,b)]-accumBy f ((x,p):ys@((y,q):xs)) | f x y = accumBy f ((x,p+q):xs)- | otherwise = (x,p):accumBy f ys-accumBy _ xs = xs--norm :: Ord a => Dist a -> Dist a-norm = normBy (==)--norm' :: Ord a => [(a,ProbRep)] -> [(a,ProbRep)]-norm' = accumBy (==) . sort----- pretty printing-instance (Ord a,Show a) => Show (Dist a) where- show (D []) = "Impossible"- show (D xs) = concatMap (\(x,p)->showR w x++' ':showP p++"\n") (sortP (norm' xs))- where w = maximum (map (length.show.fst) xs)----- *** Operations on distributions---- | product of independent distributions, identical to 'Monad.liftM2'-joinWith :: (a -> b -> c) -> Dist a -> Dist b -> Dist c-joinWith f (D d) (D d') = D [ (f x y,p*q) | (x,p) <- d, (y,q) <- d']--prod :: Dist a -> Dist b -> Dist (a,b)-prod = joinWith (,)----- ** Spread: functions to convert a list of values into a distribution---- | distribution generators-type Spread a = [a] -> Dist a--certainly :: Trans a-certainly = return--impossible :: Dist a-impossible = mzero--choose :: ProbRep -> a -> a -> Dist a-choose p x y = enum [p,1-p] [x,y]--enum :: [ProbRep] -> Spread a-enum ps xs = mkD $ zip xs ps--enumPC :: [ProbRep] -> Spread a-enumPC ps = enum (map (/100) ps)--relative :: [Int] -> Spread a-relative ns = enum (map (\n->fromIntegral n/fromIntegral (sum ns)) ns)--shape :: (Float -> Float) -> Spread a-shape _ [] = impossible-shape f xs = scale (zip xs ps)- where incr = 1 / fromIntegral ((length xs) - 1)- ps = map f (iterate (+incr) 0)--linear :: Float -> Spread a-linear c = shape (c*)--uniform :: Spread a-uniform = shape (const 1)--negexp :: Spread a-negexp = shape (\x -> exp (-x))--normal :: Spread a-normal = shape (normalCurve 0.5 0.5)--normalCurve :: Float -> Float -> Float -> Float-normalCurve mean dev x = 1 / sqrt (2 * pi) * exp (-1/2 * u^(2::Int))- where u = (x - mean) / dev----- | extracting and mapping the domain of a distribution-extract :: Dist a -> [a]-extract = map fst . unD--mapD :: (a -> b) -> Dist a -> Dist b-mapD = fmap----- | unfold a distribution of distributions into one distribution-unfoldD :: Dist (Dist a) -> Dist a-unfoldD (D d) = D [ (x,p*q) | (d',q) <- d, (x,p) <- unD d' ]----- | conditional distribution-cond :: Dist Bool -> Dist a -> Dist a -> Dist a-cond b d d' = unfoldD $ choose p d d'- where P p = truth b--truth :: Dist Bool -> Probability-truth (D ((b,p):_:[])) = P (if b then p else 1-p)-truth (D _) = error "Probability.truth: corrupt boolean random variable"----- | conditional probability-(|||) :: Dist a -> Event a -> Dist a-(|||) = flip filterD----- | filtering distributions-data Select a = Case a | Other- deriving (Eq,Ord,Show)--above :: Ord a => ProbRep -> Dist a -> Dist (Select a)-above p (D d) = D (map (\(x,q)->(Case x,q)) d1++[(Other,sumP d2)])- where (d1,d2) = span (\(_,q)->q>=p) (sortP (norm' d))--scale :: [(a,ProbRep)] -> Dist a-scale xs = D (map (\(x,p)->(x,p/q)) xs)- where q = sumP xs--filterD :: (a -> Bool) -> Dist a -> Dist a-filterD p = scale . filter (p . fst) . unD----- | selecting from distributions-selectP :: Dist a -> ProbRep -> a-selectP (D d) p = scanP p d--scanP :: ProbRep -> [(a,ProbRep)] -> a-scanP p ((x,q):ps) =- if p<=q || null ps- then x- else scanP (p-q) ps-scanP _ [] = error "Probability.scanP: distribution must be non-empty"--infix 8 ??--(??) :: Event a -> Dist a -> Probability-(??) p = P . sumP . filter (p . fst) . unD----- TO DO: generalize Float to arbitrary Num type----class ToFloat a where- toFloat :: a -> Float--instance ToFloat Float where toFloat = id-instance ToFloat Int where toFloat = fromIntegral-instance ToFloat Integer where toFloat = fromIntegral--class FromFloat a where- fromFloat :: Float -> a--instance FromFloat Float where fromFloat = id-instance FromFloat Int where fromFloat = round-instance FromFloat Integer where fromFloat = round---- expected :: ToFloat a => Dist a -> Float--- expected = sum . map (\(x,p)->toFloat x*p) . unD--class Expected a where- expected :: a -> Float---- instance ToFloat a => Expected a where--- expected = toFloat-instance Expected Float where expected = id-instance Expected Int where expected = toFloat-instance Expected Integer where expected = toFloat--instance Expected a => Expected [a] where- expected xs = sum (map expected xs) / toFloat (length xs)--instance Expected a => Expected (Dist a) where- expected = sum . map (\(x,p)->expected x*p) . unD----- | statistical analyses-variance :: Expected a => Dist a -> Float-variance d@(D ps) = sum $ map (\(x,p)->p*sqr (expected x - ex)) ps- where sqr x = x * x- ex = expected d--stddev :: Expected a => Dist a -> Float-stddev = sqrt . variance------ * Randomized values----- ** R random value---- | Random values-type R a = IO a--printR :: Show a => R a -> R ()-printR = (>>= print)---- instance Show (IO a) where--- show _ = ""--pick :: Dist a -> R a--- pick d = do {p <- Random.randomRIO (0,1); return (selectP p d)}-pick d = Random.randomRIO (0,1) >>= return . selectP d----- ** RDist random distribution---- | Randomized distributions-type RDist a = R (Dist a)--rAbove :: Ord a => ProbRep -> RDist a -> RDist (Select a)-rAbove p rd = do D d <- rd- let (d1,d2) = span (\(_,q)->q>=p) (sortP (norm' d))- return (D (map (\(x,q)->(Case x,q)) d1++[(Other,sumP d2)]))------ * Deterministic and probabilistic generators---- ** Transitions----- | deterministic generator-type Change a = a -> a---- | probabilistic generator-type Trans a = a -> Dist a--idT :: Trans a-idT = certainlyT id----- mapT maps a change function to the result of a transformation--- (mapT is somehow a lifted form of mapD)--- The restricted type of f results from the fact that the--- argument to t cannot be changed to b in the result Trans type.----mapT :: Change a -> Trans a -> Trans a-mapT f t = mapD f . t----- unfold a distribution of transitions into one transition------ NOTE: The argument transitions must be independent----unfoldT :: Dist (Trans a) -> Trans a-unfoldT (D d) x = D [ (y,p*q) | (f,p) <- d, (y,q) <- unD (f x) ]----- ** Spreading changes into transitions---- | functions to convert a list of changes into a transition-type SpreadC a = [Change a] -> Trans a--certainlyT :: Change a -> Trans a-certainlyT f = certainly . f--- certainlyT = (certainly .)--- certainlyT = maybeC 1--maybeT :: ProbRep -> Change a -> Trans a-maybeT p f = enumT [p,1-p] [f,id]--liftC :: Spread a -> [Change a] -> Trans a-liftC s cs x = s [f x | f <- cs]--- liftC s cs x = s $ map ($ x) cs--uniformT :: [Change a] -> Trans a-uniformT = liftC uniform--normalT :: [Change a] -> Trans a-normalT = liftC normal--linearT :: Float -> [Change a] -> Trans a-linearT c = liftC (linear c)--enumT :: [ProbRep] -> [Change a] -> Trans a-enumT xs = liftC (enum xs)----- ** Spreading transitions into transitions---- | functions to convert a list of transitions into a transition-type SpreadT a = [Trans a] -> Trans a--liftT :: Spread (Trans a) -> [Trans a] -> Trans a-liftT s = unfoldT . s--uniformTT :: [Trans a] -> Trans a-uniformTT = liftT uniform--normalTT :: [Trans a] -> Trans a-normalTT = liftT normal--linearTT :: Float -> [Trans a] -> Trans a-linearTT c = liftT (linear c)--enumTT :: [ProbRep] -> [Trans a] -> Trans a-enumTT xs = liftT (enum xs)------ * Randomized generators---- ** Randomized changes---- | random change-type RChange a = a -> R a--random :: Trans a -> RChange a-random t = pick . t--- random = (pick .)----- ** Randomized transitions---- | random transition-type RTrans a = a -> RDist a-type ApproxDist a = R [a]---{- |-'rDist' converts a list of randomly generated values into-a distribution by taking equal weights for all values--}-rDist :: Ord a => [R a] -> RDist a-rDist = fmap (norm . uniform) . sequence------ * Iteration and simulation----- Iterate class defining *.--- Sim class defining ~.---{- |--Naming convention:-- * @*@ takes @n :: Int@ and a generator and iterates the generator n times-- * @.@ produces a single result-- * @..@ produces a trace-- * @~@ takes @k :: Int@ [and @n :: Int@] and a generator and simulates- the [n-fold repetition of the] generator k times---There are the following functions:-- * @n *. t@ iterates t and produces a distribution-- * @n *.. t@ iterates t and produces a trace-- * @k ~. t@ simulates t and produces a distribution-- * @(k,n) ~*. t@ simulates the n-fold repetition of t and produces a distribution-- * @(k,n) ~.. t@ simulates the n-fold repetition of t and produces a trace---Iteration captures three iteration strategies:-iter builds an n-fold composition of a (randomized) transition-while and until implement conditional repetitions--The class Iterate allows the overloading of iteration for different-kinds of generators, namely transitions and random changes:-- * @Trans a = a -> Dist a ==> c = Dist@-- * @RChange a = a -> R a ==> c = R = IO@---}-class Iterate c where- (*.) :: Int -> (a -> c a) -> (a -> c a)- while :: (a -> Bool) -> (a -> c a) -> (a -> c a)- until :: (a -> Bool) -> (a -> c a) -> (a -> c a)- until p = while (not.p)--infix 8 *.---- iteration of transitions----instance Iterate Dist where- n *. t = head . (n *.. t)- while p t x = if p x then t x >>= while p t else certainly x---- iteration of random changes----instance Iterate IO where- n *. r = (>>= return . head) . rWalk n r- while p t x = do {l <- t x; if p l then while p t l else return l}----{- |-Simulation means to repeat a random chage many times and-to accumulate all results into a distribution. Therefore,-simulation can be regarded as an approximation of distributions-through randomization.--The Sim class allows the overloading of simulation for different-kinds of generators, namely transitions and random changes:-- * @Trans a = a -> Dist a ==> c = Dist@-- * @RChange a = a -> R a ==> c = R = IO@--}-class Sim c where- -- | returns the final randomized transition- (~.) :: Ord a => Int -> (a -> c a) -> RTrans a- -- | returns the whole trace- (~..) :: Ord a => (Int,Int) -> (a -> c a) -> RExpand a- (~*.) :: Ord a => (Int,Int) -> (a -> c a) -> RTrans a--infix 6 ~.-infix 6 ~..---- simulation for transitions----instance Sim Dist where- (~.) x = (~.) x . random- (~..) x = (~..) x . random- (~*.) x = (~*.) x . random----- simulation for random changes----instance Sim IO where- (~.) n t = rDist . replicate n . t- (~..) (k,n) t = mergeTraces . replicate k . rWalk n t- (~*.) (k,n) t = k ~. n *. t--infix 8 ~*.----(~*.) :: (Iterate c,Sim c,Ord a) => (Int,Int) -> (a -> c a) -> RTrans a---(k,n) ~*. t =----- * Tracing--type Trace a = [a]-type Space a = Trace (Dist a)-type Walk a = a -> Trace a-type Expand a = a -> Space a---{- |-@(>>:)@ composes the result of a transition with a space-(transition is composed on the left)--@(a -> m a) -> (a -> [m a]) -> (a -> [m a])@--}-(>>:) :: Trans a -> Expand a -> Expand a-f >>: g = \x -> let ds@(D d:_)=g x in- D [ (z,p*q) | (y,p) <- d, (z,q) <- unD (f y)]:ds--infix 6 >>:---- | walk is a bounded version of the predefined function iterate-walk :: Int -> Change a -> Walk a-walk n f = take n . iterate f--{- |-@(*..)@ is identical to @(*.)@,-but returns the list of all intermediate distributions--}-(*..) :: Int -> Trans a -> Expand a-0 *.. _ = singleton . certainly-1 *.. t = singleton . t-n *.. t = t >>: (n-1) *.. t--infix 8 *..---type RTrace a = R (Trace a)-type RSpace a = R (Space a)-type RWalk a = a -> RTrace a-type RExpand a = a -> RSpace a---- (a -> m a) -> (a -> m [a]) -> (a -> m [a])-composelR :: RChange a -> RWalk a -> RWalk a-composelR f g x = do {rs@(r:_) <- g x; s <- f r; return (s:rs)}---{- |-'rWalk' computes a list of values by randomly selecting-one value from a distribution in each step.--}-rWalk :: Int -> RChange a -> RWalk a-rWalk 0 _ = return . singleton-rWalk 1 t = (>>= return . singleton) . t-rWalk n t = composelR t (rWalk (n-1) t)---{- |-'mergeTraces' converts a list of 'RTrace's-into a list of randomized distributions, i.e., an 'RSpace',-by creating a randomized distribution for each list position across all traces--}-mergeTraces :: Ord a => [RTrace a] -> RSpace a-mergeTraces = fmap (zipListWith (norm . uniform)) . sequence- where- zipListWith :: ([a] -> b) -> [[a]] -> [b]- zipListWith f = map f . transpose--{--for quickCheck--LAWS-- const . pick = random . const---}
− Queuing.hs
@@ -1,144 +0,0 @@-{- |--Model:-- one server serving customers from one queue---}--module Queuing where---import Probability (Dist, Trans, RDist, R, pick, rDist, mapD, )-import List (nub,sort)--type Time = Int---- | (servingTime, nextArrival)-type Profile = (Time, Time)--type Event a = (a,Profile)---- | customers and their individual serving times-type Queue a = [(a,Time)]---- | (customers waiting,validity period of that queue)-type State a = (Queue a,Time)--type System a = [([a],Time)]--type Events a = [Queuing.Event a]---event :: Time -> Events a -> Queue a -> [State a]-event = mEvent 1----event _ [] [] = []---event 0 ((c,(s,a)):es) q = event a es (q++[(c,s)])---event a es [] = ([],a):event 0 es []---event a [] (q@((c,s):q')) = (q,s):event a [] q'---event a es (q@((c,s):q')) | a<s = (q,a):event 0 es ((c,s-a):q')--- | True = (q,s):event (a-s) es q'--system :: Events a -> System a---system es = map (\(q,t)->(map fst q,t)) $ event 0 es []-system = mSystem 1----- | multiple servers--mEvent :: Int -> Time -> Events a -> Queue a -> [State a]-mEvent _ _ [] [] = []-mEvent n 0 ((c,(s,a)):es) q = mEvent n a es (q++[(c,s)])-mEvent n a es [] = ([],a):mEvent n 0 es []-mEvent n _ [] q = (q,s):mEvent n 0 [] (mServe n s q)- where s = mTimeStep n q-mEvent n a es q =- if a < s- then (q,a) : mEvent n 0 es (mServe n a q)- else (q,s) : mEvent n (a-s) es (mServe n s q)- where s = mTimeStep n q----- | decrease served customers remaining time by specified amount-mServe :: Int -> Int -> Queue a -> Queue a-mServe _ _ [] = []-mServe 0 _ x = x-mServe n c ((a,t):es) =- if t > c- then (a,t-c) : mServe (n-1) c es- else mServe (n-1) c es---- | time until next completion-mTimeStep :: Int -> Queue a -> Int-mTimeStep _ ((_,t):[]) = t-mTimeStep 1 ((_,t):_) = t-mTimeStep n ((_,t):es) = min t (mTimeStep (n-1) es)-mTimeStep _ _ = error "Queuing.mTimeStep: queue must be non-empty"--mSystem :: Int -> Events a -> System a-mSystem n es = map (\(q,t)->(map fst q,t)) $ mEvent n 0 es []----- * random--type RProfile = (Dist Time, Trans Time)--type REvent a = (a, RProfile)--type REvents a = [REvent a]--rSystem :: Int -> REvents a -> R (System a)-rSystem n re = do- e <- rBuildEvents re- return (mSystem n e)--rBuildEvents :: REvents a -> R (Events a)-rBuildEvents ((a,(dt,tt)):ex) = do- rest <- rBuildEvents ex- t <- pick dt- nt <- pick $ tt t- return ((a,(t,nt)):rest)-rBuildEvents [] = return []--rmSystem :: Ord a => Int -> Int -> REvents a -> RDist (System a)-rmSystem c n re = rDist $ replicate c (rSystem n re)--evalSystem :: Ord a => Int -> Int -> REvents a -> (System a -> b) -> RDist b-evalSystem c n re ef = do- rds <- rmSystem c n re- return (mapD ef rds)--unit :: b -> ((), b)-unit = (\p->((),p)) -- mapD (\p->((),p))----- * evaluation--maxQueue :: Ord a => System a -> Int-maxQueue s = maximum [length q | (q,_) <- s]--allWaiting :: Ord a => Int -> System a -> [a]-allWaiting n s = nub $ sort $ concat [ drop n q | (q,_) <- s]---countWaiting :: Ord a => Int -> System a -> Int-countWaiting n = length . (allWaiting n)--waiting :: Int -> System a -> Time-waiting n s = sum [ t*length q' | (q,t) <- s, let q' = drop n q]--inSystem :: System a -> Time-inSystem s = sum [ t*length q | (q,t) <- s]--total :: System a -> Time-total = sum . map snd--server :: Int -> System a -> Time-server n s = sum [ t*length q' | (q,t) <- s, let q' = take n q]--idle :: Int -> System a -> Time-idle n s = sum [ t*(n - length q) | (q,t) <- s, length q <= n]--idleAvgP :: Int -> System a -> Float-idleAvgP n s = (fromIntegral $ idle n s) / (fromIntegral $ server n s)
− Show.hs
@@ -1,16 +0,0 @@-module Show where--showL :: Show a => Int -> a -> String-showL n x = s++rep (n-length s) ' '- where s=show x--showR :: Show a => Int -> a -> String-showR n x = rep (n-length s) ' '++s- where s=show x----showP :: Float -> String---showP f = showR 3 (round (f*100))++"%"--rep :: Int -> a -> [a]-rep n x = take n (repeat x)-
ToDo view
@@ -1,14 +1,28 @@+Examples:+ Election and prognoses++generalize fixed Dist to (Distribution prob) whereever possible++Use pretty printer in PrintList? Which one?+ current pretty function is nice for single Dist values,+ but not for (Dist a, Dist b) et.al.+++Collection.draw using StateT++QuickCheck properties+ use a non-empty list structure for the distribution-more efficient data structure,- we will run into the 'monad instance for Data.Set' problem- see http://www.randomhacks.net/articles/2007/03/15/data-set-monad-haskell-macros- it's certainly better to provide a 'collaps' function for removing duplicates-use monad functions instead of custom Dist functions- check where 'collaps' must be applied-separate module name space++ generalize ToFloat class to Num-use pretty printer in PrintList?- current Show instance is nice for single Dist values, but not for (Dist a, Dist b) et.al.-simplify examples (boys, monty hall et.al.)-replace RandomIO by Random and State monad-QuickCheck properties+ Need for multi-parameter type classes?+++create export list++new data type for continuous distributions:++ C (Float -> Float)++prove correctness of >>=? (What did Martin mean with this comment?)
− TreeGrowth.hs
@@ -1,143 +0,0 @@-module TreeGrowth where--import qualified Probability-import Probability- (Dist, R, Space, mapD, normal, unfoldT, certainly, printR,- Trans, RTrans, Expand, RExpand, (*.), (*..), (~..), (~*.), enumPC, )-import Visualize (- Vis, Color(Green, Red, Blue), Plot,- fig, figP, figure, title,- xLabel, yLabel, plotD, color, label,- )---type Height = Int--data Tree = Alive Height | Hit Height | Fallen- deriving (Ord,Eq,Show)--grow :: Trans Tree-grow (Alive h) = normal [Alive k | k <- [h+1..h+5]]-grow _ = error "TreeGrowth.grow: only alive trees can grow"--hit :: Trans Tree-hit (Alive h) = certainly (Hit h)-hit _ = error "TreeGrowth.hit: only alive trees can be hit"--fall :: Trans Tree-fall _ = certainly Fallen--evolve :: Trans Tree-evolve t@(Alive _) = unfoldT (enumPC [90,4,6] [grow,hit,fall]) t-evolve t = certainly t--- evolve t@(Alive _) = unfoldT (enum [0.9,0.04,0.06] [grow,hit,fall]) t--{- |-tree growth simulation:- start with seed and run for n generations--}-seed :: Tree-seed = Alive 0----- * exact results---- | @tree n@ : tree distribution after n generations-tree :: Int -> Tree -> Dist Tree-tree n = n *. evolve---- | @hist n@ : history of tree distributions for n generations-hist :: Int -> Expand Tree-hist n = n *.. evolve----- * simulation results--{- |-Since '(*.)' is overloaded for Trans and RChange,-we can run the simulation ~. directly to @n *. live@.--}----simTree k n = k ~. tree n-simTree :: Int -> Int -> RTrans Tree-simTree k n = (k,n) ~*. evolve--simHist :: Int -> Int -> RExpand Tree-simHist k n = (k,n) ~.. evolve--t2 :: Dist Tree-t2 = tree 2 seed--h2 :: Space Tree-h2 = hist 2 seed--sh2, st2 :: R ()-st2 = printR $ simTree 2000 2 seed-sh2 = printR $ simHist 2000 2 seed----- Alternatives:------ simTree k n = k ~. n *. random evolve--- simTree k n = (k,n) ~*. evolve----- take a trace---height :: Tree -> Int-height Fallen = 0-height (Hit h) = h-height (Alive h) = h-{---myPlot = plotD ((5 *. evolve) (Alive 0) >>= height)--myPlot2 = figP figure{title="Tree Growth",xLabel="Height (m)",- yLabel="Probability"}- (autoColor [- plotD ((5 *. evolve) (Alive 0) >>= height)- ])----}--p1, p2, p3, p4, p5, p6 :: Vis--p1 = fig [plotD $ normal ([1..20]::[Int])]--p2 = fig [plotD $ mapD height (tree 5 seed)]--p3 = figP figure{title="Tree Growth",- xLabel="Height (ft)",- yLabel="Probability"}- [plotD $ mapD height (tree 5 seed)]---p4 = figP figure{title="Tree Growth",- xLabel="Height (ft)",- yLabel="Probability"}- [heightAtTime 5, heightAtTime 10,heightAtTime 15]--heightAtTime :: Int -> Plot-heightAtTime y = plotD $ mapD height (tree y seed)--p5 = figP figure{title="Tree Growth",- xLabel="Height (ft)",- yLabel="Probability"}- (map heightAtTime [3,5,7])--heightCurve :: (Int,Color) -> Plot-heightCurve (n,c) = (heightAtTime n){color=c,label=show n++" Years"}--p6 = figP figure{title="Tree Growth",- xLabel="Height (ft)",- yLabel="Probability"}- (map heightCurve- [(3,Blue),(5,Green),(7,Red)])---done :: Tree -> Bool-done (Alive x) = x >= 5-done _ = True--ev5 :: Tree -> Dist Tree-ev5 = Probability.until done evolve
− Visualize.hs
@@ -1,241 +0,0 @@-module Visualize where--import Probability- (Dist, R, RDist, mapD, unD, norm,- ToFloat, FromFloat, toFloat, fromFloat, )-import PrintList (asTuple, )-import List (nub, sort, sortBy, )---{- TO DO:--* Change function representation in Plot to- xs :: [Float]- ys :: [Float]- and add functions to create this representation from- functions, distributions, and lists- (i.e. plotF, plotD, plotL)---}----- | global settings for one figure----data FigureEnv = FE { fileName :: String,- title :: String,- xLabel :: String,- yLabel :: String }- deriving Show---- | default settings for figure environment----figure :: FigureEnv-figure = FE { fileName = "FuSE.R",- title = "Output",- xLabel = "x",- yLabel = "f(x)" }----- * types to represent settings for individual plots----data Color = Black | Blue | Green | Red | Brown | Gray- | Purple | DarkGray | Cyan | LightGreen | Magenta- | Orange | Yellow | White | Custom Int Int Int- deriving Eq--instance Show Color where- show Black = "\"black\""- show Blue = "\"blue\""- show Green = "\"green\""- show Red = "\"red\""- show Brown = "\"brown\""- show Gray = "\"gray\""- show Purple = "\"purple\""- show DarkGray = "\"darkgray\""- show Cyan = "\"cyan\""- show LightGreen = "\"lightgreen\""- show Magenta = "\"magenta\""- show Orange = "\"orange\""- show Yellow = "\"yellow\""- show White = "\"white\""- show (Custom r g b) = "rgb("++(show r)++", "++(show g)++", "++(show b)++")"--data LineStyle = Solid | Dashed | Dotted | DotDash | LongDash | TwoDash- deriving Eq--instance Show LineStyle where- show Solid = "1"- show Dashed = "2"- show Dotted = "3"- show DotDash = "4"- show LongDash = "5"- show TwoDash = "6"--type PlotFun = Float -> Float----- | settings for individual plots----data Plot = Plot { ys :: [Float],- xs :: [Float],- color :: Color,- lineStyle :: LineStyle,- lineWidth :: Int,- label :: String }--{--instance Show Plot where- show _ = "Individual plots cannot be printed.\nPlease use plots \- \ as arguments to the fig function."--}----- | default plotting environment----plot :: Plot-plot = Plot { ys = [0],- xs = [0],- color = Black,- lineStyle = Solid,- lineWidth = 1,- label = "" }--colors :: [Color]-colors = [Blue,Green,Red,Purple,Black,Orange,Brown,Yellow]--setColor :: Plot -> Color -> Plot-setColor p c = p{color=c}--autoColor :: [Plot] -> [Plot]-autoColor ps | length ps <= n = zipWith setColor ps colors- | otherwise = error ("autoColor works for no more than "++- show n++" plots.")- where n=length colors---- | create a plot from a distribution----plotD :: ToFloat a => Dist a -> Plot---plotD d = plot{ys = map (\x->(dp $ prob' x d')) (extract d'),--- xs = extract d'}-plotD d = plot{xs = tfl, ys = pdl}- where d' = mapD toFloat d- d'' = norm d'- pl = unD d''- pl' = sortBy (\(a,_) (a',_) -> compare (toFloat a) (toFloat a')) pl- (tfl, pdl) = unzip pl'- -- dp (P p) = p- -- pl'' = map dp pdl--plotRD :: ToFloat a => RDist a -> IO Plot-plotRD a = fmap plotD a---- | create a plot from a function----plotF :: (FromFloat a,ToFloat b) => (Float,Float,Float) -> (a -> b) -> Plot-plotF xd g = plot{ys = map (\x->toFloat (g (fromFloat x))) (xvals xd),xs = xvals xd}- where xvals (a,b,d) =- if a > b then [] else a:xvals (a+d,b,d)---- | create a plot from a list----plotL :: ToFloat a => [a] -> Plot-plotL vs = plot{ys = map toFloat vs, xs = map toFloat [1..length vs]}---plotRL :: ToFloat a => R [a] -> IO Plot-plotRL a = fmap plotL a-----yls :: ToFloat a => [a] -> [Plot] -> [[Float]]---yls xs (p:ps) = [f p (toFloat v) | v <- xs ]:yls xs ps---yls _ [] = []--yls :: [Float] -> Plot -> Plot-yls xl p = p{xs=x', ys=y'}- where t = zip (xs p) (ys p)- t' = metaTuple xl t- (x', y') = unzip t'--metaTuple :: [Float] -> [(Float,Float)] -> [(Float,Float)]-metaTuple (x:xl) ((p,v):px) | p == x = (p,v):(metaTuple xl px)-metaTuple (x:xl) p'@( (p,_):_ ) | p > x = (x,0):(metaTuple xl p')-metaTuple x [] = map (\v->(v,0)) x-metaTuple x y = error $ (show x)++(show y)---- | we want to increase the bounds absolutely, account for negative numbers----incr, decr :: (Ord a, Fractional a) => a -> a-incr x =- if x > 0- then x * 1.05- else x * 0.95--decr x =- if x > 0- then x * 0.95- else x * 1.05---- | Visualization output----type Vis = IO ()----- * creating figures----fig :: [Plot] -> Vis-fig = figP figure--figP :: FigureEnv -> [Plot] -> Vis-figP fe ps = do let xl = sort $ nub $ concatMap xs ps- let minx = minimum xl--- let maxx = maximum xl- let n = length xl- let ys' = map ys (map (yls xl) ps) -- yls xl ps- let miny = minimum (map minimum ys')- let maxy = maximum (map maximum ys')- let out0' = out0 (fileName fe)- let out1' = out1 (fileName fe)- out0' ("x <- "++(vec xl))- out1' ("y <- "++(vec $ (decr miny):(replicate (n-1) (incr maxy))))- out1' ("plot(x,y,type=\"n\",main=\""++- 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\")")---{--define:- * autoLabel- * showParams--}--showParams :: Show a => [a] -> [String] -> String-showParams xs0 ss =- asTuple id (zipWith (\x s-> show x++":"++s) xs0 ss)--legend :: Float -> Float -> [Plot] -> String-legend x y ps = "legend("++(show x)++", "++(show y)++","++- "lty="++vec (map lineStyle ps)++","++- "col="++vec (map color ps)++","++- "lwd="++vec (map lineWidth ps)++","++- "legend="++vec (map label ps)++")"--drawy :: ToFloat a => Int -> Plot -> [a] -> String-drawy yn p fl = "y"++(show yn)++" <- "++(vec (map toFloat fl))++"\n"++- "lines(x,y"++(show yn)++",col="++(show $ color p)++","++- "lty="++(show $ lineStyle p)++",lwd="++(show $ lineWidth p)++")"---vec :: Show a => [a] -> String-vec xs0 = "c"++asTuple show xs0--out0 :: String -> String -> IO ()-out0 f s = writeFile (f) (s++"\n")--out1 :: String -> String -> IO ()-out1 f s = appendFile (f) (s++"\n")
probability.cabal view
@@ -1,12 +1,12 @@ Name: probability-Version: 0.1+Version: 0.2 License: BSD3-Author: Martin Erwig <erwig@eecs.oregonstate.edu>+Author: Martin Erwig <erwig@eecs.oregonstate.edu>, Steve Kollmansberger Maintainer: Henning Thielemann <haskell@henning-thielemann.de> Homepage: http://darcs.haskell.org/probability Category: Math, Monads, Graphics-Build-Depends: base, haskell98-Synopsis: Computations with discrete random variables+Build-Depends: base, mtl+Synopsis: Probabilistic Functional Programming Description: The Library allows exact computation with discrete random variables in terms of their distributions by using a monad.@@ -16,26 +16,36 @@ Tested-With: GHC==6.4 Build-Type: Simple License-File: COPYRIGHT-Hs-Source-Dirs: .+Hs-Source-Dirs: src+GHC-Options: -Wall Exposed-Modules:- Alarm- Barber- Bayesian- Boys- Collection- Dice- MontyHall- NBoys- Predator- Probability- Queuing- TreeGrowth- Visualize+ Numeric.Probability.Visualize+ Numeric.Probability.Expectation+ Numeric.Probability.Percentage+ Numeric.Probability.Distribution+ Numeric.Probability.Transition+ Numeric.Probability.Random+ Numeric.Probability.Shape+ Numeric.Probability.Trace+ Numeric.Probability.Simulation+ Numeric.Probability.Object+ Numeric.Probability.Example.Alarm+ Numeric.Probability.Example.Barber+ Numeric.Probability.Example.Bayesian+ Numeric.Probability.Example.Boys+ Numeric.Probability.Example.Collection+ Numeric.Probability.Example.Diagnosis+ Numeric.Probability.Example.Dice+ Numeric.Probability.Example.DiceAccum+ Numeric.Probability.Example.MontyHall+ Numeric.Probability.Example.NBoys+ Numeric.Probability.Example.Predator+ Numeric.Probability.Example.Queuing+ Numeric.Probability.Example.TreeGrowth Other-Modules:- ListUtils- PrintList- Show+ Numeric.Probability.Monad+ Numeric.Probability.PrintList+ Numeric.Probability.Show Extra-Source-Files: README ToDo-GHC-Options: -Wall -O2
+ src/Numeric/Probability/Distribution.hs view
@@ -0,0 +1,342 @@+-- | Deterministic and probabilistic values++module Numeric.Probability.Distribution where++import Numeric.Probability.Show (showR)+import qualified Numeric.Probability.Shape as Shape++import Control.Monad (liftM, liftM2, join, )++import qualified Data.Map as Map+import qualified Data.List as List++import Prelude hiding (map, filter)+++-- * Events+type Event a = a -> Bool++oneOf :: Eq a => [a] -> Event a+oneOf = flip elem++just :: Eq a => a -> Event a+just = (==)++++-- * Distributions++{- |+Probability disribution++The underlying data structure is a list.+Unfortunately we cannot use a more efficient data structure+because the key type must be of class 'Ord',+but the 'Monad' class does not allow constraints for result types.+The Monad instance is particularly useful+because many generic monad functions make sense here,+monad transformers can be used+and the monadic design allows to simulate probabilistic games in an elegant manner.++We have the same problem like making "Data.Set" an instance of 'Monad',+see <http://www.randomhacks.net/articles/2007/03/15/data-set-monad-haskell-macros>++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.+-}+newtype T prob a = Cons {decons :: [(a,prob)]}++certainly :: Num prob => a -> T prob a+certainly x = Cons [(x,1)]++instance Num prob => Monad (T prob) where+ return = certainly+ d >>= f = Cons [(y,q*p) | (x,p) <- decons d, (y,q) <- decons (f x)]+ fail _ = Cons []++{-+Dist cannot be an instance of MonadPlus,+because there is no mzero+(it would be an empty list of events, but their probabilities do not sum up to 1)+and thus it breaks the normalization for the >>= combinator.+See for instance the Boys example:++ do f <- family+ guard (existsBoy f)+ return f++mplus is not associative because we have to normalize the sum of probabilities to 1.++instance MonadPlus Dist where+ mzero = Cons []+ mplus d d' =+ if isZero d || isZero d'+ then mzero+ else unfoldD $ choose 0.5 d d'++isZero :: Dist a -> Bool+isZero (Cons d) = null d+-}+++instance Functor (T prob) where+ fmap f (Cons d) = Cons [(f x,p) | (x,p) <- d]++++errorMargin :: RealFloat prob => prob+errorMargin =+ let eps = 10 * fromInteger (floatRadix eps) ^ (- floatDigits eps)+ in eps++{- |+Check whether two distributions are equal when neglecting rounding errors.+We do not want to put this into an 'Eq' instance,+since it is not exact equivalence+and it seems to be too easy to mix it up with @liftM2 (==) x y@.+-}+approx :: (RealFloat prob, Ord a) =>+ T prob a -> T prob a ->+ Bool+approx (Cons xs) (Cons ys) =+ let (xse, xsp) = unzip (norm' xs)+ (yse, ysp) = unzip (norm' ys)+ in xse == yse &&+ all (\p -> abs p < errorMargin) (zipWith (-) xsp ysp)+++-- ** Auxiliary functions for constructing and working with distributions+lift :: (Num prob) =>+ ([(a,prob)] -> [(a,prob)]) ->+ T prob a -> T prob a+lift f = Cons . f . decons++size :: T prob a -> Int+size = length . decons++check :: RealFloat 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 = 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))++sortElem :: Ord a => [(a,prob)] -> [(a,prob)]+sortElem = List.sortBy (\x y -> compare (fst y) (fst x))+++-- ** Normalization = grouping+norm :: (Num prob, Ord a) => T prob a -> T prob a+norm = lift norm'++norm' :: (Num prob, Ord a) => [(a,prob)] -> [(a,prob)]+norm' =+ Map.toAscList . Map.fromListWith (+)++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+++-- | pretty printing+pretty :: (Ord a, Show a, Num prob, Ord prob) =>+ (prob -> String) -> T prob a -> String+pretty _ (Cons []) = "Impossible"+pretty showProb (Cons xs) =+ let w = maximum (List.map (length.show.fst) xs)+ in concatMap+ (\(x,p) -> showR w x++' ': showProb p++"\n")+ (sortP (norm' xs))++infix 0 //%++(//%) :: (Ord a, Show a) => T Rational a -> () -> IO ()+(//%) x () = putStr (pretty show x)++instance (Num prob, Ord prob, Ord a, Show a) =>+ Show (T prob a) where+ showsPrec p (Cons xs) =+ showParen (p>10)+ (showString "fromFreqs " . showsPrec 11 (sortP (norm' xs)))++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."++{-+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.+This is essential for performance.++Thus @sum $ replicate 10 d@ is significantly faster+than @fmap sum $ replicateM 10 d@+-}+instance (Num prob, Ord prob, Ord a, Num a) => Num (T prob a) where+ fromInteger = return . fromInteger+ x + y = norm (liftM2 (+) x y)+ x - y = norm (liftM2 (-) x y)+ x * y = norm (liftM2 (*) x y)+ abs x = norm (liftM abs x)+ signum x = norm (liftM signum x)+ negate x = liftM negate x++instance (Num prob, Ord prob, Ord a, Fractional a) =>+ Fractional (T prob a) where+ fromRational = return . fromRational+ recip x = liftM recip x+ x / y = norm (liftM2 (/) x y)++++-- * Spread: functions to convert a list of values into a distribution++-- | distribution generators+type Spread prob a = [a] -> T prob a++{- not a valid distribution+impossible :: T prob a+impossible = mzero+-}++choose :: Num prob => prob -> a -> a -> T prob a+choose p x y = Cons $ zip [x,y] [p,1-p]++enum :: Fractional prob => [Int] -> Spread prob a+enum = relative . List.map fromIntegral++{- |+Give a list of frequencies, they do not need to sum up to 1.+-}+relative :: Fractional prob => [prob] -> Spread prob a+relative ns = fromFreqs . flip zip ns++shape :: Fractional prob =>+ (prob -> prob) -> Spread prob a+shape _ [] = error "Probability.shape: empty list"+shape f xs =+ let incr = 1 / fromIntegral (length xs - 1)+ ps = List.map f (iterate (+incr) 0)+ in fromFreqs (zip xs ps)++linear :: Fractional prob => Spread prob a+linear = shape Shape.linear++uniform :: Fractional prob => Spread prob a+uniform = shape Shape.uniform++negExp :: Floating prob => Spread prob a+negExp = shape Shape.negExp++normal :: Floating prob => Spread prob a+normal = shape Shape.normal++++-- | extracting and mapping the domain of a distribution+extract :: T prob a -> [a]+extract = List.map fst . decons++-- | 'fmap' with normalization+map :: (Num prob, Ord b) =>+ (a -> b) -> T prob a -> T prob b+map f = norm . fmap f+++{- |+unfold a distribution of distributions into one distribution,+this is 'Control.Monad.join' with normalization.+-}+unfold :: (Num prob, Ord a) =>+ T prob (T prob a) -> T prob a+unfold = norm . join+++-- | conditional distribution+cond :: (Num prob) =>+ T prob Bool ->+ T prob a {-^ True -} ->+ T prob a {-^ False -} ->+ T prob a+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"+++infixl 1 >>=?+infixr 1 ?=<<++-- | conditional probability, identical to 'Dist.filter'+(?=<<) :: (Fractional prob) =>+ (a -> Bool) -> T prob a -> T prob a+(?=<<) = filter++{- |+'Dist.filter' in infix form.+Can be considered an additional monadic combinator,+which can be used where you would want 'Control.Monad.guard' otherwise.+-}+(>>=?) :: (Fractional prob) =>+ T prob a -> (a -> Bool) -> T prob a+(>>=?) = flip filter+++-- | filtering distributions+data Select a = Case a | Other+ deriving (Eq,Ord,Show)++above :: (Num prob, Ord prob, Ord a) =>+ prob -> T prob a -> T prob (Select a)+above p (Cons d) =+ let (d1,d2) = span (\(_,q)->q>=p) (sortP (norm' d))+ in Cons (List.map (\(x,q)->(Case x,q)) d1++[(Other,sumP d2)])++fromFreqs :: (Fractional prob) => [(a,prob)] -> T prob a+fromFreqs xs = Cons (List.map (\(x,p)->(x,p/q)) xs)+ where q = sumP xs++filter :: (Fractional prob) =>+ (a -> Bool) -> T prob a -> T prob a+filter p = fromFreqs . List.filter (p . fst) . decons+++-- | selecting from distributions+selectP :: (Num prob, Ord prob) => T prob a -> prob -> a+selectP (Cons d) p = scanP p d++scanP :: (Num prob, Ord prob) => prob -> [(a,prob)] -> a+scanP p ((x,q):ps) =+ if p<=q || null ps+ then x+ else scanP (p-q) ps+scanP _ [] = error "Probability.scanP: distribution must be non-empty"++infixr 1 ??++(??) :: Num prob => Event a -> T prob a -> prob+(??) p = sumP . List.filter (p . fst) . decons+++-- | expectation value+expected :: (Num a) => T a a -> a+expected = sum . List.map (\(x,p) -> x * p) . decons++-- | statistical analyses+variance :: (Num a) => T a a -> a+variance x =+ expected (fmap ((^(2::Int)) . subtract (expected x)) x)++stdDev :: (Floating a) => T a a -> a+stdDev = sqrt . variance
+ src/Numeric/Probability/Example/Alarm.hs view
@@ -0,0 +1,62 @@+module Numeric.Probability.Example.Alarm where++import qualified Numeric.Probability.Distribution as Dist+import Numeric.Probability.Distribution ((??), (?=<<), )+++type Probability = Rational+type Dist a = Dist.T Probability a+type PBool = Dist.T Probability Bool+++flp :: Probability -> PBool+flp p = Dist.choose p True False+++-- * Numeric.Probability.Example.Alarm network++-- | prior burglary 1%+b :: PBool+b = flp 0.01++-- | prior earthquake 0.1%+e :: PBool+e = flp 0.001++-- | conditional probability of alarm given burglary and earthquake+a :: Bool -> Bool -> PBool+a b0 e0 =+ case (b0,e0) of+ (False,False) -> flp 0.01+ (False,True) -> flp 0.1+ (True,False) -> flp 0.7+ (True,True) -> flp 0.8+++-- | conditional probability of john calling given alarm+j :: Bool -> PBool+j a0 = if a0 then flp 0.8 else flp 0.05++-- | conditional probability of mary calling given alarm+m :: Bool -> PBool+m a0 = if a0 then flp 0.9 else flp 0.1++-- | calculate the full joint distribution+data Burglary = B { burglary :: Bool,+ earthquake :: Bool,+ alarm :: Bool,+ john :: Bool,+ mary :: Bool }+ deriving (Eq, Ord, Show)++bJoint :: Dist Burglary+bJoint = do b' <- b -- burglary+ e' <- e -- earthquake+ a' <- a b' e' -- alarm+ j' <- j a' -- john+ m' <- m a' -- mary+ return (B b' e' a' j' m')++-- | what is the probability that mary calls given that john calls?+pmj :: Probability+pmj = mary ?? john ?=<< bJoint
+ src/Numeric/Probability/Example/Barber.hs view
@@ -0,0 +1,69 @@+module Numeric.Probability.Example.Barber where++import qualified Numeric.Probability.Distribution as Dist+import Numeric.Probability.Example.Queuing+ (Time, System, unit, evalSystem, idleAvgP, waiting)++import Numeric.Probability.Percentage+ (Dist, RDist, Trans, )++{- no Random instance for Rational+type Probability = Rational+type Dist a = Dist.T Probability a+type RDist a = Rnd.Distribution Probability a+type Trans a = Transition Probability a+-}+++-- * barber shop++custServ :: Dist Time+custServ = Dist.normal [5..10]++nextCust :: Trans Time -- not dependant on serving time+nextCust _ = Dist.normal [3..6]++barbers :: Int+barbers = 1++customers :: Int+customers = 20++runs :: Int+runs = 50++barberEvent :: ((), (Dist Time, Time -> Dist Time))+barberEvent = unit (custServ, nextCust)++barberEvents :: [((), (Dist Time, Time -> Dist Time))]+barberEvents = replicate customers barberEvent++barberSystem :: (Ord b) => (System () -> b) -> RDist b+barberSystem eval = evalSystem runs barbers barberEvents eval+++-- * category++data Category = ThreeOrLess | FourToTen | MoreThanTen+ deriving (Eq,Ord,Show)++cat :: Time -> Category+cat n | n <= 3 = ThreeOrLess+cat n | n <= 10 = FourToTen+cat _ = MoreThanTen++perc :: Float -> String+perc n | n <= 0.25 = "0% to 25%"+perc n | n <= 0.5 = "25% to 50%"+perc n | n <= 0.75 = "50% to 75%"+perc _ = "75% to 100%"++-- * evaluation++-- | avg barber idle time+barberIdle :: RDist String+barberIdle = barberSystem (perc . idleAvgP barbers)++-- | avg customer waiting time (unserved customers)+customerWait :: RDist Category+customerWait = barberSystem (cat . (`div` customers) . waiting barbers)
+ src/Numeric/Probability/Example/Bayesian.hs view
@@ -0,0 +1,101 @@+{- |++Approach: model a node with k predecessors as a function with k+ parameters++-}+module Numeric.Probability.Example.Bayesian where++import qualified Numeric.Probability.Distribution as Dist+import qualified Numeric.Probability.Transition as Trans+import qualified Numeric.Probability.Monad as MonadExt+import Numeric.Probability.Distribution ((??), (?=<<), )++++-- * Abbreviations, smart constructors++type Probability = Rational+type Dist a = Dist.T Probability a++type State a = [a]+type PState a = Dist (State a)+type STrans a = State a -> PState a+type SPred a = a -> State a -> Bool++event :: Probability -> a -> STrans a+event p e0 = Trans.maybe p (e0:)++happens :: Eq a => SPred a+happens = elem++network :: [STrans a] -> PState a+network = flip MonadExt.compose []+++source :: Probability -> a -> STrans a+source = event++bin :: Eq a =>+ a -> a -> Probability -> Probability -> Probability -> Probability ->+ a -> STrans a+bin x y a b c d z s | elem x s && elem y s = event a z s+ | elem x s = event b z s+ | elem y s = event c z s+ | otherwise = event d z s+++-- | Two possible causes for one effect++data Nodes = A | B | E deriving (Eq,Ord,Show)++g :: PState Nodes+g = network [source 0.1 A,+ source 0.2 B,+ bin A B 1 1 0.5 0 E]++-- * queries++e, aE, bE :: Probability+e = happens E ?? g+aE = happens A ?? happens E ?=<< g+bE = happens B ?? happens E ?=<< g+++{-+data State = State {causeA :: Bool, causeB :: Bool, effect :: Bool}+ deriving (Eq,Ord,Show)++nCauseA s = s{causeA=True}+-}++--+-- Wet grass example+--+-- cloudy = true 0.5+--+-- sprinkler c = dep c 0.1 0.5+--+-- rain c = dep c 0.8 0.2+--+-- wetGrass s r = bin s r 0.99 0.9 0.9 0+--+-- c = cloudy+-- s = sprinkler cloudy+-- r = rain cloudy+-- w = wetGrass s r+++-- alarm :: Prob -> Prob -> Prob+-- alarm b e = cond b (pTrue 0.8)+-- (cond e (pTrue 0.1) (pTrue 0.01))+--+-- john :: Prob -> Prob+-- john a = cond a (pTrue 0.7) (pTrue 0.1)+--+-- mary :: Prob -> Prob+-- mary a = cond a (pTrue 0.6) (pTrue 0.2)+--+--+-- maryWhenJohn = mary a ?? john a+-- where a = alarm (pTrue 0.5) (pTrue 0.1)
+ src/Numeric/Probability/Example/Boys.hs view
@@ -0,0 +1,55 @@+{- |+Consider a family of two children. Given that there is a boy in the family,+what is the probability that there are two boys in the family?+-}++module Numeric.Probability.Example.Boys where++import qualified Numeric.Probability.Distribution as Dist+import Numeric.Probability.Distribution ((??), (?=<<), )++import Control.Monad (liftM2, )+++type Probability = Rational+type Dist a = Dist.T Probability a++data Child = Boy | Girl+ deriving (Eq,Ord,Show)++type Family = (Child, Child)++birth :: Dist Child+birth = Dist.uniform [Boy, Girl]++family :: Dist Family+family = liftM2 (,) birth birth++allBoys :: Dist.Event Family+allBoys (c0, c1) = (c0 == Boy && c1 == Boy)++existsBoy :: Dist.Event Family+existsBoy (c0, c1) = (c0 == Boy || c1 == Boy)++familyWithBoy :: Dist Family+familyWithBoy = existsBoy ?=<< family+{-+familyWithBoy =+ do f <- family+ guard (existsBoy f)+ return f+-}++twoBoys :: Probability+twoBoys = allBoys ?? familyWithBoy+++countBoy :: Child -> Int+countBoy Boy = 1+countBoy Girl = 0++countBoys :: Family -> Int+countBoys (c0,c1) = countBoy c0 + countBoy c1++numBoys :: Dist Int+numBoys = Dist.map countBoys familyWithBoy
+ src/Numeric/Probability/Example/Collection.hs view
@@ -0,0 +1,133 @@+module Numeric.Probability.Example.Collection where++import qualified Numeric.Probability.Distribution as Dist+import qualified Numeric.Probability.Random as Rnd+import Numeric.Probability.Distribution ((??), )+import Numeric.Probability.Simulation ((~.), )++import Numeric.Probability.Percentage (Dist)++import Numeric.Probability.Monad (doWhile, )+import Control.Monad.State+ (StateT(StateT, runStateT), evalStateT, liftM2, replicateM)++import qualified Data.List as List+import System.Random (Random)++++type Collection a = [a]++type Probability = Rational+++{- |+see also the proposal+ <http://www.haskell.org/pipermail/libraries/2008-February/009270.html>+-}+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)++select1 :: (Fractional prob) => Collection a -> Dist.T prob a+select1 = evalStateT selectOne++select2 :: (Fractional prob) => Collection a -> Dist.T prob (a,a)+select2 = evalStateT (liftM2 (,) selectOne selectOne)++select :: (Fractional prob) => Int -> Collection a -> Dist.T prob [a]+select n = evalStateT (replicateM n selectOne)+++-- * Example collections++-- ** marbles++data Marble = R | G | B deriving (Eq,Ord,Show)++bucket :: Collection Marble+bucket = [R,R,R,R,R, G,G,G, B,B]++jar :: Collection Marble+jar = [R,R,G,G,B]++-- pRGB = prob (just [R,G,B]) (select 3 bucket)+pRGB :: Probability+pRGB = Dist.just [R,G,B] ?? select 3 jar+pRG :: Probability+pRG = Dist.oneOf [[R,G],[G,R]] ?? select 2 jar++-- ** cards++data Suit = Club | Spade | Heart | Diamond+ deriving (Eq,Ord,Show,Enum)++data Rank = Plain Int | Jack | Queen | King | Ace+ deriving (Eq,Ord,Show)++type Card = (Rank,Suit)++plains :: [Rank]+plains = map Plain [2..10]++faces :: [Rank]+faces = [Jack,Queen,King,Ace]++isFace :: Card -> Bool+isFace (r,_) = r `elem` faces+-- isFace = (`elem` faces) . fst++isPlain :: Card -> Bool+isPlain (r,_) = r `elem` plains++ranks :: [Rank]+ranks = plains ++ faces++suits :: [Suit]+suits = [Club,Spade,Heart,Diamond]++deck :: Collection Card+deck = liftM2 (,) ranks suits+++-- * Example++{- | mini-blackjack:+draw 2 cards, and if value is less than 14, continue drawing+until value equals or exceeds 14. if values exceeds 21,+you lose, otherwise you win.+-}++value :: Card -> Int+value ((Plain n),_) = n+value (Ace,_) = 11+value _ = 10++totalValue :: Collection Card -> Int+totalValue cards = sum (map value cards)++-- this can be made with StateT, too, I think+draw :: (Fractional prob) =>+ ([Card], Collection Card) -> Dist.T prob ([Card], Collection Card)+draw (cards,cl) =+ runStateT (fmap (:cards) selectOne) cl++drawF :: ([Card], Collection Card) -> Dist ([Card], Collection Card)+drawF = draw+++drawTo16 :: Rnd.T ([Card], Collection Card)+drawTo16 =+ doWhile+ (\(cards,_) -> totalValue cards < 16)+ (Rnd.change drawF) ([], deck)++win :: ([Card], b) -> Bool+win (cards,_) = totalValue cards <= 21++chanceWin :: (Fractional prob, Ord prob, Random prob) =>+ Rnd.T (Dist.T prob Bool)+chanceWin = fmap (Dist.map win) ((100 ~. const drawTo16) undefined)
+ src/Numeric/Probability/Example/Diagnosis.hs view
@@ -0,0 +1,60 @@+{- |+You take part in a screening test for a disease+that you have with a probability 'pDisease'.+The test can fail in two ways:+If you are ill,+the test says with probability 'pFalseNegative' that you are healthy.+If you are healthy,+it says with probability 'pFalsePositive' that you are ill.++Now consider the test is positive -+what is the probability that you are indeed ill?+-}+module Numeric.Probability.Example.Diagnosis where++import qualified Numeric.Probability.Distribution as Dist+import Numeric.Probability.Distribution ((??), (?=<<), )+++type Probability = Rational+type Dist a = Dist.T Probability a+++data State = Healthy | Ill+ deriving (Eq, Ord, Show, Enum)++data Finding = Negative | Positive+ deriving (Eq, Ord, Show, Enum)+++pDisease, pFalseNegative, pFalsePositive :: Probability+pDisease = 0.001+pFalseNegative = 0.01+pFalsePositive = 0.01+++dist :: Dist (State, Finding)+dist =+ do s <- Dist.choose pDisease Ill Healthy+ f <- case s of+ Ill -> Dist.choose pFalseNegative Negative Positive+ Healthy -> Dist.choose pFalsePositive Positive Negative+ return (s,f)+++{- |+Alternative way for computing the distribution.+It is usually more efficient because we do not need to switch on the healthy state.+-}+distAlt :: Dist (State, Finding)+distAlt =+ do (s,fr) <-+ Dist.choose pDisease+ (Ill, Dist.choose pFalseNegative Negative Positive)+ (Healthy, Dist.choose pFalsePositive Positive Negative)+ f <- fr+ return (s,f)+++p :: Probability+p = (Dist.just Ill . fst) ?? (Dist.just Positive . snd) ?=<< dist
+ src/Numeric/Probability/Example/Dice.hs view
@@ -0,0 +1,43 @@+module Numeric.Probability.Example.Dice where++import qualified Numeric.Probability.Distribution as Dist+import Numeric.Probability.Distribution ((??), )+import Control.Monad (liftM2, replicateM)+++type Die = Int++type Probability = Rational+type Dist = Dist.T Probability++die :: Dist Die+die = Dist.uniform [1..6]++-- | product of independent distributions+twoDice :: Dist (Die,Die)+twoDice = liftM2 (,) die die++dice :: Int -> Dist [Die]+dice = flip replicateM die+++twoSixes :: Probability+twoSixes = (==(6,6)) ?? twoDice++{- |+@sixes p n@ computes the probability of getting+p sixes (@>1@, @==2@, ...) when rolling n dice+-}+sixes :: (Int -> Bool) -> Int -> Probability+sixes p n = (p . length . filter (==6)) ?? dice n++droll :: Dist Die+droll =+ liftM2 (+) (Dist.uniform [0,1]) die++g3 :: Probability+g3 = (>3) ?? die++addTwo :: Dist Die+addTwo =+ liftM2 (+) die die
+ src/Numeric/Probability/Example/DiceAccum.hs view
@@ -0,0 +1,77 @@+{- |+We play the following game:+We roll a die until we stop or we get three spots.+In the first case we own all spots obtained so far,+in the latter case we own nothing.++What is the strategy for maximizing the expected score?+-}+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+import qualified Numeric.Probability.Monad as MonadExt+import Numeric.Probability.Trace (Trace)++import Numeric.Probability.Example.Dice (Die, )+++type Score = Int+++die :: Fractional prob => Dist.T prob Die+die = Dist.uniform [1..6]++roll :: Fractional prob => Trans.T prob (Maybe Score)+roll =+ maybe+ (return Nothing)+ (\score -> flip fmap die $+ \spots ->+ -- where is my beloved 'toMaybe' ?+ if spots == 3+ then Nothing+ else Just (score + spots))++continue :: Score -> Bool+continue scoreInt =+ let score = fromIntegral scoreInt :: Rational+ in Dist.expected+ (Dist.uniform (0 : map (score+) [1,2,4,5,6])) > score++-- | optimal strategy+strategy :: Fractional prob => Trans.T prob (Maybe Score)+strategy s0 =+ maybe+ (return Nothing)+ (\score ->+ if continue score+ then roll s0+ else return s0) s0++-- | distribution of the scores that are achieved with the optimal strategy+game :: Fractional prob => Dist.T prob (Maybe Score)+game =+ Trans.compose (replicate 18 strategy) (Just 0)+ -- MonadExt.compose (replicate 8 turn) (Just 0)+++{- too inefficient+game :: Fractional prob => Dist.T prob Score+game =+ let turn score =+ if continue score+ then roll score >>= \s -> if s==0 then return 0 else turn s+ else return score+ in turn 0+-}+++walk :: Int -> IO (Trace (Maybe Score))+walk n =+ Rnd.run $+ MonadExt.walk n+ (Rnd.change (roll :: Trans.T Double (Maybe Score)))+ (Just 0)
+ src/Numeric/Probability/Example/MontyHall.hs view
@@ -0,0 +1,112 @@+module Numeric.Probability.Example.MontyHall where++import qualified Numeric.Probability.Distribution as Dist+import qualified Numeric.Probability.Transition as Trans+import Numeric.Probability.Simulation ((~.), )++import Numeric.Probability.Percentage+ (Dist, RDist, Trans, )++import qualified Numeric.Probability.Monad as MonadExt++import Data.List ( (\\) )+++{- no Random instance for Rational+type Probability = Rational+type Dist a = Dist.T Probability a+type RDist a = Rnd.Distribution Probability a+type Trans a = Transition Probability a+-}++data Door = A | B | C+ deriving (Eq,Ord,Show)++doors :: [Door]+doors = [A,B,C]++data State = Doors {prize :: Door, chosen :: Door, opened :: Door}+ deriving (Eq,Ord,Show)+++-- | initial configuration of the game status+start :: State+start = Doors {prize=u,chosen=u,opened=u} where u=undefined+++{- |+Steps of the game:++ (1) hide the prize++ (2) choose a door++ (3) open a non-open door, not revealing the prize++ (4) apply strategy: switch or stay+-}+hide :: Trans State+hide s = Dist.uniform [s {prize = d} | d <- doors]++choose :: Trans State+choose s = Dist.uniform [s {chosen = d} | d <- doors]++open :: Trans State+open s = Dist.uniform [s {opened = d} | d <- doors \\ [prize s,chosen s]]++type Strategy = Trans State++switch :: Strategy+switch s = Dist.uniform [s {chosen = d} | d <- doors \\ [chosen s,opened s]]++stay :: Strategy+stay = Trans.id++game :: Strategy -> Trans State+game s = MonadExt.compose [hide,choose,open,s]+++-- * Playing the game++data Outcome = Win | Lose+ deriving (Eq,Ord,Show)++result :: State -> Outcome+result s = if chosen s==prize s then Win else Lose++eval :: Strategy -> Dist Outcome+eval s = Dist.map result (game s start)++simEval :: Int -> Strategy -> RDist Outcome+simEval k s = Dist.map result `fmap` (k ~. game s) start+++-- * Alternative modeling++firstChoice :: Dist Outcome+firstChoice = Dist.uniform [Win,Lose,Lose]++switch' :: Trans Outcome+switch' Win = Dist.certainly Lose+switch' Lose = Dist.certainly Win+++-- * Play the game the monadic way++type StrategyM = Door -> Door -> Door++stayM :: StrategyM+stayM chosenDoor _openedDoor = chosenDoor++switchM :: StrategyM+switchM chosenDoor openedDoor =+ let [finalDoor] = doors \\ [chosenDoor, openedDoor]+ in finalDoor++evalM :: StrategyM -> Dist Outcome+evalM chooseFinalDoor =+ do prizeDoor <- Dist.uniform doors+ chosenDoor <- Dist.uniform doors+ openedDoor <- Dist.uniform (doors \\ [prizeDoor, chosenDoor])+ return (if chooseFinalDoor chosenDoor openedDoor == prizeDoor+ then Win else Lose)
+ src/Numeric/Probability/Example/NBoys.hs view
@@ -0,0 +1,41 @@+{- |+Ceneralization of "Numeric.Probability.Example.Boys"++Consider a family of n children. Given that there are k boys in the family,+what is the probability that there are m boys in the family?+-}++module Numeric.Probability.Example.NBoys where++import qualified Numeric.Probability.Distribution as Dist+import Numeric.Probability.Distribution (Event, (??), (?=<<), )++import Numeric.Probability.Example.Boys+ (Dist, Probability, Child(Boy), birth, )++import Control.Monad (replicateM)+++type Family = [Child]++family :: Int -> Dist Family+family n = replicateM n birth++countBoys :: Family -> Int+countBoys = length . filter (==Boy)++boys :: Int -> Event Family+boys k f = countBoys f >= k++nBoys :: Int -> Int -> Int -> Probability+nBoys n k m = boys m ?? boys k ?=<< family n++numBoys :: Int -> Int -> Dist Int+numBoys n k = Dist.map countBoys (boys k ?=<< family n)+++-- * Special cases++-- | only boys in a family that has one boy+onlyBoys1 :: Int -> Probability+onlyBoys1 n = nBoys n 1 n
+ src/Numeric/Probability/Example/Predator.hs view
@@ -0,0 +1,87 @@+{- |+Lotka-Volterra predator-prey model++parameters++ * @g@ : victims' growth factor++ * @d@ : predators' death factor++ * @s@ : search rate++ * @e@ : energetic efficiency+-}++module Numeric.Probability.Example.Predator where++import Numeric.Probability.Visualize (+ Vis, Color(Green, Red),+ figP, figure, title,+ showParams, xLabel, yLabel, plotL, color, label,+ )+++-- try: n>=500+-- g = 1.05+-- d = 0.95+-- s = 0.01+-- e = 0.01+++g, d, s, e :: Float+g = 1.02+d = 0.98+s = 0.01+e = 0.01+++-- 'direct' function-over-time approach -- very inefficient due to recursion+--+-- v :: Int -> Float+-- v 0 = 20+-- v t = ((1 + r - a*p(t-1)) * v (t-1)) `max` 0+--+-- p :: Int -> Float+-- p 0 = 15+-- p t = ((1 - d + a*b*v(t-1)) * p (t-1)) `max` 0+--+--+-- fig1 = figP figure{title="Predator/Prey Simulation "+++-- showParams [r,d,a,b] ["r","d","a","b"],+-- xLabel="Time (generation)",+-- yLabel="Population"}+-- [(plotF (0,15,1) v){color=Green,label="Victim"},+-- (plotF (0,15,1) p){color=Red,label="Prey"}]++v0 :: Float+v0 = 1++p0 :: Float+p0 = 1++dv :: (Float,Float) -> Float+dv (v,p) = (g*v - s*v*p) `max` 0++dp :: (Float,Float) -> Float+dp (v,p) = (d*p + e*v*p) `max` 0++dvp :: (Float, Float) -> (Float, Float)+dvp vp' = (dv vp', dp vp')++vp :: [(Float, Float)]+vp = (v0,p0):map dvp vp++vs :: [Float]+vs = map fst vp++ps :: [Float]+ps = map snd vp+++fig1 :: Int -> Vis+fig1 n = figP figure{title="Predator/Prey Simulation "+++ showParams [g,d,s,e] ["g","d","s","e"],+ xLabel="Time (generation)",+ yLabel="Population"}+ [(plotL (take n vs)){color=Green,label="Victim"},+ (plotL (take n ps)){color=Red,label="Prey"}]
+ src/Numeric/Probability/Example/Queuing.hs view
@@ -0,0 +1,159 @@+{- |++Model:++ one server serving customers from one queue++-}++module Numeric.Probability.Example.Queuing where++import qualified Numeric.Probability.Distribution as Dist+import qualified Numeric.Probability.Random as Rnd++import Numeric.Probability.Percentage+ (Dist, RDist, Trans, )++import Data.List (nub,sort)+++{- no Random instance for Rational+type Probability = Rational+type Dist a = Dist.T Probability a+type RDist a = Rnd.Distribution Probability a+type Trans a = Transition Probability a+-}+++type Time = Int++-- | (servingTime, nextArrival)+type Profile = (Time, Time)++type Event a = (a,Profile)++-- | customers and their individual serving times+type Queue a = [(a,Time)]++-- | (customers waiting,validity period of that queue)+type State a = (Queue a,Time)++type System a = [([a],Time)]++type Events a = [Event a]+++event :: Time -> Events a -> Queue a -> [State a]+event = mEvent 1++--event _ [] [] = []+--event 0 ((c,(s,a)):es) q = event a es (q++[(c,s)])+--event a es [] = ([],a):event 0 es []+--event a [] (q@((c,s):q')) = (q,s):event a [] q'+--event a es (q@((c,s):q')) | a<s = (q,a):event 0 es ((c,s-a):q')+-- | True = (q,s):event (a-s) es q'++system :: Events a -> System a+--system es = map (\(q,t)->(map fst q,t)) $ event 0 es []+system = mSystem 1+++-- | multiple servers++mEvent :: Int -> Time -> Events a -> Queue a -> [State a]+mEvent _ _ [] [] = []+mEvent n 0 ((c,(s,a)):es) q = mEvent n a es (q++[(c,s)])+mEvent n a es [] = ([],a):mEvent n 0 es []+mEvent n _ [] q = (q,s):mEvent n 0 [] (mServe n s q)+ where s = mTimeStep n q+mEvent n a es q =+ if a < s+ then (q,a) : mEvent n 0 es (mServe n a q)+ else (q,s) : mEvent n (a-s) es (mServe n s q)+ where s = mTimeStep n q+++-- | decrease served customers remaining time by specified amount+mServe :: Int -> Int -> Queue a -> Queue a+mServe _ _ [] = []+mServe 0 _ x = x+mServe n c ((a,t):es) =+ if t > c+ then (a,t-c) : mServe (n-1) c es+ else mServe (n-1) c es++-- | time until next completion+mTimeStep :: Int -> Queue a -> Int+mTimeStep _ ((_,t):[]) = t+mTimeStep 1 ((_,t):_) = t+mTimeStep n ((_,t):es) = min t (mTimeStep (n-1) es)+mTimeStep _ _ = error "Queuing.mTimeStep: queue must be non-empty"++mSystem :: Int -> Events a -> System a+mSystem n es = map (\(q,t)->(map fst q,t)) $ mEvent n 0 es []+++-- * random++type RProfile = (Dist Time, Trans Time)++type REvent a = (a, RProfile)++type REvents a = [REvent a]++rSystem :: Int -> REvents a -> Rnd.T (System a)+rSystem n re = do+ e <- rBuildEvents re+ return (mSystem n e)++rBuildEvents :: REvents a -> Rnd.T (Events a)+rBuildEvents ((a,(dt,tt)):ex) = do+ rest <- rBuildEvents ex+ t <- Rnd.pick dt+ nt <- Rnd.pick $ tt t+ return ((a,(t,nt)):rest)+rBuildEvents [] = return []++rmSystem :: Ord a => Int -> Int -> REvents a -> RDist (System a)+rmSystem c n re = Rnd.dist $ replicate c (rSystem n re)++evalSystem :: (Ord a, Ord b) =>+ Int -> Int -> REvents a -> (System a -> b) -> RDist b+evalSystem c n re ef =+ do+ rds <- rmSystem c n re+ return (Dist.map ef rds)++unit :: b -> ((), b)+unit = (\p->((),p)) -- Dist.map (\p->((),p))+++-- * evaluation++maxQueue :: Ord a => System a -> Int+maxQueue s = maximum [length q | (q,_) <- s]++allWaiting :: Ord a => Int -> System a -> [a]+allWaiting n s = nub $ sort $ concat [ drop n q | (q,_) <- s]+++countWaiting :: Ord a => Int -> System a -> Int+countWaiting n = length . allWaiting n++waiting :: Int -> System a -> Time+waiting n s = sum [ t*length (drop n q) | (q,t) <- s]++inSystem :: System a -> Time+inSystem s = sum [ t*length q | (q,t) <- s]++total :: System a -> Time+total = sum . map snd++server :: Int -> System a -> Time+server n s = sum [ t*length (take n q) | (q,t) <- s]++idle :: Int -> System a -> Time+idle n s = sum [ t*(n - length q) | (q,t) <- s, length q <= n]++idleAvgP :: Int -> System a -> Float+idleAvgP n s = (fromIntegral $ idle n s) / (fromIntegral $ server n s)
+ src/Numeric/Probability/Example/TreeGrowth.hs view
@@ -0,0 +1,151 @@+module Numeric.Probability.Example.TreeGrowth where++import qualified Numeric.Probability.Distribution as Dist+import qualified Numeric.Probability.Transition as Trans+import qualified Numeric.Probability.Random as Rnd+import qualified Numeric.Probability.Trace as Trace+import Numeric.Probability.Simulation ((~..), (~*.), )+import Numeric.Probability.Percentage+ (Dist, Trans, RTrans, Expand, RExpand, Space, )++import Numeric.Probability.Visualize (+ Vis, Color(Green, Red, Blue), Plot,+ fig, figP, figure, title,+ xLabel, yLabel, plotD, color, label,+ )++import qualified Numeric.Probability.Monad as MonadExt+++type Height = Int++data Tree = Alive Height | Hit Height | Fallen+ deriving (Ord,Eq,Show)++grow :: Trans Tree+grow (Alive h) = Dist.normal (map Alive [h+1..h+5])+grow _ = error "TreeGrowth.grow: only alive trees can grow"++hit :: Trans Tree+hit (Alive h) = Dist.certainly (Hit h)+hit _ = error "TreeGrowth.hit: only alive trees can be hit"++fall :: Trans Tree+fall _ = Dist.certainly Fallen++evolve :: Trans Tree+evolve t =+ case t of+ (Alive _) -> Trans.unfold (Dist.enum [90,4,6] [grow,hit,fall]) t+-- (Alive _) -> Trans.unfold (Dist.relative [0.9,0.04,0.06] [grow,hit,fall]) t+ _ -> Dist.certainly t++{- |+tree growth simulation:+ start with seed and run for n generations+-}+seed :: Tree+seed = Alive 0+++-- * exact results++-- | @tree n@ : tree distribution after n generations+tree :: Int -> Tree -> Dist Tree+tree n = MonadExt.iterate n evolve++-- | @hist n@ : history of tree distributions for n generations+hist :: Int -> Expand Tree+hist n = Trace.walk n (evolve =<<) . return+++-- * simulation results++{- |+Since '(*.)' is overloaded for Trans and RChange,+we can run the simulation ~. directly to @n *. live@.+-}++--simTree k n = k ~. tree n+simTree :: Int -> Int -> RTrans Tree+simTree k n = (k,n) ~*. evolve++simHist :: Int -> Int -> RExpand Tree+simHist k n = (k,n) ~.. evolve++t2 :: Dist Tree+t2 = tree 2 seed++h2 :: Space Tree+h2 = hist 2 seed++sh2, st2 :: IO ()+st2 = Rnd.print $ simTree 2000 2 seed+sh2 = Rnd.print $ simHist 2000 2 seed+++-- Alternatives:+--+-- simTree k n = k ~. n *. random evolve+-- simTree k n = (k,n) ~*. evolve+++-- take a trace+++height :: Tree -> Int+height Fallen = 0+height (Hit h) = h+height (Alive h) = h+{--+myPlot = plotD ((5 *. evolve) (Alive 0) >>= height)++myPlot2 = figP figure{title="Tree Growth",xLabel="Height (m)",+ yLabel="Probability"}+ (autoColor [+ plotD ((5 *. evolve) (Alive 0) >>= height)+ ])++--}++p1, p2, p3, p4, p5, p6 :: Vis++p1 = fig [plotD $ Dist.normal ([1..20]::[Int])]++p2 = fig [plotD $ Dist.map height (tree 5 seed)]++p3 = figP figure{title="Tree Growth",+ xLabel="Height (ft)",+ yLabel="Probability"}+ [plotD $ Dist.map height (tree 5 seed)]+++p4 = figP figure{title="Tree Growth",+ xLabel="Height (ft)",+ yLabel="Probability"}+ [heightAtTime 5, heightAtTime 10,heightAtTime 15]++heightAtTime :: Int -> Plot+heightAtTime y = plotD $ Dist.map height (tree y seed)++p5 = figP figure{title="Tree Growth",+ xLabel="Height (ft)",+ yLabel="Probability"}+ (map heightAtTime [3,5,7])++heightCurve :: (Int,Color) -> Plot+heightCurve (n,c) = (heightAtTime n){color=c,label=show n++" Years"}++p6 = figP figure{title="Tree Growth",+ xLabel="Height (ft)",+ yLabel="Probability"}+ (map heightCurve+ [(3,Blue),(5,Green),(7,Red)])+++done :: Tree -> Bool+done (Alive x) = x >= 5+done _ = True++ev5 :: Tree -> Dist Tree+ev5 = MonadExt.while (not . done) evolve
+ src/Numeric/Probability/Expectation.hs view
@@ -0,0 +1,57 @@+module Numeric.Probability.Expectation where++import qualified Numeric.Probability.Distribution as Dist+import qualified Numeric.Probability.Percentage as Probability++-- TO DO: generalize Float to arbitrary Num type+--+class ToFloat a where+ toFloat :: a -> Float++instance ToFloat Float where toFloat = id+instance ToFloat Int where toFloat = fromIntegral+instance ToFloat Integer where toFloat = fromIntegral+instance ToFloat Probability.T+ where toFloat (Probability.Cons x) = x++class FromFloat a where+ fromFloat :: Float -> a++instance FromFloat Float where fromFloat = id+instance FromFloat Int where fromFloat = round+instance FromFloat Integer where fromFloat = round++-- expected :: ToFloat a => Prob.Dist a -> Float+-- expected = sum . map (\(x,p)->toFloat x*p) . Dist.decons++class Expected a where+ expected :: a -> Float++-- instance ToFloat a => Expected a where+-- expected = toFloat+instance Expected Float where expected = id+instance Expected Int where expected = toFloat+instance Expected Integer where expected = toFloat++instance Expected a => Expected [a] where+ expected = Dist.expected . Dist.uniform . map expected+-- expected xs = sum (map expected xs) / toFloat (length xs)++floatDist :: (ToFloat prob, Expected a) =>+ Dist.T prob a -> Dist.T Float Float+floatDist =+ Dist.Cons .+ map (\(x,p) -> (expected x, toFloat p)) .+ Dist.decons++instance (ToFloat prob, Expected a) => Expected (Dist.T prob a) where+ expected = Dist.expected . floatDist+-- expected = Dist.expected . fmap expected+++-- | statistical analyses+variance :: Expected a => Probability.Dist a -> Float+variance = Dist.variance . floatDist++stdDev :: Expected a => Probability.Dist a -> Float+stdDev = sqrt . variance
+ src/Numeric/Probability/Monad.hs view
@@ -0,0 +1,53 @@+-- | Monad helper functions+module Numeric.Probability.Monad where++import Control.Monad (liftM, )++-- | 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+++iterate :: Monad m => Int -> (a -> m a) -> (a -> m a)+iterate n f = compose $ replicate n f+++-- | like 'iterate' but returns all intermediate values+walk :: (Monad m) => Int -> (a -> m a) -> (a -> m [a])+walk n f =+ let recurse 0 _ = return []+ recurse m x = liftM (x:) (recurse (pred m) =<< f x)+ in recurse n+++{- |+While loop with a posteriori check.+Loops at least once.+-}+doWhile :: Monad m => (a -> Bool) -> (a -> m a) -> (a -> m a)+doWhile p t =+ let recurse x = t x >>= \l -> if p l then recurse l else return l+ in recurse++{- |+While loop with a priori check.+Can loop zero times.+-}+while :: Monad m => (a -> Bool) -> (a -> m a) -> (a -> m a)+while p t =+ let recurse x = if p x then t x >>= recurse else return x+ in recurse+++whileTrace :: Monad m => (a -> m Bool) -> m a -> m [a]+whileTrace p t =+ do x <- t+ b <- p x+ liftM (x:) $+ if b+ then whileTrace p t+ else return []
+ src/Numeric/Probability/Object.hs view
@@ -0,0 +1,63 @@+{-# OPTIONS_GHC -fglasgow-exts #-}+{- |+Portability: Multi-parameter type class with functional dependency++Abstract interface to probabilistic objects+like random generators and probability distributions.+It allows to use the same code+both for computing complete distributions+and for generating random values according to the distribution.+The latter one is of course more efficient+and may be used for approximation of the distribution by simulation.+-}+module Numeric.Probability.Object where++import qualified Numeric.Probability.Distribution as Dist+import qualified Numeric.Probability.Random as Rnd+import qualified Numeric.Probability.Shape as Shape++import qualified Data.List as List+++class Monad obj => C prob obj | obj -> prob where+ fromFrequencies :: [(a,prob)] -> obj a+++instance C Double Rnd.T where+ fromFrequencies = Rnd.pick . Dist.fromFreqs++instance Fractional prob => C prob (Dist.T prob) where+ fromFrequencies = Dist.fromFreqs++++type Spread obj a = [a] -> obj a++shape :: (C prob obj, Fractional prob) =>+ (prob -> prob) -> Spread obj a+shape _ [] = error "Probability.shape: empty list"+shape f xs =+ let incr = 1 / fromIntegral (length xs - 1)+ ps = List.map f (iterate (+incr) 0)+ in fromFrequencies (zip xs ps)++linear :: (C prob obj, Fractional prob) => Spread obj a+linear = shape Shape.linear++uniform :: (C prob obj, Fractional prob) => Spread obj a+uniform = shape Shape.uniform++negExp :: (C prob obj, Floating prob) => Spread obj a+negExp = shape Shape.negExp++normal :: (C prob obj, Floating prob) => Spread obj a+normal = shape Shape.normal++enum :: (C prob obj, Floating prob) => [Int] -> Spread obj a+enum = relative . List.map fromIntegral++{- |+Give a list of frequencies, they do not need to sum up to 1.+-}+relative :: (C prob obj, Floating prob) => [prob] -> Spread obj a+relative ns = fromFrequencies . flip zip ns
+ src/Numeric/Probability/Percentage.hs view
@@ -0,0 +1,127 @@+{- |+Number type based on Float with formatting in percents.+-}+module Numeric.Probability.Percentage where++import qualified Numeric.Probability.Distribution as Dist+import qualified Numeric.Probability.Random as Rnd++import Numeric.Probability.Show (showR)+import Numeric.Probability.Trace (Trace)++import qualified System.Random as Random+++-- ** Probabilities+newtype T = Cons Float+ deriving (Eq, Ord)+++percent :: Float -> T+percent x = Cons (x/100)+++showPfix :: (RealFrac prob) => Int -> prob -> String+showPfix precision f =+ if precision==0+ then showR 3 (round (f*100) :: Integer)++"%"+ else showR (4+precision) (roundRel precision (f*100))++"%"++roundRel :: (RealFrac a) => Int -> a -> a+roundRel p x =+ let d = 10^p+ in fromIntegral (round (x*d) :: Integer)/d++-- -- mixed precision+-- --+-- showP :: ProbRep -> String+-- showP f | f>=0.1 = showR 3 (round (f*100))++"%"+-- | otherwise = show (f*100)++"%"++-- fixed precision+--+-- showP :: ProbRep -> String+-- showP = showPfix 1+++instance Show T where+ show (Cons p) = showPfix 1 p++++infix 0 //++{- |+Print distribution as table with configurable precision.+-}+(//) :: (Ord a, Show a) => Dist a -> Int -> IO ()+(//) x prec = putStr (Dist.pretty (\(Cons p) -> showPfix prec p) x)+++++liftP :: (Float -> Float) -> T -> T+liftP f (Cons x) = Cons (f x)++liftP2 :: (Float -> Float -> Float) -> T -> T -> T+liftP2 f (Cons x) (Cons y) = Cons (f x y)++instance Num T where+ fromInteger = Cons . fromInteger+ (+) = liftP2 (+)+ (-) = liftP2 (-)+ (*) = liftP2 (*)+ abs = liftP abs+ signum = liftP signum+ negate = liftP negate++instance Fractional T where+ fromRational = Cons . fromRational+ recip = liftP recip+ (/) = liftP2 (/)++instance Floating T where+ pi = Cons pi+ exp = liftP exp+ sqrt = liftP sqrt+ log = liftP log+ (**) = liftP2 (**)+ logBase = liftP2 logBase+ sin = liftP sin+ tan = liftP tan+ cos = liftP cos+ asin = liftP asin+ atan = liftP atan+ acos = liftP acos+ sinh = liftP sinh+ tanh = liftP tanh+ cosh = liftP cosh+ asinh = liftP asinh+ atanh = liftP atanh+ acosh = liftP acosh++instance Random.Random T where+ randomR (Cons l, Cons r) =+ (\(x,g) -> (Cons x, g)) . Random.randomR (l,r)+ random =+ (\(x,g) -> (Cons x, g)) . Random.random+ randomRIO (Cons l, Cons r) = fmap Cons $ Random.randomRIO (l,r)+ randomIO = fmap Cons $ Random.randomIO+++type Dist a = Dist.T T a+++type Spread a = [a] -> Dist a++type RDist a = Rnd.T (Dist a)++type Trans a = a -> Dist a++type Space a = Trace (Dist a)+type Expand a = a -> Space a++type RTrans a = a -> RDist a++type RSpace a = Rnd.T (Space a)+type RExpand a = a -> RSpace a
+ src/Numeric/Probability/PrintList.hs view
@@ -0,0 +1,54 @@+-- | Utilities for printing lists+module Numeric.Probability.PrintList where+++import Data.List (intersperse)+++----------------------------------------------------------------------+-- PRINT UTILITIES+----------------------------------------------------------------------++newtype Lines a = Lines [a]++instance Show a => Show (Lines a) where+ show (Lines xs) = printList ("","\n","") show xs++asLines :: [a] -> Lines a+asLines = Lines+++showNQ :: Show a => a -> String+showNQ = filter ('"'/=) . show++indent :: Int -> Int -> [Char]+indent i l = take (i*l) (repeat ' ')++printList :: ([a],[a],[a]) -> (b -> [a]) -> [b] -> [a]+printList (sep0,sep1,sep2) f xs =+ sep0++concat (intersperse sep1 (map f xs))++sep2+++asTuple, asSeq, asList, asSet, asLisp,+ asString, asPlain, asPlain' :: (a -> [Char]) -> [a] -> [Char]++asTuple = printList ("(",",",")")+asSeq = printList ("",",","")+asList = printList ("[",",","]")+asSet = printList ("{",",","}")+asLisp = printList ("("," ",")")+asPlain f xs = if null xs then "" else printList (" "," ","") f xs+asPlain' f xs = if null xs then "" else printList (""," ","") f xs+asString = printList ("","","")+-- asLines = printList ["","\n",""]++asCases :: Int -> (a -> [Char]) -> [a] -> [Char]+asCases l =+ let ind = indent 4 l+ in printList ("\n"++ind++" ","\n"++ind++" | ","")++asDefs :: [Char] -> (a -> [Char]) -> [a] -> [Char]+asDefs n = printList ("\n"++n,"\n"++n,"\n")++asParagraphs :: (a -> [Char]) -> [a] -> [Char]+asParagraphs = printList ("\n","\n\n","\n")
+ src/Numeric/Probability/Random.hs view
@@ -0,0 +1,102 @@+-- | Randomized values+module Numeric.Probability.Random where++import qualified Numeric.Probability.Distribution as Dist+import qualified Numeric.Probability.Transition as Trans++import qualified System.Random as Random+import System.Random (Random, )++import Control.Monad.State (State(State), evalState, )++import qualified System.IO as IO+import Prelude hiding (print)+++-- * random generator++-- | Random values+newtype T a = Cons {decons :: State Random.StdGen a}++instance Monad T where+ return x = Cons (return x)+ Cons x >>= y =+ Cons (decons . y =<< x)++instance Functor T where+ fmap f = Cons . fmap f . decons+++randomR :: Random.Random a => (a, a) -> T a+randomR rng =+ Cons (State (Random.randomR rng))++{- |+Run random action in 'IO' monad.+-}+run :: T a -> IO a+run r = fmap (evalState (decons r)) Random.getStdGen++{- |+Run random action without 'IO' using a seed.+-}+runSeed :: Random.StdGen -> T a -> a+runSeed g r = evalState (decons r) g+++print :: Show a => T a -> IO ()+print r = IO.print =<< run r++-- instance Show (IO a) where+-- show _ = ""++pick :: (Num prob, Ord prob, Random prob) =>+ Dist.T prob a -> T a+pick d = return . Dist.selectP d =<< randomR (0,1)+++-- * random distribution++-- | Randomized distribution+type Distribution prob a = T (Dist.T prob a)++above :: (Num prob, Ord prob, Ord a) =>+ prob -> Distribution prob a -> Distribution prob (Dist.Select a)+above p rd = fmap (Dist.above p) rd++{- |+'dist' converts a list of randomly generated values into+a distribution by taking equal weights for all values+-}+dist :: (Fractional prob, Ord a) => [T a] -> Distribution prob a+dist = fmap (Dist.norm . Dist.uniform) . sequence+++-- * Randomized changes++-- | random change+type Change a = a -> T a++change :: (Num prob, Ord prob, Random prob) =>+ Trans.T prob a -> Change a+change t = pick . t+++-- * Randomized transitions++-- | random transition+type Transition prob a = a -> Distribution prob a++type ApproxDist a = T [a]+++++{-+for quickCheck++LAWS++ const . pick = random . const++-}
+ src/Numeric/Probability/Shape.hs view
@@ -0,0 +1,32 @@+{- |+Collection of some shapes of distribution.+-}+module Numeric.Probability.Shape where++{- |+A shape is a mapping from the interval @[0,1]@ to non-negative numbers.+They need not to be normalized (sum up to 1)+because this is done by subsequent steps.+(It would also be impossible to normalize the function in a way+that each discretization is normalized as well.)+-}+type T prob = prob -> prob+++linear :: Fractional prob => T prob+linear = id++uniform :: Fractional prob => T prob+uniform = const 1++negExp :: Floating prob => T prob+negExp x = exp (-x)++normal :: Floating prob => T prob+normal = normalCurve 0.5 0.5++normalCurve :: Floating prob =>+ prob -> prob -> prob -> prob+normalCurve mean dev x =+ let u = (x - mean) / dev+ in exp (-1/2 * u^(2::Int)) / sqrt (2 * pi)
+ src/Numeric/Probability/Show.hs view
@@ -0,0 +1,16 @@+module Numeric.Probability.Show where++showL :: Show a => Int -> a -> String+showL n x = s++rep (n-length s) ' '+ where s=show x++showR :: Show a => Int -> a -> String+showR n x = rep (n-length s) ' '++s+ where s=show x++--showP :: Float -> String+--showP f = showR 3 (round (f*100))++"%"++rep :: Int -> a -> [a]+rep n x = take n (repeat x)+
+ src/Numeric/Probability/Simulation.hs view
@@ -0,0 +1,96 @@+-- | Simulation+module Numeric.Probability.Simulation where++import qualified Numeric.Probability.Distribution as Dist+import qualified Numeric.Probability.Random as Rnd+import qualified Numeric.Probability.Trace as Trace++import System.Random (Random, )++import qualified Numeric.Probability.Monad as MonadExt+++{-++Naming convention:++ * @*@ takes @n :: Int@ and a generator and iterates the generator n times++ * @.@ produces a single result++ * @..@ produces a trace++ * @~@ takes @k :: Int@ [and @n :: Int@] and a generator and simulates+ the [n-fold repetition of the] generator k times+++There are the following functions:++ * @n *. t@ iterates t and produces a distribution++ * @n *.. t@ iterates t and produces a trace++ * @k ~. t@ simulates t and produces a distribution++ * @(k,n) ~*. t@ simulates the n-fold repetition of t and produces a distribution++ * @(k,n) ~.. t@ simulates the n-fold repetition of t and produces a trace+++Iteration captures three iteration strategies:+iter builds an n-fold composition of a (randomized) transition+while and until implement conditional repetitions++The class Iterate allows the overloading of iteration for different+kinds of generators, namely transitions and Rnd.change changes:++ * @Trans a = a -> Dist a ==> c = Dist@++ * @RChange a = a -> Rnd.T a ==> c = Rnd.T = IO@++-}+++{- |+Simulation means to repeat a Rnd.change change many times and+to accumulate all results into a distribution. Therefore,+simulation can be regarded as an approximation of distributions+through randomization.++The Sim class allows the overloading of simulation for different+kinds of generators, namely transitions and Rnd.change changes:++ * @Trans a = a -> Dist a ==> c = Dist@++ * @RChange a = a -> Rnd.T a ==> c = Rnd.T = IO@+-}+class C c where+ -- | returns the final randomized transition+ (~.) :: (Fractional prob, Ord prob, Random prob, Ord a) =>+ Int -> (a -> c a) -> Rnd.Transition prob a+ -- | returns the whole trace for a k-fold simulation+ (~..) :: (Fractional prob, Ord prob, Random prob, Ord a) =>+ (Int,Int) -> (a -> c a) -> Trace.RExpand prob a+ -- | returns the whole trace for a single simulation+ (~*.) :: (Fractional prob, Ord prob, Random prob, Ord a) =>+ (Int,Int) -> (a -> c a) -> Rnd.Transition prob a++infix 6 ~. , ~..++infix 8 ~*.+++-- simulation for transitions+--+instance (Num prob, Ord prob, Random prob) => C (Dist.T prob) where+ (~.) x = (~.) x . Rnd.change+ (~..) x = (~..) x . Rnd.change+ (~*.) x = (~*.) x . Rnd.change+++-- simulation for Rnd.change changes+--+instance C Rnd.T where+ (~.) n t = Rnd.dist . replicate n . t+ (~..) (k,n) t = Trace.merge . replicate k . MonadExt.walk n t+ (~*.) (k,n) t = k ~. MonadExt.iterate n t
+ src/Numeric/Probability/Trace.hs view
@@ -0,0 +1,50 @@+-- | Tracing+module Numeric.Probability.Trace where++import qualified Numeric.Probability.Distribution as Dist+import qualified Numeric.Probability.Transition as Trans+import qualified Numeric.Probability.Random as Rnd++import Data.List (transpose)+++-- * traces of distributions++type Trace a = [a]+type Walk a = a -> Trace a++type Space prob a = Trace (Dist.T prob a)+type Expand prob a = a -> Space prob a++++-- for ListUtils+-- | walk is a bounded version of the predefined function iterate+walk :: Int -> Trans.Change a -> Walk a+walk n f = take n . iterate f++++-- * traces of random experiments++type RTrace a = Rnd.T (Trace a)+type RWalk a = a -> RTrace a++type RSpace prob a = Rnd.T (Space prob a)+type RExpand prob a = a -> RSpace prob a+++{- |+'merge' converts a list of 'RTrace's+into a list of randomized distributions, i.e., an 'RSpace',+by creating a randomized distribution for each list position across all traces+-}+merge :: (Fractional prob, Ord a) =>+ [RTrace a] -> RSpace prob a+merge =+ fmap (zipListWith (Dist.norm . Dist.uniform)) . sequence+++-- for ListUtils+zipListWith :: ([a] -> b) -> [[a]] -> [b]+zipListWith f = map f . transpose
+ src/Numeric/Probability/Transition.hs view
@@ -0,0 +1,108 @@+-- | Deterministic and probabilistic generators+module Numeric.Probability.Transition where++import qualified Numeric.Probability.Distribution as Dist++import qualified Data.List as List+import Prelude hiding (map, maybe, id, )+++-- * Transitions+++-- | deterministic generator+type Change a = a -> a++-- | probabilistic generator+type T prob a = a -> Dist.T prob a+++id :: (Num prob) => T prob a+id = Dist.certainly+++{- |+'map' maps a change function to the result of a transformation+('map' is somehow a lifted form of 'Dist.map')+The restricted type of @f@ results from the fact that the+argument to @t@ cannot be changed to @b@ in the result 'T' type.+-}+map :: (Num prob, Ord a) =>+ Change a -> T prob a -> T prob a+map f t = Dist.map f . t+++{- |+unfold a distribution of transitions into one transition++NOTE: The argument transitions must be independent+-}+unfold :: (Num prob, Ord a) =>+ Dist.T prob (T prob a) -> T prob a+unfold d x = Dist.unfold (fmap ($x) d)++{- |+Composition of transitions similar to 'Numeric.Probability.Monad.compose'+but with intermediate duplicate elimination.+-}+compose :: (Num prob, Ord a) =>+ [T prob a] -> T prob a+compose = foldl (\acc x v -> Dist.norm (acc v >>= x)) return++++-- * Spreading changes into transitions++-- | functions to convert a list of changes into a transition+type SpreadC prob a = [Change a] -> T prob a++apply :: (Num prob) =>+ Change a -> T prob a+apply f = id . f+++maybe :: (Num prob) => prob -> Change a -> T prob a+maybe p f x = Dist.choose p (f x) x++lift :: Dist.Spread prob a -> SpreadC prob a+lift s cs x = s $ List.map ($ x) cs++uniform :: (Fractional prob) => SpreadC prob a+uniform = lift Dist.uniform++linear :: (Fractional prob) => SpreadC prob a+linear = lift Dist.linear++normal :: (Floating prob) => SpreadC prob a+normal = lift Dist.normal++enum :: (RealFloat prob) => [Int] -> SpreadC prob a+enum xs = lift (Dist.enum xs)++relative :: (RealFloat prob) => [prob] -> SpreadC prob a+relative xs = lift (Dist.relative xs)+++-- * Spreading transitions into transitions++-- | functions to convert a list of transitions into a transition+type SpreadT prob a = [T prob a] -> T prob a++liftT :: (Num prob, Ord a) =>+ Dist.Spread prob (T prob a) -> SpreadT prob a+liftT s = unfold . s++uniformT :: (Fractional prob, Ord a) => SpreadT prob a+uniformT = liftT Dist.uniform++linearT :: (Fractional prob, Ord a) => SpreadT prob a+linearT = liftT Dist.linear++normalT :: (Floating prob, Ord a) => SpreadT prob a+normalT = liftT Dist.normal++enumT :: (RealFloat prob, Ord a) => [Int] -> SpreadT prob a+enumT xs = liftT (Dist.enum xs)++relativeT :: (RealFloat prob, Ord a) => [prob] -> SpreadT prob a+relativeT xs = liftT (Dist.relative xs)
+ src/Numeric/Probability/Visualize.hs view
@@ -0,0 +1,244 @@+module Numeric.Probability.Visualize where++import qualified Numeric.Probability.Random as Rnd+import Numeric.Probability.Expectation+ (ToFloat, FromFloat, toFloat, fromFloat, )+import Numeric.Probability.Percentage+ (Dist, RDist, )+import Numeric.Probability.PrintList (asTuple, )++import qualified Numeric.Probability.Distribution as Dist++import Data.List (nub, sort, )+++{- TO DO:++* Change function representation in Plot to+ xs :: [Float]+ ys :: [Float]+ and add functions to create this representation from+ functions, distributions, and lists+ (i.e. plotF, plotD, plotL)++-}+++-- | global settings for one figure+--+data FigureEnv = FE { fileName :: String,+ title :: String,+ xLabel :: String,+ yLabel :: String }+ deriving Show++-- | default settings for figure environment+--+figure :: FigureEnv+figure = FE { fileName = "FuSE.R",+ title = "Output",+ xLabel = "x",+ yLabel = "f(x)" }+++-- * types to represent settings for individual plots+--+data Color = Black | Blue | Green | Red | Brown | Gray+ | Purple | DarkGray | Cyan | LightGreen | Magenta+ | Orange | Yellow | White | Custom Int Int Int+ deriving Eq++instance Show Color where+ show Black = "\"black\""+ show Blue = "\"blue\""+ show Green = "\"green\""+ show Red = "\"red\""+ show Brown = "\"brown\""+ show Gray = "\"gray\""+ show Purple = "\"purple\""+ show DarkGray = "\"darkgray\""+ show Cyan = "\"cyan\""+ show LightGreen = "\"lightgreen\""+ show Magenta = "\"magenta\""+ show Orange = "\"orange\""+ show Yellow = "\"yellow\""+ show White = "\"white\""+ show (Custom r g b) = "rgb("++(show r)++", "++(show g)++", "++(show b)++")"++data LineStyle = Solid | Dashed | Dotted | DotDash | LongDash | TwoDash+ deriving Eq++instance Show LineStyle where+ show Solid = "1"+ show Dashed = "2"+ show Dotted = "3"+ show DotDash = "4"+ show LongDash = "5"+ show TwoDash = "6"++type PlotFun = Float -> Float+++-- | settings for individual plots+--+data Plot = Plot { ys :: [Float],+ xs :: [Float],+ color :: Color,+ lineStyle :: LineStyle,+ lineWidth :: Int,+ label :: String }++{-+instance Show Plot where+ show _ = "Individual plots cannot be printed.\nPlease use plots \+ \ as arguments to the fig function."+-}+++-- | default plotting environment+--+plot :: Plot+plot = Plot { ys = [0],+ xs = [0],+ color = Black,+ lineStyle = Solid,+ lineWidth = 1,+ label = "" }++colors :: [Color]+colors = [Blue,Green,Red,Purple,Black,Orange,Brown,Yellow]++setColor :: Plot -> Color -> Plot+setColor p c = p{color=c}++autoColor :: [Plot] -> [Plot]+autoColor ps | length ps <= n = zipWith setColor ps colors+ | otherwise = error ("autoColor works for no more than "+++ show n++" plots.")+ where n=length colors++-- | create a plot from a distribution+--+plotD :: ToFloat a => Dist a -> Plot+--plotD d = plot{ys = map (\x->(dp $ prob' x d')) (extract d'),+-- xs = extract d'}+plotD d =+ let (tfl, pdl) =+ unzip $ Dist.sortElem $+ Dist.norm' (map (\(x,p) -> (toFloat x, toFloat p)) (Dist.decons d))+ in plot{xs = tfl, ys = pdl}+++plotRD :: ToFloat a => RDist a -> IO Plot+plotRD a = Rnd.run (fmap plotD a)++-- | create a plot from a function+--+plotF :: (FromFloat a,ToFloat b) => (Float,Float,Float) -> (a -> b) -> Plot+plotF xd g = plot{ys = map (\x->toFloat (g (fromFloat x))) (xvals xd),xs = xvals xd}+ where xvals (a,b,d) =+ if a > b then [] else a:xvals (a+d,b,d)++-- | create a plot from a list+--+plotL :: ToFloat a => [a] -> Plot+plotL vs = plot{ys = map toFloat vs, xs = map toFloat [1..length vs]}+++plotRL :: ToFloat a => Rnd.T [a] -> IO Plot+plotRL a = Rnd.run (fmap plotL a)+++--yls :: ToFloat a => [a] -> [Plot] -> [[Float]]+--yls xs (p:ps) = [f p (toFloat v) | v <- xs ]:yls xs ps+--yls _ [] = []++yls :: [Float] -> Plot -> Plot+yls xl p = p{xs=x', ys=y'}+ where t = zip (xs p) (ys p)+ t' = metaTuple xl t+ (x', y') = unzip t'++metaTuple :: [Float] -> [(Float,Float)] -> [(Float,Float)]+metaTuple (x:xl) ((p,v):px) | p == x = (p,v):(metaTuple xl px)+metaTuple (x:xl) p'@( (p,_):_ ) | p > x = (x,0):(metaTuple xl p')+metaTuple x [] = map (\v->(v,0)) x+metaTuple x y = error $ (show x)++(show y)++-- | we want to increase the bounds absolutely, account for negative numbers+--+incr, decr :: (Ord a, Fractional a) => a -> a+incr x =+ if x > 0+ then x * 1.05+ else x * 0.95++decr x =+ if x > 0+ then x * 0.95+ else x * 1.05++-- | Visualization output+--+type Vis = IO ()+++-- * creating figures+--+fig :: [Plot] -> Vis+fig = figP figure++figP :: FigureEnv -> [Plot] -> Vis+figP fe ps = do let xl = sort $ nub $ concatMap xs ps+ let minx = minimum xl+-- let maxx = maximum xl+ let n = length xl+ let ys' = map ys (map (yls xl) ps) -- yls xl ps+ let miny = minimum (map minimum ys')+ let maxy = maximum (map maximum ys')+ let out0' = out0 (fileName fe)+ let out1' = out1 (fileName fe)+ out0' ("x <- "++(vec xl))+ out1' ("y <- "++(vec $ (decr miny):(replicate (n-1) (incr maxy))))+ out1' ("plot(x,y,type=\"n\",main=\""+++ 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\")")+++{-+define:+ * autoLabel+ * showParams+-}++showParams :: Show a => [a] -> [String] -> String+showParams xs0 ss =+ asTuple id (zipWith (\x s-> show x++":"++s) xs0 ss)++legend :: Float -> Float -> [Plot] -> String+legend x y ps = "legend("++(show x)++", "++(show y)++","+++ "lty="++vec (map lineStyle ps)++","+++ "col="++vec (map color ps)++","+++ "lwd="++vec (map lineWidth ps)++","+++ "legend="++vec (map label ps)++")"++drawy :: ToFloat a => Int -> Plot -> [a] -> String+drawy yn p fl = "y"++(show yn)++" <- "++(vec (map toFloat fl))++"\n"+++ "lines(x,y"++(show yn)++",col="++(show $ color p)++","+++ "lty="++(show $ lineStyle p)++",lwd="++(show $ lineWidth p)++")"+++vec :: Show a => [a] -> String+vec xs0 = "c"++asTuple show xs0++out0 :: String -> String -> IO ()+out0 f s = writeFile (f) (s++"\n")++out1 :: String -> String -> IO ()+out1 f s = appendFile (f) (s++"\n")