%hide
%auto
##
## Illustration of Bayes Theorem for introductory Probability and Statistics
##
## John Travis
## Mississippi College
## 06/17/11
##
## The user inputs:
## - A list of probabilities corresponding to each sector of a partition
## - A list of prior conditional probabilities P( A | B_1), etc.
## Output:
## - Two Venn diagrams. The first describes the input data. The second describes the posterior probabilities.
## - Numerical values for everything, if desired.
##
## An updated version of this worksheet may be available at http://sagenb.mc.edu
##
# This function is used to convert an input string into separate entries
def g(s): return str(s).replace(',',' ').replace('(',' ').replace(')',' ').split()
@interact
def _(Partition_Probabilities=input_box('0.35,0.25,0.40',label="$P(B_1),P(B_2),...$"),
Conditional_Probabilities=input_box('0.02,0.01,0.03',label='$P(A|B_1),P(A|B_2),...$'),
print_numbers=checkbox(True,label='Numerical Results on Graphs?'),
auto_update=False):
Partition_Probabilities = g(Partition_Probabilities)
Conditional_Probabilities = g(Conditional_Probabilities)
n = len(Partition_Probabilities)
n0 = len(Conditional_Probabilities)
if n<>n0:
print html("<font size=+2 color=red>You must have the same number of partition probabilities and conditional probabilities.</font>")
else: # input data streams now are the same size!
colors = rainbow(n)
accum = float(0) # to test whether partition probs sum to one
ends = [0] # where the graphed partition sectors change in pie chart
mid = [] # middle of each pie chart sector used for placement of text
p_Bk_given_A = [] # P( B_k | A )
pA = 0 # P(A)
PP=[] # array to hold the numerical Partition Probabilities
CP=[] # array to hold the numerical Conditional Probabilities
for k in range(n):
PP.append(float(Partition_Probabilities[k]))
CP.append(float(Conditional_Probabilities[k]))
p_Bk_given_A.append(PP[k]*CP[k] )
pA += p_Bk_given_A[k]
accum = accum + PP[k]
ends.append(accum)
mid.append((ends[k]+accum)/2)
#
# Marching along from 0 to 1, saving angles for each partition sector boundary.
# Later, we will multiple these by 2*pi to get actual sector boundary angles.
#
if abs(accum-float(1))>0.0000001: # Due to roundoff issues, this should be close enough.
print html("<font size=+1 color=red>Sum of probabilities should equal 1.</font>")
else: # probability data is sensible
#
# Draw the Venn diagram by drawing sectors from the angles determined above
# First, create a circle of radius 1 to illustrate the the sample space S
# Then draw each sector with varying colors and print out their names on the edge
#
G = circle((0,0), 1, rgbcolor='black',fill=False, alpha=0.4,aspect_ratio=True,axes=False,thickness=5)
for k in range(n):
G += disk((0,0), 1, (ends[k]*2*pi, ends[k+1]*2*pi), color=colors[mod(k,10)],alpha = 0.2)
G += text('$B_'+str(k+1)+'$',(1.1*cos(mid[k]*2*pi), 1.1*sin(mid[k]*2*pi)), rgbcolor='black')
G += circle((0,0), 0.6, facecolor='yellow', fill = True, alpha = 0.1, thickness=5,edgecolor='black')
# Print the probabilities corresponding to each particular region as a list and on the graphs
if print_numbers:
html("$P(A) = %s$"%(str(pA),))
for k in range(n):
html("$P(B_{%s} | A)$"%(str(k+1))+"$ = %s$"%str(p_Bk_given_A[k]/pA))
G += text(str(p_Bk_given_A[k]),(0.4*cos(mid[k]*2*pi), 0.4*sin(mid[k]*2*pi)), rgbcolor='black')
G += text(str(PP[k] - p_Bk_given_A[k]),(0.8*cos(mid[k]*2*pi), 0.8*sin(mid[k]*2*pi)), rgbcolor='black')
# This is essentially a repeat of some of the above code but focused only on creating the smaller inner circle dealing
# with the set A so that the sectors now correspond in area to the Bayes Theorem probabilities
accum = float(0)
ends = [0] # where the graphed partition sectors change in pie chart
mid = [] # middle of each pie chart sector used for placement of text
for k in range(n):
accum += float(p_Bk_given_A[k]/pA)
ends.append(accum)
mid.append((ends[k]+accum)/2)
H = circle((0,0), 1, rgbcolor='black',fill=False, alpha=0,aspect_ratio=True,axes=False,thickness=0)
H += circle((0,0), 0.6, facecolor='yellow',fill=True, alpha=0.1,aspect_ratio=True,axes=False,thickness=5,edgecolor='black')
for k in range(n):
H += disk((0,0), 0.6, (ends[k]*2*pi, ends[k+1]*2*pi), color=colors[mod(k,10)],alpha = 0.2)
H += text('$B_'+str(k+1)+'|A$',(0.7*cos(mid[k]*2*pi), 0.7*sin(mid[k]*2*pi)), rgbcolor='black')
# Now, print out the bayesian probabilities using the smaller set A only
if print_numbers:
for k in range(n):
H += text(str( N(p_Bk_given_A[k]/pA,digits=4) ),(0.4*cos(mid[k]*2*pi), 0.4*sin(mid[k]*2*pi)), rgbcolor='black')
html.table([[G,H],['Venn diagram of partition with A in middle','Venn diagram presuming A has occured']])
|
Click to the left again to hide and once more to show the dynamic interactive window
|