packages feed

probability (empty) → 0.1

raw patch · 21 files changed

+2011/−0 lines, 21 filesdep +basedep +haskell98setup-changed

Dependencies added: base, haskell98

Files

+ Alarm.hs view
@@ -0,0 +1,58 @@+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 view
@@ -0,0 +1,56 @@+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 view
@@ -0,0 +1,95 @@+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 view
@@ -0,0 +1,47 @@+{- |+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+
@@ -0,0 +1,28 @@+Copyright (c) 2005, Martin Erwig and Steve Kollmansberger+All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++1. Redistributions of source code must retain the above copyright notice,+   this list of conditions and the following disclaimer.++2. Redistributions in binary form must reproduce the above copyright+   notice, this list of conditions and the following disclaimer in the+   documentation and/or other materials provided with the distribution.++3. Neither the name of the author nor the names of its contributors may be+   used to endorse or promote products derived from this software without+   specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE+ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE+LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR+CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF+SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS+INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN+CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)+ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE+POSSIBILITY OF SUCH DAMAGE.
+ Collection.hs view
@@ -0,0 +1,112 @@+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 view
@@ -0,0 +1,39 @@+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 view
@@ -0,0 +1,16 @@+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 view
@@ -0,0 +1,79 @@+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 view
@@ -0,0 +1,46 @@+{- |+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 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 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 view
@@ -0,0 +1,54 @@+-- | 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 view
@@ -0,0 +1,655 @@+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 view
@@ -0,0 +1,144 @@+{- |++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)
+ README view
@@ -0,0 +1,37 @@+Probabilistic Functional Programming in Haskell++Contact:+Martin Erwig, Oregon State University, erwig@eecs.oregonstate.edu+++These files have been tested with GHC 6.4++Core Library files:++Show.hs		Pretty Printing+ListUtils.hs	+PrintList.hs	+Probability.hs	Core probabilistic module+Visualize.hs	Visualization system for use with R++Examples:++Barber.hs		An example of the queueing system+BayesianNetwork.hs	Implementing Bayesian networks+Boys.hs			A statistical examples+NBoys.hs		A generalized version of the previous+Collection.hs		Collections and two examples:+			Marbles and cards+Dice.hs			Rolling dice+MontyHall.hs		The "Monty Hall" Game (statistical)+Predator.hs		Non-probabilistic, demonstrates visualization+TreeGrowth.hs		A simple tree growth example++++Visualize output is placed in the file FuSE.R which can be loaded into the +R statistical program to see visualizations.++Randomized values can be displayed to the console using the printR +function, which shows the value from a IO monad function.+
+ Setup.lhs view
@@ -0,0 +1,3 @@+#! /usr/bin/env runhaskell+> import Distribution.Simple+> main = defaultMain
+ Show.hs view
@@ -0,0 +1,16 @@+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
@@ -0,0 +1,14 @@+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
+ TreeGrowth.hs view
@@ -0,0 +1,143 @@+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 view
@@ -0,0 +1,241 @@+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
@@ -0,0 +1,41 @@+Name:               probability+Version:            0.1+License:            BSD3+Author:             Martin Erwig <erwig@eecs.oregonstate.edu>+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+Description:+   The Library allows exact computation with discrete random variables+   in terms of their distributions by using a monad.+   The monad is similar to the List monad for non-deterministic computations,+   but extends the List monad by a measure of probability.+   Small interface to R plotting.+Tested-With:        GHC==6.4+Build-Type:         Simple+License-File:       COPYRIGHT+Hs-Source-Dirs:     .+Exposed-Modules:+    Alarm+    Barber+    Bayesian+    Boys+    Collection+    Dice+    MontyHall+    NBoys+    Predator+    Probability+    Queuing+    TreeGrowth+    Visualize+Other-Modules:+    ListUtils+    PrintList+    Show+Extra-Source-Files:+    README+    ToDo+GHC-Options:        -Wall -O2