#h6.txt:Binomial distributions
Help:=proc(): print(`  HT(n,p,x), CheckHT(n,p,x,N):`): end:


with(Statistics):

#Bp(n,p,k): The prob. mass function of the Binomial distribution, tossing a coin n times
#Pr(H)=p, k=0..n, Pr(X=k)
Bp:=proc(n,p,k)  binomial(n,k)*p^k*(1-p)^(n-k): end:


#HT(n,p,x): Frequentist Hypothesis testing. If the owner of the coin claims that Pr(H) is at least
#p what is the cutoff of number of Heads below which you can be confident with level of confidence
#x that hs is a lier.
HT:=proc(n,p,x) local k1,k:

for k from 0 to n while add(Bp(n,p,k1),k1=0..k)<=1-x do od:
k-1:
end:




#CheckHT(n,p,x,N):
CheckHT:=proc(n,p,x,N) local k,co,i:

k:=HT(n,p,x):

co:=0:

for i from 1 to N do
 if convert(Sample(RandomVariable(BernoulliDistribution(p)),n), `+` )<=k then
  co:=co+1:
 fi:
od:

evalf(co/N):

end: