%matplotlib inline

px = linspace(0,1,50)
plot( px, px**5)
xlabel(r'$\theta$',fontsize=18)
ylabel(r'$\beta$',fontsize=18)

from sympy.abc import p,k # get some variable symbols
import sympy as S

expr=S.Sum(S.binomial(5,k)*p**(k)*(1-p)**(5-k),(k,3,5)).doit()
p0=S.plot(p**5,(p,0,1),xlabel='p',show=False,line_color='b')
p1=S.plot(expr,(p,0,1),xlabel='p',show=False,line_color='r',legend=True,ylim=(0,1.5))
p1.append(p0[0])
p1.show()

# implement 1st test where all five must be heads to declare underlying parameter > 0.5
from scipy import stats
b = stats.bernoulli(0.8) # true parameter is 0.8 (i.e. hypothesis H_1)
samples = b.rvs(1000).reshape(-1,5) # -1 means let numpy figure out the other dimension i.e. (200,5)
print 'prob of detection = %0.3f'%mean(samples.sum(axis=1)==5) # approx 0.8**3

# here's the false alarm case
b = stats.bernoulli(0.3) # true parameter is 0.3 (i.e. hypothesis H_0)
samples = b.rvs(1000).reshape(-1,5)
print 'prob of false alarm = %0.3f'%mean(samples.sum(axis=1)==5)

# implement majority vote test where three of five must be heads to declare underlying parameter > 0.5
b = stats.bernoulli(0.8) # true parameter is 0.8 (i.e. hypothesis H_1)
samples = b.rvs(1000).reshape(-1,5)
print 'prob of detection = %0.3f'%mean(samples.sum(axis=1)>=3) 

# here's the false alarm case
b = stats.bernoulli(0.3) # true parameter is 0.3 which means it's hypothesis H_0
samples = b.rvs(1000).reshape(-1,5)
print 'prob of false alarm = %0.3f'%mean(samples.sum(axis=1)>=3)