Instantly share code, notes, and snippets.

Embed
What would you like to do?
# -*- 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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment