Skip to content

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
You can’t perform that action at this time.
You signed in with another tab or window. Reload to refresh your session. You signed out in another tab or window. Reload to refresh your session.