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