Beginning of version with Data.Map for PD.

IntegerLineMPD/prob_map.lhs

+> import Prelude hiding (seq)
+> 
+> import System.Random
+
+> import qualified Data.Map.Strict as Map
+
+> import Data.List
+
+ import Data.Char
+ import Data.Array
+
+----------------------------------------------------------------------
+
+Standard combinators
+
+> pair  (f, g) a      = (f a, g a)
+> cross (f, g)        = pair (f . fst, g . snd)
+
+
+1. Probability distributions
+----------------------------
+
+> newtype PD a = PD (Map.Map a Double) deriving (Eq, Ord, Show)
+
+> fmapPD :: Ord b => (a -> b) -> PD a -> PD b
+> fmapPD f (PD pd) = PD (Map.mapKeysWith (+) f pd)
+
+> pd2Map :: PD a -> Map.Map a Double
+> pd2Map (PD pd) = pd
+
+> probs :: PD pd -> [Double]
+> probs (PD pd) = Map.elems pd
+
+> results  :: PD a -> [a]
+> results (PD pd) = Map.keys pd
+
+What is the probability of ~a~?  (If ~a~ is not in the support of pd, we "round" it up.)
+
+> get_prob a (PD pd) = if g == Nothing
+>                         then get_val s 
+>                         else get_val g
+>   where
+>   s = Map.lookupLE a pd
+>   g = Map.lookupGE a pd
+>   fromMaybe (Just x) = x
+>   get_val = snd . fromMaybe
+
+Some distributions
+
+> bernoulli :: Double -> PD Double
+> bernoulli p = PD (Map.fromList [(0.0, 1-p), (1.0, p)])
+
+> uniform   :: (Ord a) => [a] -> PD a
+> uniform as  =  PD (Map.fromList (zip as (repeat p)))
+>   where p = 1.0 / (toEnum (length as))
+
+> softMax :: [Double] -> PD Double
+> softMax xs = PD (Map.fromList (zip xs (map (/tot) exps)))
+>   where
+>   exps = map exp xs
+>   tot  = sum exps
+
+Expected value:
+
+> expectation :: PD Double -> Double
+> expectation (PD pd) = Map.foldlWithKey f 0 pd
+>   where
+>   f sum x p = sum + x * p
+
+
+Cumulative probabilities:
+
+> cprobs = scanl1 (+) . probs
+
+ > type CPD a = Map.Map a Double
+
+> cPD :: Ord a => PD a -> [(a, Double)]
+> cPD pd = zip (results pd) (cprobs pd)
+
+
+> {--
+
+Probability of an event:
+
+> type Event a = a -> Bool
+
+> prob :: Event a -> PD a -> Double
+> prob e (PD aps) = sum [p | (a, p) <- aps, e a == True]
+
+Conditional probability:
+
+> condProb :: Event a -> Event a -> PD a -> Maybe Double
+> condProb e1 e2 pd = let pe2 = prob e2 pd in
+>                  if pe2 == 0.0
+>                  then Nothing
+>                  else Just (prob (\ a -> e1 a && e2 a) pd / pe2)
+
+
+> --}
+
+2. Random values (sampling)
+---------------------------
+
+> type RNG = StdGen
+> 
+> rng0 :: RNG
+> rng0 = mkStdGen 123321
+
+The random value datatype:
+
+> newtype RV a = RV (RNG -> (a, RNG))
+>
+> app :: RV a -> RNG -> (a, RNG)
+> app (RV rv) rng = rv rng
+
+Monad implementation:
+
+> class Monad1 m where
+>   eta :: a -> m a
+>   seq :: (a -> m b) -> (b -> m c) -> a -> m c
+>   bind :: m a -> (a -> m b) -> m b
+>   bind ma f = ((const ma) `seq` f)()
+>   mu  :: m (m a) -> m a
+>   mu mma = mma `bind` id
+>
+> instance Monad1 RV where
+>   eta a = RV (\ rng -> (a, rng))
+>   seq f g a = RV (uncurry app . cross (g, id) . app (f a))
+
+> instance Functor RV where
+>   fmap f rva = rva `bind` (eta . f)
+
+> instance Applicative RV where
+>   pure = eta
+>   rvf <*> rv = rvf `bind` (\ f -> fmap f rv)
+>
+> instance Monad RV where
+>   return = eta
+>   (>>=)  = bind
+
+3. Sampling and inference
+-------------------------
+
+Roulette sampling:
+
+> sample :: Ord a => PD a -> RV a
+> sample pd = RV rv
+>   where
+>   rv rng = let (p, rng') = random rng in 
+>            let acs = dropWhile (\ (a, cp) -> cp < p) (cPD pd) in
+>            (fst (head acs), rng')
+
+Example: tossing a biased coin
+
+> toss :: Double -> RV Double
+> toss p = sample (bernoulli p)
+
+
+Infer: construct a PD by sampling a given number of times (should be strictly positive)
+
+> type Nat = Integer
+
+> infer :: Ord a => Nat -> RV a -> RV (PD a)
+> infer n rv = do ans <- infer' n rv (Map.empty)
+>                 return (freq ans)
+
+> infer' :: Ord a => Nat -> RV a -> Map.Map a Nat -> RV (Map.Map a Nat)
+> infer' n rv m  = if n == 0
+>                   then return m
+>                   else do a  <- rv
+>                           infer' (n-1) rv (Map.alter inc a m)
+>   where
+>   inc Nothing = Just 1
+>   inc (Just seenSoFar) = Just (seenSoFar + 1)
+
+
+> freq :: Map.Map a Nat -> PD a
+> freq ans = PD (Map.map f ans)
+>   where
+>   tot = sum (Map.elems ans)
+>   f n = (fromInteger n)/(fromInteger tot)
+
+
+The relationship between "sample" and "infer" (and RNG) is
+
+a) "Frequentist" semantics: sampling frequencies approximate the pd
+
+    infer bignum (sample pd) =?= return pd
+
+b) "Bayesian" semantics: any random variable comes from a pd
+
+   forall rv: exists pd: 
+    forall rng: fst ((infer bignum rv rng)) ~ pd
+
+Nicola: these requirements are actually a specification for the number generator.
+
+Alternative requirements: vary the random number generator in "infer".
+
+Helper function: Get the probability distribution representing a random variable
+
+> get_pd numsamples rv = fst (app (infer numsamples rv) rng0)
+
+
+Example: Monty Hall
+
+> doors = ['A', 'B', 'C']
+> monty_hall switch = do car <- sample (uniform doors)
+>                        guess <- sample (uniform doors)
+>                        hint <- sample(uniform (doors \\ [car, guess]))
+>                        if switch then return (car ==  head (doors \\ [guess, hint]))
+>                                  else return (car == guess)
+
+> keep_strategy_results = get_pd 100 (monty_hall False)
+> switch_strategy_results = get_pd 100 (monty_hall True)
+
+
+3.1 Conditioning
+
+What is the probability of the first flip being a head given that there are at least two heads?
+
+> twoOutOfThree =
+>   do a <- toss 0.5
+>      b <- toss 0.5
+>      c <- toss 0.5
+>      if a + b + c >= 2.0
+>         then return a 
+>         else twoOutOfThree
+
+> {--
+
+
+> pdProd :: PD a -> PD b -> PD (a, b)
+> pdProd (PD aps) (PD bps) = PD [((a, b), pa*pb) | (a, pa) <- aps, (b, pb) <- bps]
+
+Nicola:
+
+> toss3 :: PD (Double, Double, Double)
+> toss3 = arrange (pdProd coin (pdProd coin coin))
+>   where coin = bernoulli 0.5
+>         arrange = fmap f
+>         f (x, (y, z)) = (x, y, z)
+
+
+> fstH :: (Double, Double, Double) -> Bool
+> fstH (x, _, _) = x == 1.0
+
+> twoHs :: (Double, Double, Double) -> Bool
+> twoHs (x, y, z) = x + y + z >= 2.0
+
+> lala = condProb fstH twoHs toss3          
+
+> --}
+
+3.2 Factor
+----------
+
+From https://agentmodels.org/chapters/3-agents-as-programs.html:
+
+> factor :: RV (PD Double) -> Double -> RV Bool
+> factor inferred v = do pd <- inferred
+>                        r <- (RV random)
+>                        if r <= get_prob  v (softMax (results pd))
+>                            then return True
+>                            else return False
+
+Examples:
+
+a) "Prefer" tosses where there are more heads:
+
+> softHeads = softHeads'
+>   where
+>     softHeads' = do (a, b, c) <- rabc
+>                     s <- sumabc (a, b, c)
+>                     accept <- factor inferred s
+>                     if accept
+>                       then return a
+>                       else softHeads'
+>                      
+>     rabc = do x <- toss 0.5
+>               y <- toss 0.5
+>               z <- toss 0.5
+>               return (x, y, z)
+>     rsum = rabc >>= sumabc
+>     sumabc (x, y, z) = return (x+y+z)
+>     inferred = infer 100 rsum
+
+
+b) "Prefer" numbers proportionally to their square:
+
+> nsqr  = nsqr'
+>   where
+>   rn       = sample (uniform [0, 1, 2])
+>   sqr_n n  = return (n * n)
+>   inferred = infer 100 (rn >>= sqr_n)
+>   nsqr'    = do n <- rn
+>                 -- nn <- sqr_n n
+>                 accept <- factor inferred (n * n)
+>                 if accept
+>                   then return n
+>                   else nsqr'
+
+
+4. Applications
+---------------
+
+4.1 One-shot decisions
+
+- S states
+- A actions
+- transition : A -> PD S
+- utility : S -> Double
+
+Choose "a" that maximizes expected utility (https://agentmodels.org/chapters/3-agents-as-programs.html)
+
+> data A1 = Italian | French deriving (Eq, Ord, Show, Enum)
+> data S1 = Bad | Good | Spectacular deriving (Eq, Ord, Show, Enum)
+
+> t1 Italian = PD (Map.fromList (zip [Bad ..] [0.2, 0.6, 0.2]))
+> t1 French  = PD (Map.fromList (zip [Bad ..] [0.05, 0.9, 0.05]))
+
+> u1 Bad = -10
+> u1 Good = 6
+> u1 Spectacular = 8
+
+> softMaxAgent = softMaxAgent'
+>   where
+>   random_action = sample (uniform [Italian ..])
+>   f a = return (expectation (fmapPD u1 (t1 a)))
+>   inferred = infer 100 (random_action >>= f)
+>   softMaxAgent' = do a <- random_action
+>                      e <- f a
+>                      accept <- factor inferred e
+>                      if accept then return a else softMaxAgent'
+
+4.2 Integer Line MDP
+
+From https://agentmodels.org/chapters/3a-mdp.html.
+
+> type State  = Int
+> type Action = Int
+
+> actions = [-1, 0, 1]
+
+> transition :: State -> Action -> State
+> transition = (+)
+
+> utility :: State -> Double
+> utility s = if s == 3 then 1.0 else 0.0
+
+> random_actions = sample (uniform actions)
+
+> numsamp = 5
+
+> act :: State -> Int -> RV (PD Action)
+> act s n = infer numsamp ra
+>   where
+>   f a = do e <- expectedUtility s n a
+>            return (100 * e)
+>   inferred = infer numsamp (random_actions >>= f)
+>   ra = do a <- random_actions
+>           e <- f a
+>           accept <- factor inferred e
+>           if accept then return a else ra
+
+> expectedUtility s n a = if n == 0
+>                           then return u
+>                           else do rexp <- inferred
+>                                   return (u + expectation rexp)
+>   where
+>   u = utility s
+>   n' = n - 1
+>   s' = transition s a
+>   re = do pa' <- act s' n'
+>           a'  <- sample pa'
+>           expectedUtility s' n' a'
+>   inferred = infer numsamp re
+ 
+> simulate s n = do if n == 0
+>                     then return []
+>                     else do pa <- act s n
+>                             a  <- sample pa
+>                             s' <- return (transition s a)
+>                             ss <- simulate s' (n-1)
+>                             return (s':ss)
+
+Without memoisation, we cannot compute, e.g,
+
+  λ> app (act 2 3) rng0
+ (PD (fromList *** Exception: stack overflow
+
+
+> {--
+
+
+
+
+Using memoisation:
+
+("tabulate" taken from ADWH, p. 314)
+
+> tabulate :: Ix i => (i -> e) -> (i, i) -> Array i e
+> tabulate f bounds = array bounds [(x, f x) | x <- range bounds]
+
+> simMemo state steps = sim' state steps
+>   where
+>   sim' s n = do if n == 0
+>                   then return []
+>                   else do a <- actArray!(s, n)
+>                           s' <- return (transition s a)
+>                           ss <- sim' s' (n-1)
+>                           return (s':ss)
+>   actArray = tabulate act ((state-steps, 0), (state+steps, steps))
+>   act (s, n) = do a <- random_actions
+>                   e <- f a
+>                   accept <- factor inferred e
+>                   if accept then return a else actArray!(s, n)
+>     where
+>       f a = do e <- expUtilArray!(s, n, a)
+>                return (100*e)
+>       inferred = infer numsamp (random_actions >>= f)
+>   expUtilArray = tabulate expUtil ((state-steps, 0, -1), (state+steps, steps, 1))
+>   expUtil (s, n, a) = if n == 0
+>                         then return u
+>                         else do rexp <- inferred
+>                                 return (u + expectation rexp)
+>     where
+>       u = utility s
+>       n' = n - 1
+>       s' = transition s a
+>       re = do a' <- actArray!(s', n')
+>               expUtilArray!(s', n', a')
+>       inferred = infer numsamp re
+>       
+
+
+
+> --}