HaskCalc/cp2223t/Probability.hs
2022-11-19 11:35:14 +00:00

697 lines
17 KiB
Haskell

-- Credits: Erwig, Martin and Kollmannsberger, Steve
-- FUNCTIONAL PEARLS: Probabilistic functional programming in Haskell,
-- JFP, 2006
-- DOI: 10.1017/S0956796805005721
module Probability where
import qualified System.Random
import Data.List (sort,sortBy,transpose)
import Control.Monad
import System.IO.Unsafe (unsafePerformIO)
import ListUtils
import Show
{- TO DO:
* create export list
* extend Dist by a constructor for continuous distributions:
C (Float -> Float)
* prove correctness of |||
-}
--------- jno (to be checked) --------------
import Control.Applicative
instance Applicative Dist where
pure = return
(<*>) = ap
instance Alternative Dist where
empty = D []
(<|>) = (>>)
--------- end jno --------------------------
------------------------------------------------------------------------------
-- CONTENTS:
--
-- 0 AUXILIARY DEFINITIONS
-- 1 DETERMINISTIC AND PROBABILISTIC VALUES
-- 2 RANDOMIZED VALUES
-- 3 DETERMINISTIC AND PROBABILISTIC GENERATORS
-- 4 RANDOMIZED GENERATORS
-- 5 ITERATORS AND SIMULATORS
-- 6 TRACING
------------------------------------------------------------------------------
------------------------------------------------------------------------------
-- 0 AUXILIARY DEFINITIONS
--
-- Event
-- Probability
------------------------------------------------------------------------------
--
-- 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 | precision==0 = showR 3 (round (f*100))++"%"
| otherwise = showR (4+precision) (fromIntegral (round (f*100*d))/d)++"%"
where d = 10^precision
-- -- 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
-- sequ composition of a list of monadic functions
--
(>@>) :: Monad m => (a -> m b) -> (b -> m c) -> a -> m c
f >@> g = (>>= g) . f
sequ :: Monad m => [a -> m a] -> a -> m a
sequ = foldl (>@>) return
------------------------------------------------------------------------------
-- 1 DETERMINISTIC AND PROBABILISTIC VALUES
--
-- Dist probability disribution
-- Spread functions to convert a list of values into a distribution
------------------------------------------------------------------------------
--
-- Distributions
--
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)]
-- note: mzero is a zero for >>= and a unit for mplus
--
instance MonadPlus Dist where
mzero = D []
mplus d d' | isZero d || isZero d' = mzero
| otherwise = unfoldD $ choose 0.5 d d'
isZero :: Dist a -> Bool
isZero (D d) = null d
instance Functor Dist where
fmap f (D d) = D [(f x,p) | (x,p) <- d]
instance (Ord a,Eq a) => Eq (Dist a) where
D xs == D ys = map fst (norm' xs)==map fst (norm' ys) &&
all (\((_,p),(_,q))->abs (p-q)<errorMargin) (zip (norm' xs) (norm' ys))
-- auxiliary functions for constructing and working with distributions
--
onD :: ([(a,ProbRep)] -> [(a,ProbRep)]) -> Dist a -> Dist a
onD f = D . f . unD
sizeD :: Dist a -> Int
sizeD = length . unD
checkD :: Dist a -> Dist a
checkD (D d) | abs (1-sumP d) < errorMargin = D d
| otherwise = error ("Illegal distribution: total probability = "++show (sumP d))
mkD :: [(a,ProbRep)] -> Dist a
mkD = checkD . D
sumP :: [(a,ProbRep)] -> ProbRep
sumP = sum . map snd
sortP :: [(a,ProbRep)] -> [(a,ProbRep)]
sortP = sortBy (\x y->compare (snd y) (snd x))
-- normalization = grouping
--
normBy :: Ord a => (a -> a -> Bool) -> Dist a -> Dist a
normBy f = onD $ accumBy f . sort
accumBy :: Num b => (a -> a -> Bool) -> [(a,b)] -> [(a,b)]
accumBy f ((x,p):ys@((y,q):xs)) | f x y = accumBy f ((x,p+q):xs)
| otherwise = (x,p):accumBy f ys
accumBy _ xs = xs
norm :: Ord a => Dist a -> Dist a
norm = normBy (==)
norm' :: Ord a => [(a,ProbRep)] -> [(a,ProbRep)]
norm' = accumBy (==) . sort
-- pretty printing
--
instance (Ord a,Show a) => Show (Dist a) where
show (D []) = "Impossible"
show (D xs) = concatMap (\(x,p)->showR w x++' ':showP p++"\n") (sortP (norm' xs))
where w = maximum (map (length.show.fst) xs)
--
-- Operations on distributions
--
-- product of independent distributions
--
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 (,)
-- 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 stddev x = 1 / sqrt (2 * pi) * exp (-1/2 * u^2)
where u = (x - mean) / stddev
-- 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)
-- 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) | p<=q || null ps = x
| otherwise = scanP (p-q) ps
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
instance Expected a => Expected (IO a) where
expected r = expected (System.IO.Unsafe.unsafePerformIO r)
-- 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
------------------------------------------------------------------------------
-- 2 RANDOMIZED VALUES
--
-- R random value
-- RDist random distribution
------------------------------------------------------------------------------
--
-- 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 = System.Random.randomRIO (0,1) >>= return . selectP d
--
-- 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)]))
------------------------------------------------------------------------------
-- 3 DETERMINISTIC AND PROBABILISTIC GENERATORS
--
-- Change deterministic generator
-- Trans probabilistic generator
-- SpreadC functions to convert a list of changes into a transition
-- SpreadT functions to convert a list of transitions into a transition
------------------------------------------------------------------------------
--
-- transitions
--
type Change a = a -> a
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
--
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 = liftC uniform
normalT = liftC normal
linearT c = liftC (linear c)
enumT xs = liftC (enum xs)
-- spreading transitions into transitions
--
type SpreadT a = [Trans a] -> Trans a
liftT :: Spread (Trans a) -> [Trans a] -> Trans a
liftT s = unfoldT . s
uniformTT = liftT uniform
normalTT = liftT normal
linearTT c = liftT (linear c)
enumTT xs = liftT (enum xs)
------------------------------------------------------------------------------
-- 4 RANDOMIZED GENERATORS
--
-- RChange random change
-- RTrans random transition
------------------------------------------------------------------------------
--
-- Randomized changes
--
type RChange a = a -> R a
random :: Trans a -> RChange a
random t = pick . t
-- random = (pick .)
--
-- Randomized transitions
--
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
------------------------------------------------------------------------------
-- 5 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
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 contains two functions:
--
-- ~. returns the final randomized transition
-- ~.. returns the whole trace
--
-- 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
(~.) :: Ord a => Int -> (a -> c a) -> RTrans a
(~..) :: 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 =
------------------------------------------------------------------------------
-- 7 TRACING
--
-- (R)Trace
-- (R)Space
-- (R)Walk
-- (R)Expand
------------------------------------------------------------------------------
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 RTraces, 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
{-
LAWS
const . pick = random . const
-}