#Aniket Shah
#2/14/15
#Ok to Post

###########C7.txt

Help:=proc(): print(`N(x), ANHn(p,n,k,eps), HMsd(eps), CoInt(p,n,eps), PB(p,n,k)`): end:

N:=proc(x) local t:
  int(exp(-t^2/2)/sqrt(2*Pi),t=-infinity..x):
end:

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:

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

evalb(evalf(int(exp(-t^2/2)/sqrt(2*Pi),t=-infinity..-Y))>=evalf(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:

#The Bayesian "degree of belief" that the coin has probability p
#of heads if, in an experiment with n tosses, you got k heads
PB:=proc(p,n,k) local P,H,C:
#with prob. 1/2, Pr(H)=p, but the remaining prbabilities are equally likely
#If P=p then binomial(n,k)*p^k*(1-p)^(n-k) (prob. 1/2)
#If P<>p, 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:

############end of C7

#1
PBg:=proc(p,n,k,c0,a,b) local P,H,C:
  H:=binomial(n,k)*p^k*(1-p)^(n-k)*(c0):
  C:=(1-c0)*int(binomial(n,k)*P^k*(1-P)^(n-k),P=0..1):
  H/(H+C):
end:

#2
OneRun:=proc(p,n) local numheads,i:

  numheads:=0:
  for i from 1 to n do
    if rand(1..denom(p))() <= numer(p) then
      numheads:=numheads+1:
    fi:
  od:

numheads:
end:

#3
ENHnS:=proc(p,n,eps,K) local numaccepts, i:

  numaccepts:=0:
  for i from 1 to K do
    if ENHn(p,n,OneRun(p,n),eps)=true then
      numaccepts:=numaccepts+1:
    fi:
  od:
 
evalf(numaccepts/K):
end:

#5
SimulatePoisson:=proc(n,a,K) local i, k,L:

  L:=[0$(n+1)]:

  for i from 1 to K do
    k:=(OneRun(a/n,n)):
    L[k+1]:=L[k+1]+1:

  od:

evalf(L/K):
end: