exercise42jaynes.py

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

#Jaynes' PT:TLoS exercise 4.2

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 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

F = []

for x in range(100):
    F.append([])
    priorC = x/99.0
    priorA = 1.0/11*(1.0 - priorC)
    priorB = 10.0/11*(1.0 - priorC)
    for y in range(100):
        boxC = y/99.0
        F[x].append(0)
        for f in range(100):
            A, B, C = real_posterior(f, 100)
            if C >= A and C >= B:
                F[x][y] = f
                break
                
                
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
x = y = range(100)
X, Y = np.meshgrid(x, y)
F = np.array(F)

ax.plot_surface(X, Y, F)

ax.set_xlabel('Box defects')
ax.set_ylabel('Prior')
ax.set_zlabel('Ft threshold')

plt.show()