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