{-# LANGUAGE GeneralizedNewtypeDeriving, TypeSynonymInstances, FlexibleInstances #-}
module Probability (module Perhaps,
PerhapsT(..),
Dist, weighted, uniform, certainly, (??)
) where
import Perhaps
-- All of these are exclusively used to print
-- histograms of probability distributions
import Text.PrettyPrint
import Data.List (groupBy, sortBy)
import Data.Function (on)
-----------------------------------------------------------------------
-- Monad transformer version of monad Perhaps:
-- any monad m layered on top of the probability monad Perhaps
newtype PerhapsT m a = PerhapsT { runPerhapsT :: m (Perhaps a) }
instance Monad m => Monad (PerhapsT m) where
return x = PerhapsT $ return (return x)
mo >>= g = PerhapsT $ do
ph <- runPerhapsT mo
case ph of
Perhaps _ 0 -> return never
Perhaps x p1 -> do Perhaps y p2 <- runPerhapsT (g x)
return $ Perhaps y (p1 * p2)
-----------------------------------------------------------------------
-- Probability distribution:
-- A list of values with associated probability
type Dist = PerhapsT []
-- Rendering a probability distribution (as a histogram), e.g.
--
-- 2 2.8% ##
-- 3 5.6% ####
-- 4 8.3% ######
-- 5 11.1% ########
-- 6 13.9% ##########
-- 7 16.7% ############
-- 8 13.9% ##########
-- 9 11.1% ########
-- 10 8.3% ######
-- 11 5.6% ####
-- 12 2.8% ##
instance (Show a, Ord a) => Show (Dist a) where
show x = show $ histogram $ merge $ runPerhapsT x
where
-- Merge entries in the distribution for the same value
merge :: (Ord a) => [Perhaps a] -> [Perhaps a]
merge = map (foldr1 add) .
groupBy ((==) `on` value) .
sortBy (compare `on` value)
add (Perhaps v1 p1) (Perhaps _ p2) = Perhaps v1 (p1 + p2)
histogram :: (Show a) => [Perhaps a] -> Doc
histogram phs = vcat . map bucket $ phs
where
bucket (Perhaps v p) =
hsep [pad width (show v),
pad 6 (show p),
hcat (replicate (round (p * height)) (char '#'))]
width = maximum . map (length . show . value) $ phs
height = 70
pad w s = text $ replicate (w - length s) ' ' ++ s
-- A few useful combinators over probability distributions
weighted :: [(a, Prob)] -> Dist a
weighted [] = error "empty distribution"
weighted xws = PerhapsT $ map weight xws
where
weight (v, w) = Perhaps v (w / total)
total = sum . map snd $ xws
uniform :: [a] -> Dist a
uniform xs = weighted [ (x, 1.0) | x <- xs ]
-- uniform xs = weighted $ zip xs (repeat 1.0)
certainly :: a -> Dist a
certainly v = uniform [v]
-- certainly = return
infix 8 ??
-- Probability that predicate (event) q is satisfied by
-- the values in distribution ds
(??) :: (a -> Bool) -> Dist a -> Prob
q ?? ds = sum [ p | Perhaps v p <- runPerhapsT ds, q v ]