exercise44jaynes.py

# -*- coding: utf-8 -*-
import numpy as np
from scipy.special import binom
import random

#Jaynes' PT:TLoS exercise 4.4

def db(x): # Decibel transform
    return 10.0*np.log10(x)

def evidence(x): # Equation 4.8
    return db(1E-100 + x/(1.0 - x + 1E-100))

def b(r, n, f): # Equation 3.86
    return binom(n, r)*np.power(f, r)*np.power(1.0 - f, n - r)
    
# Boxes prior - equation 4.31
priorA = 1.0/11*(1 - 1e-6)
priorB = 10.0/11*(1 - 1e-6)
priorC = 1.0 - priorA - priorB

# Boxes defect rate
boxA = 1.0/3.0
boxB = 1.0/6.0    
boxC = 99.0/100.0
    
def real_posterior(r, n): # Equations 4.33 and 4.39
    plA = priorA*b(r, n, boxA)
    plB = priorB*b(r, n, boxB)    
    plC = priorC*b(r, n, boxC)
    
    eA = evidence(priorA) + db(b(r, n, boxA)*(priorC + priorB)/(plC + plB))
    eB = evidence(priorB) + db(b(r, n, boxB)*(priorA + priorC)/(plA + plC))
    eC = evidence(priorC) + db(b(r, n, boxC)*(priorA + priorB)/(plA + plB))
    
    return eA, eB, eC

realAplot = [evidence(priorA)]
realBplot = [evidence(priorB)]
realCplot = [evidence(priorC)]

Bcount = Blen = Acount = Alen = 0

N = 100000 # Monte Carlo experiment - do as many iterations as you can!

for i in range(N):
    r = n = 0
    while True:
        n += 1
        if random.random() <= boxB: # Suppose hypothesis B is true
            r += 1
        realA, realB, realC = real_posterior(r, n)
        # Stop when 20dB of evidence is accumulated for either A or B
        if realB >= 20:
            Bcount += 1
            Blen += n
            break
        elif realA >= 20:
            Acount += 1
            Alen += n
            break

print "Number of choices A and B:", Acount, Bcount
print "Percentage of choices A and B:", Acount/float(N), Bcount/float(N)
if Acount > 0:
    print "A avg experiment length:", Alen/float(Acount)
if Bcount > 0 :
    print "B avg experiment length:", Blen/float(Bcount)