#Nathan Fox
#Homework 7
#I give permission for this work to be posted online

#Read packages
with(`combinat`):

#Read procedures from class
read(`C7.txt`):

#Help procedure
Help:=proc():
    print(` BPg(p, n, k, c0, a, b) , OneRun(p, n) `):
    print(` ENHnS(p, n, eps, K) , SimulatePoisson(n, a, K) `):
end:

##PROBLEM 1##

#BPg(p, n, k, c0, a, b): c0 is the prior probability that the
#creator of the coin actually made the probability of a Head p,
#and if not, with prob. 1-c0, he chooses, uniformaly at random,
#a probability from the interval [a,b].
#For example BPg(p, n, k, 1/2, 0, 1) should give the same output as
#BP(p, n, k) in file C7.txt
BPg:=proc(p, n, k, c0, a, b) local P, H, C:
    #Honest
    H:=c0*binomial(n, k)*p^k*(1-p)^(n-k):
    #Cheat
    C:=(1-c0)*int(binomial(n, k)*P^k*(1-P)^(n-k), P=a..b):
    #Posterior probability
    return H/(H+C):
end:

##PROBLEM 2##

#OneRun(p, n): simulates tossing, n times, a loaded coin whose
#probability of Heads in p (assume that p is a rational number a/b,
#so that you can use the rand(1..b) command to simulate a coin
#toss with Prob. a/b), and returns the number of times it got Heads.
OneRun:=proc(p, n) local a, b, r, toss, i:
    a:=numer(p):
    b:=denom(p):
    r:=rand(1..b):
    toss:=proc()
        if r() <= a then
            return 1:
        else
            return 0:
        fi:
    end:
    return add(toss(), i=1..n):
end:

##PROBLEM 3##

#ENHnS(p, n, eps, K): runs OneRun(p, n) K times,
#looks at its output, and, using ENHn(p, n, k, eps), decides
#whether to accept the Null Hypothesis, and record the ratio of
#those where you accepted it (correctly, of course).
ENHnS:=proc(p, n, eps, K) local count, i:
    count:=0:
    for i from 1 to K do
        if ENHn(p, n, OneRun(p, n), eps) then
            count:=count+1:
        fi:
    od:
    return count / K:
end:

#The values are consistently near, but slightly below, 0.95
#for ENHnS(1/3,300,1/20,1000).  Values around 0.91 are more common

##PROBLEM 5##

#SimulatePoisson(n, a, K): n is a positive integer, K is a positive
#integer, and a is a small integer
#It runs the procedure OneRun(a/n, n) K times and outputs the list
#L=[L0,L1,L2,L3], where L0 is the ratio (out of the K trials) that
#yielded 0, L1 is the ratio (out of the K trials) that yielded 1,
#etc.
SimulatePoisson:=proc(n, a, K) local L, i, trial:
    L:=[0$(n+1)]:
    for i from 1 to K do
        trial:=OneRun(a/n, n):
        L[trial+1]:=L[trial+1] + 1:
    od:
    return [seq(L[i]/K, i=1..nops(L))]:
end:

#The results of this procedure agree with the theoretical values