"""This file contains code for use with "Think Bayes", by Allen B. Downey, available from greenteapress.com Copyright 2012 Allen B. Downey License: GNU GPLv3 http://www.gnu.org/licenses/gpl.html """ """This file contains a partial solution to a problem from MacKay, "Information Theory, Inference, and Learning Algorithms." Exercise 3.15 (page 50): A statistical statement appeared in "The Guardian" on Friday January 4, 2002: When spun on edge 250 times, a Belgian one-euro coin came up heads 140 times and tails 110. 'It looks very suspicious to me,' said Barry Blight, a statistics lecturer at the London School of Economics. 'If the coin were unbiased, the chance of getting a result as extreme as that would be less than 7%.' MacKay asks, "But do these data give evidence that the coin is biased rather than fair?" """ import thinkbayes class Euro(thinkbayes.Suite): """Represents hypotheses about the probability of heads.""" def Likelihood(self, data, hypo): """Computes the likelihood of the data under the hypothesis. hypo: integer value of x, the probability of heads (0-100) data: tuple of (number of heads, number of tails) """ x = hypo / 100.0 heads, tails = data like = x**heads * (1-x)**tails return like def TrianglePrior(): """Makes a Suite with a triangular prior.""" suite = Euro(xrange(0, 101)) for x in range(0, 51): suite.Set(x, x) for x in range(51, 101): suite.Set(x, 100-x) suite.Normalize() return suite def SuiteLikelihood(suite, data): """Computes the weighted average of likelihoods for sub-hypotheses. suite: Suite that maps sub-hypotheses to probability data: some representation of the data returns: float likelihood """ total = 0 for hypo, prob in suite.Items(): like = suite.Likelihood(data, hypo) total += prob * like return total def Main(): data = 140, 110 data = 8, 12 suite = Euro() like_f = suite.Likelihood(data, 50) print 'p(D|F)', like_f actual_percent = 100.0 * 140 / 250 likelihood = suite.Likelihood(data, actual_percent) print 'p(D|B_cheat)', likelihood print 'p(D|B_cheat) / p(D|F)', likelihood / like_f like40 = suite.Likelihood(data, 40) like60 = suite.Likelihood(data, 60) likelihood = 0.5 * like40 + 0.5 * like60 print 'p(D|B_two)', likelihood print 'p(D|B_two) / p(D|F)', likelihood / like_f b_uniform = Euro(xrange(0, 101)) b_uniform.Remove(50) b_uniform.Normalize() likelihood = SuiteLikelihood(b_uniform, data) print 'p(D|B_uniform)', likelihood print 'p(D|B_uniform) / p(D|F)', likelihood / like_f b_tri = TrianglePrior() b_tri.Remove(50) b_tri.Normalize() likelihood = b_tri.Update(data) print 'p(D|B_tri)', likelihood print 'p(D|B_tri) / p(D|F)', likelihood / like_f if __name__ == '__main__': Main()