diff --git a/Alarm.hs b/Alarm.hs new file mode 100644 --- /dev/null +++ b/Alarm.hs @@ -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 diff --git a/Barber.hs b/Barber.hs new file mode 100644 --- /dev/null +++ b/Barber.hs @@ -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) ) diff --git a/Bayesian.hs b/Bayesian.hs new file mode 100644 --- /dev/null +++ b/Bayesian.hs @@ -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) diff --git a/Boys.hs b/Boys.hs new file mode 100644 --- /dev/null +++ b/Boys.hs @@ -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 + diff --git a/COPYRIGHT b/COPYRIGHT new file mode 100644 --- /dev/null +++ b/COPYRIGHT @@ -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. diff --git a/Collection.hs b/Collection.hs new file mode 100644 --- /dev/null +++ b/Collection.hs @@ -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) diff --git a/Dice.hs b/Dice.hs new file mode 100644 --- /dev/null +++ b/Dice.hs @@ -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 diff --git a/ListUtils.hs b/ListUtils.hs new file mode 100644 --- /dev/null +++ b/ListUtils.hs @@ -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' diff --git a/MontyHall.hs b/MontyHall.hs new file mode 100644 --- /dev/null +++ b/MontyHall.hs @@ -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 diff --git a/NBoys.hs b/NBoys.hs new file mode 100644 --- /dev/null +++ b/NBoys.hs @@ -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 + diff --git a/Predator.hs b/Predator.hs new file mode 100644 --- /dev/null +++ b/Predator.hs @@ -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"}] diff --git a/PrintList.hs b/PrintList.hs new file mode 100644 --- /dev/null +++ b/PrintList.hs @@ -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") diff --git a/Probability.hs b/Probability.hs new file mode 100644 --- /dev/null +++ b/Probability.hs @@ -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) [(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 + +-} diff --git a/Queuing.hs b/Queuing.hs new file mode 100644 --- /dev/null +++ b/Queuing.hs @@ -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 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) diff --git a/README b/README new file mode 100644 --- /dev/null +++ b/README @@ -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. + diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,3 @@ +#! /usr/bin/env runhaskell +> import Distribution.Simple +> main = defaultMain diff --git a/Show.hs b/Show.hs new file mode 100644 --- /dev/null +++ b/Show.hs @@ -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) + diff --git a/ToDo b/ToDo new file mode 100644 --- /dev/null +++ b/ToDo @@ -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 diff --git a/TreeGrowth.hs b/TreeGrowth.hs new file mode 100644 --- /dev/null +++ b/TreeGrowth.hs @@ -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 diff --git a/Visualize.hs b/Visualize.hs new file mode 100644 --- /dev/null +++ b/Visualize.hs @@ -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") diff --git a/probability.cabal b/probability.cabal new file mode 100644 --- /dev/null +++ b/probability.cabal @@ -0,0 +1,41 @@ +Name: probability +Version: 0.1 +License: BSD3 +Author: Martin Erwig +Maintainer: Henning Thielemann +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