Initial.

IntegerLineMPD/prob1.lhs

+> import Prelude hiding (seq)
+> 
+> import System.Random
+> import Data.List
+> import Data.Char
+> import Data.Array
+
+1. Probability distributions
+----------------------------
+
+> newtype PD a = PD [(a, Double)] deriving (Eq, Ord, Show)
+
+> instance Functor PD where
+>   fmap f (PD pd) = PD (map (cross(f, id)) pd)
+
+> pd2List :: PD a -> [(a, Double)]
+> pd2List (PD pd) = pd
+
+> probs :: PD a -> [Double]
+> probs (PD aps) = map snd aps
+
+> results  :: PD a -> [a]
+> results (PD aps) = map fst aps
+
+What is the probability of ~a~?  (If ~a~ is not in the support of pd, we "round" it up.)
+
+> get_prob :: (Eq a, Ord a) => a -> PD a -> Double
+> get_prob a (PD pd) = lookup a sorted_nubbed
+>   where
+>   snub [(x, p)] = [(x, p)]
+>   snub ((x1, p1):(x2, p2):rest) = if x1 == x2
+>                                     then snub ((x1, p1+p2):rest)
+>                                     else ((x1,p1):snub ((x2,p2):rest))
+>   sorted_nubbed = snub (sort pd)
+>   lookup a [(_, p)] = p
+>   lookup a ((a', p):rest) = if a <= a' then p else lookup a rest
+
+
+Cummulative probabilities:
+
+> cprobs :: PD a -> [Double]
+> cprobs = scanl1 (+) . probs
+
+> type CPD a = [(a, Double)]
+
+> cPD :: PD a -> CPD a
+> cPD pd = (zip' (pair (results, cprobs) pd)) where zip' = uncurry zip
+
+Expected value:
+
+> expectation :: PD Double -> Double
+> expectation (PD as) = sum (map (uncurry (*)) as) 
+
+Some distributions
+
+> bernoulli :: Double -> PD Double
+> bernoulli p = PD [(0.0, 1-p), (1.0, p)]
+
+> uniform   :: [a] -> PD a
+> uniform as  =  PD (zip as (repeat p))
+>   where p = 1.0 / (toEnum (length as))
+
+> softMax :: [Double] -> PD Double
+> softMax xs = PD (zip xs (map (/tot) exps))
+>   where
+>   exps = map exp xs
+>   tot  = sum exps
+
+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 :: 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 :: Eq a => Nat -> RV a -> RV (PD a)
+> infer n rv = if n == 0 then return (PD []) -- error
+>                 else do ans <- infer' n rv 
+>                         return (freq ans)
+
+> infer' :: Eq a => Nat -> RV a -> RV [(a, Nat)]
+> infer' n rv  = if n == 0
+>                   then return []
+>                   else do a  <- rv
+>                           as <- infer' (n-1) rv
+>                           return (incrementFreq a as)
+>                              
+
+> freq :: [(a, Nat)] -> PD a
+> freq ans = PD [(a, (fromInteger n)/(fromInteger tot)) | (a, n) <- ans]
+>   where
+>   tot = sum (map snd ans)
+
+> incrementFreq :: Eq a => a -> [(a, Nat)] -> [(a, Nat)]
+> incrementFreq a [] = [(a, 1)]
+> incrementFreq a ((a', n) : as) = if a == a'
+>                                     then (a, n+1):as
+>                                     else (a', n):(incrementFreq a as)
+
+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 = pd2List (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 :: RV Double
+> 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
+
+> toss3 :: PD (Double, Double, Double)
+> toss3 = PD [((x, y, z), p) | x <- s, y <- s, z <- s]
+>   where p = 0.5 * 0.5 * 0.5
+>         s = [0.0, 1.0]
+
+> prob :: (a -> Bool) -> PD a -> Double
+> prob e (PD aps) = sum [p | (a, p) <- aps, e a == True]
+
+> cprob :: (a -> Bool) -> (a -> Bool) -> PD a -> Maybe Double
+> cprob 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)
+
+> 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 = cprob fstH twoHs toss3
+
+
+3.2 Factor
+
+From https://agentmodels.org/chapters/3-agents-as-programs.html:
+
+> 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 nn
+>                 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 (zip [Bad ..] [0.2, 0.6, 0.2])
+> t1 French  = PD (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 (fmap 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 Action
+> act s n = act'
+>   where
+>     f a = do e <- expectedUtility s n a
+>              return (100 * e)
+>     inferred = infer numsamp (random_actions >>= f)
+>     act' = do a <- random_actions
+>               e <- f a
+>               accept <- factor inferred e
+>               if accept then return a else act'
+> 
+> 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 a' <- act s' n'
+>           expectedUtility s' n' a'
+>   inferred = infer numsamp re
+
+> simulate s n = do if n == 0
+>                     then return []
+>                     else do a <- act s n
+>                             s' <- return (transition s a)
+>                             ss <- simulate s' (n-1)
+>                             return (s':ss)
+
+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
+>       
+
+
+----------------------------------------------------------------------
+
+Standard combinators
+
+> pair  (f, g) a      = (f a, g a)
+> cross (f, g)        = pair (f . fst, g . snd)
+