#C7.txt, Feb. 12, 2015
#Confidence Intervals, Bayesian method
Help:=proc(): print(` N(x), ENHn(p,n,k,eps), HMsd(eps) ,CoInt(p,n,eps), BP(p,n,k)  `): end:

#N(x): The cumulative distribution function of the standard Normal Distribution N(0,1)
N:=proc(x) local t:
 int(exp(-t^2/2)/sqrt(2*Pi),t=-infinity..x):
end:

#HMsd(eps): How many standard-deviations to tolerate with confidence level eps?
HMsd:=proc(eps) local x:

fsolve(N(-x)=eps/2,x):

end:



#ENHn(p,n,k,eps): If someone claims that the Pr(H) is p, and you toss
#it n times, and get k Heads should you believe it or not with CONFIDENCE
#1-eps
ENHn:=proc(p,n,k,eps) local Y,t:

Y:=abs((k-n*p)/sqrt(n*p*(1-p))):

evalb(evalf(int(exp(-t^2/2)/sqrt(2*Pi),t=-infinity..-Y))>=eps/2):
end:

#CoInt(p,n,eps): The confidence interval of level 1-eps
CoInt:=proc(p,n,eps) local s,sd:
sd:=sqrt(n*p*(1-p)):

s:=HMsd(eps):

evalf([n*p-s*sd, n*p+s*sd]):

end:

#BP(p,n,k): The Bayesian "degree of belief" that the coin has probability p
#of Heads if in an experiment with n tosses you got k Heads
BP:=proc(p,n,k) local P,H,C:
#With prob. 1/2 Pr(H)=p, but the remaining probabilities are equally likely
#If P=p then binomial(n,k)*p^k*(1-p)^(n-k) (prob. 1/2)
#If all equally likely
#Total prob. that you got k Heads =binomial(n,k)*p^k*(1-p)^(n-k)*1/2
#1/2*int(binomial(n,k)*P^k*(1-P)^(n-k),P=0..1):
H:=binomial(n,k)*p^k*(1-p)^(n-k)*1/2:
C:=1/2*int(binomial(n,k)*P^k*(1-P)^(n-k),P=0..1):
H/(H+C):

end: