Playing with Probabilities

Posted on April 7, 2014

The Header

Imports, pragmas, and the module declaration

> {-# LANGUAGE TupleSections #-}
> module MDists where
> import Data.Map (Map, (!))
> import qualified Data.Map as M
> import Data.Function
> import Data.List
> import Data.Ratio

The Code

> type P = Rational
> type PDist a = Map a P

Hypotheses are getting represented as probability distributions with a probability attached to them

> data Hypothesis a = H {conditionals :: PDist a,
>                          chance :: P} deriving (Show)

Next, the distribution over hypotheses - let’s make strings our labels, just because strings are easy to use

> data HDist a = HD {hps :: Map String (Hypothesis a)} deriving (Show)

Helpers

First, we want to implement a few basic funcitons that will allow the rest of these to work nicely For example, we’ll want to be able to apply some function [0, 1] → [0, 1] to the probabilities we’re juggling here

> applyToProbability :: (P -> P) -> Hypothesis a -> Hypothesis a
> applyToProbability f (H d p) = H d $ f p

That one’s pretty simple, but on its own not that useful. Instead, let’s wrap it up so we don’t have to deal with the individual hypotheses - instead, we’ll map it over a the hypothesis distribution

> applyToProbs :: (P -> P) -> HDist a -> HDist a
> applyToProbs f (HD hps) = HD $ M.map (applyToProbability f) hps

Now, the main reason we’re defining those two functions is to write a function renormalize that takes a distribution and makes sure the probabilities actually sum to 1 without changing any of the ratios between probabilities

> renormalize :: HDist a -> HDist a
> renormalize h@(HD dist) = applyToProbs (/normalizer) h
>   where normalizer = M.foldl' (+) 0 . M.map chance $ dist

Our basic guess function is just a small wrapper around M.lookup - it fetches the likelihood of a hypothesis on its identifier If we don’t have a hypothesis in our distribution, it assigns it a probability of 0, but otherwise just fetches the probability from the hypothesis

> guess :: Ord a => HDist a -> String -> P
> guess (HD hps) str = case M.lookup str hps of Nothing -> 0
>                                               Just h -> chance h

Updating!

Now, our not-so-secret goal here was to model Bayesian updates If you’re not familiar with Bayes’ Theorem, here it is:


$$p(H|D)=p(D|H)\frac{p(H)}{p(D)}$$

To explain: the probability of a hypothesis being true, on some data, is the probability of that data given the hypothesis (p(DH)) multiplied by the probability we had previously assigned the hypothesis (our “prior”, or p(H)), then divided by the probability of that data under all hypotheses (p(D)). We don’t deal directly with the probability of the data under all hypotheses - instead preferring to renormalize the distribution after the fact

Of course, we want to automate this adjustment. The first step is to multiply each probability by the numerator - p(DH) ⋅ p(H)

> pAndCond :: Ord a => a -> Hypothesis a -> Hypothesis a
> pAndCond event (H dist prior) = H dist $ prior * likelihood
>   where likelihood = M.findWithDefault 0 event dist

With that done, we just renormalize, effecitvely dividing by p(D) We know this works because if we divided by anything other than the renormalizer, the ending probabilities would sum to some number other than 1, so, on pain of contradiction,
p(D) = renormalizer

> updateOnEvent :: Ord a => HDist a -> a -> HDist a
> updateOnEvent (HD dist) event = renormalize . HD . M.map (pAndCond event) $ dist

Of course, we want to be able to update on a series of events, so we fold over a list of events This fold ends up being fully strict, so instead of a foldr we use foldl’ (hooray tail recursion!)

> update :: Ord b => HDist b -> [b] -> HDist b
> update = foldl' updateOnEvent

Construction helpers

We’ll use these utilities to construct and solve some toy problems a little further down First, a few more helpers that’ll be really important for constructing the toy problems

This one just normalizes a hypothesis, so we can think of it in terms of allocating odds instead of probabilities

> normalizeHyp :: Hypothesis a -> Hypothesis a
> normalizeHyp (H m p) = if totalP == 1 then (H m p) else H (M.map (/totalP) m) p
>   where totalP = M.foldl' (+) 0 m

Next, this one just makes a hypothesis distribution from a list

> fromList :: [(String, Hypothesis a)] -> HDist a
> fromList = renormalize . HD . M.fromList

Toy Problems

Monty Hall

We’ll start out by getting a list of doors

> drs = "abc"

Next, our function for generating a hypothesis for which door gets opened based on our guess and the correct answer

> isIn :: Char -> Char -> Hypothesis Char
> isIn guess correct = normalizeHyp . flip H 1 . M.fromList $
>                      [(dr,1) | dr <- drs, dr /= guess, dr /= correct]

And, finally, the problem itself! The identifiers are each just the single character name of the door (a, b, or c) The hypotheses are each generated assuming the identifier is the correct door according to the isIn function

> montyHall :: Char -> HDist Char
> montyHall guess = fromList [([dr],isIn guess dr) | dr <- drs]
> shouldYouSwitch = guess (update (montyHall 'a') "b") "a" < 1 % 2

And there you go! The boolean shouldYouSwitch starts from you picking door a, then Monty opens b, and checks if you have worse than even chances of it being behind a

Cookies

We’ve got two bools of chocolate and vanilla cookies The first bowl is three quarters vanilla and one quarter chocolate, while the other is half and half

> bowl1 = normalizeHyp $ flip H 1 $ M.fromList [('v',3),('c',1)]
> bowl2 = normalizeHyp $ flip H 1 $ M.fromList [(a,1) | a <- "vc"]

Next we’ll make a distribution from both the bowls

> bowls = fromList  [("b1",bowl1),("b2",bowl2)]

From this, given a stream of cookies, it’ll guess the chance you were drawing from one of the bowls

M&M’s

Some background: The mixture of M&M’s has been changed a few times In 1995, blue M&M’s were introduced Given a bundle of M&M’s, what’s the probability it came from a 1994 mix vs a 1996 mix?

color percentage in ’94 percentage in ’96
brown 30 24
green 10 20
orange 10 16
yellow 20 14
red 20 13
tan 10 0
blue 0 13

Now, let’s encode that into hypotheses

> mix94 = normalizeHyp . flip H 1 $ M.fromList
>         [("brown",30),("green",10),("orange",10),("yellow",20),("red",20),("tan",10)]
> mix96 = normalizeHyp . flip H 1 $ M.fromList
>         [("brown",24),("green",20),("orange",16),("yellow",14),("red",13),("blue",13)]

And, finally, our hypothesis distribution

> bag = fromList [("94",mix94),("96",mix96)]