Bayesian update

In this notebook we experiment with the Bayesian update of the "probability of a probability" from the number of successes in repeated independent experiments.

We start from a Beta prior, which can model many different kinds of previous knowledge about a probability. For example:





 Beta(3,2) 


































 

 






















 0.1 
 0.2 
 0.3 
 0.4 
 0.5 
 0.6 
 0.7 
 0.8 
 0.9 

 p 




 0.0 
 0.2 
 0.4 
 0.6 
 0.8 
 1.0 
 1.2 
 1.4 
 1.6 
 1.8 







 P

Using the tools described in this post and the "grid" approximation used by Kruschke we can model the experiment as follows:

exper = joint (coins 5) prior

coins n p = replicateM n (coin p)

coin p = bernoulli p 'H' 'T' -- weighted [('H',p),('T',1-p)]

prior = discretize nT (0+deltaHist,1-deltaHist) (betaD 3 2)
 exper

  0.00 % ("HHHHH","0.01")
  0.00 % ("HHHHH","0.02")
  0.00 % ("HHHHH","0.03")
   ⋮
  0.00 % ("TTTTT","0.98")
  0.00 % ("TTTTT","0.99")

We are actually interested in the total number of heads:

exper2 = (length . filter (=='H') *** id) <-- exper
 exper2

  0.00 % (0,"0.01")
  0.00 % (0,"0.02")
  0.01 % (0,"0.03")
   ⋮
  0.21 % (5,"0.98")
  0.11 % (5,"0.99")

We can now compute the posterior distribution given the observed number of heads:

 snd <-- exper2 ! (0==) . fst

  0.02 % "0.01"
  0.09 % "0.02"
  0.19 % "0.03"
   ⋮
  0.00 % "0.98"
  0.00 % "0.99"

We graphically show the update process for different priors and observations:





 prior beta(3.0,2.0), 0 heads in 5 coins 







































 

 












































 0.1 
 0.2 
 0.3 
 0.4 
 0.5 
 0.6 
 0.7 
 0.8 
 0.9 

 p 




 0.0 
 0.2 
 0.4 
 0.6 
 0.8 
 1.0 
 1.2 
 1.4 
 1.6 
 1.8 
 2.0 
 2.2 
 2.4 
 2.6 
 2.8 







 P 















































 0.06 




 0.59 











 prior beta(3.0,2.0), 5 heads in 5 coins 




























 

 












































 0.1 
 0.2 
 0.3 
 0.4 
 0.5 
 0.6 
 0.7 
 0.8 
 0.9 

 p 




 0 
 1 
 2 
 3 







 P 




































 0.56 




 1.00 











 prior beta(1.0,1.0), 4 heads in 5 coins 





































 

 












































 0.1 
 0.2 
 0.3 
 0.4 
 0.5 
 0.6 
 0.7 
 0.8 
 0.9 

 p 




 0.0 
 0.2 
 0.4 
 0.6 
 0.8 
 1.0 
 1.2 
 1.4 
 1.6 
 1.8 
 2.0 
 2.2 
 2.4 







 P 













































 0.40 




 0.98 











 prior beta(1.0,1.0), 5 heads in 10 coins 







































 

 












































 0.1 
 0.2 
 0.3 
 0.4 
 0.5 
 0.6 
 0.7 
 0.8 
 0.9 

 p 




 0.0 
 0.2 
 0.4 
 0.6 
 0.8 
 1.0 
 1.2 
 1.4 
 1.6 
 1.8 
 2.0 
 2.2 
 2.4 
 2.6 
 2.8 







 P 















































 0.22 




 0.78 











 prior beta(0.4,0.4), 4 heads in 6 coins 





























 

 












































 0.1 
 0.2 
 0.3 
 0.4 
 0.5 
 0.6 
 0.7 
 0.8 
 0.9 

 p 




 0 
 1 
 2 
 3 
 4 







 P 





































 0.30 




 0.96

Obviously, for a larger number of coins there are too many outcomes in exper and this exhaustive approach breaks down. Alternately, we can directly compute the distribution of the number of heads using the binomial distribution.

(In fact, this problem has a beautiful analytic solution since the Beta and binomial are conjugate distributions, but we are interested here in a general computational approach valid for other priors and observation models.)

Using this shortcut the result is exactly the same as above and we can handle bigger experiments:





 prior beta(0.4,0.4), 4 heads in 6 coins 





























 

 












































 0.1 
 0.2 
 0.3 
 0.4 
 0.5 
 0.6 
 0.7 
 0.8 
 0.9 

 p 




 0 
 1 
 2 
 3 
 4 







 P 





































 0.30 




 0.96 











 prior beta(2.0,2.0), 150 heads in 200 coins 































 

 












































 0.1 
 0.2 
 0.3 
 0.4 
 0.5 
 0.6 
 0.7 
 0.8 
 0.9 

 p 




 0 
 2 
 4 
 6 
 8 
 10 
 12 







 P 







































 0.68 




 0.80 











 prior beta(1.0,1.0), 0 heads in 10 coins 






























 

 












































 0.1 
 0.2 
 0.3 
 0.4 
 0.5 
 0.6 
 0.7 
 0.8 
 0.9 

 p 




 0 
 2 
 4 
 6 
 8 
 10 







 P 






































 0.00 




 0.26 











 prior beta(1.0,1.0), 0 heads in 100 coins 































 

 












































 0.1 
 0.2 
 0.3 
 0.4 
 0.5 
 0.6 
 0.7 
 0.8 
 0.9 

 p 




 0 
 10 
 20 
 30 
 40 
 50 
 60 







 P 







































 0.00 




 0.03

You can try other parameters: