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

#Read packages
with(`combinat`):

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

#Help procedure
Help:=proc():
    print(` AcceptNHnormal(p, n, k, eps) , Phi(x) `):
    print(` Problem2a(n) , Problem2b(n) , Problem3(p, L, n) `):
    print(` Pr(P, x, n, a, b) , CheckCLT(P, x, n, c1, c2) `):
end:

##PROBLEM 1##

#AcceptNHnormal(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
AcceptNHnormal:=proc(p, n, k, eps) local Y, t:
    Y:=abs((k-n*p)/sqrt(n*p*(1-p))):
    return evalb(evalf(int(exp(-t^2/2)/sqrt(2*Pi), t=-infinity..-Y)) >= eps(2)):
end:

#The binomial distribution is asymptotically normal

##PROBLEM 2##

Problem2a:=proc(n) local k:
    return add(binomial(n,k)*binomial(n+k,k)*x^k,k=0..n):
end:

#evalf(ID(Problem2a(40), x, 8)); returns
#[28.26254448, 7.097306998, -.9468059113e-1, 2.956207200,
#-.9326335538, 14.44204252, -9.504014540, 98.56104812]
#evalf(ID(Problem2a(60), x, 8)); returns
#[42.40477438, 10.63272270, -.7712402503e-1, 2.970714819,
#-.7635338678, 14.62435663, -7.859059156, 100.6272295]
#These are close to the moments of the standard normal distribution

Problem2b:=proc(n) local k:
    return add(binomial(n,k)^3*x^k,k=0..n):
end:

#evalf(ID(Problem2b(40), x, 8)); returns
#[20., 3.389505939, 0., 2.983618787, 0., 14.75535645, 0.,
#101.5992776]
#evalf(ID(Problem2b(60), x, 8)); returns
#[30., 5.055967009, 0., 2.989014625, 0., 14.83570234, 0.,
#102.7107728]
#These are close to the moments of the standard normal distribution

##PROBLEM 3##

#Probability of rolling L[i] on a die is p[i]
#Check up to kth moment
Problem3:=proc(p, L, k) local P, x, M, i:
    P:=add(p[i]*x^(L[i]), i=1..nops(L)):
    M:=ID(P^n, x, k):
    return [seq(limit(M[i], n=infinity), i=1..nops(M))]:
end:

#Everything I tried returned [infinity, infinity, 0, 3, 0, 15],
#as expected for an asymptotically normal random variable

##PROBLEM 5##

#Pr(P, x, n, a, b): computes the probability that if the die whose
#probability generating function (for one roll) is the polynomial
#P(x) in x, is rolled n times, the total sum of the dots is
#between a and b (inclusive)
Pr:=proc(P, x, n, a, b) local Pton, i:
    Pton:=expand(P^n):
    return add(coeff(Pton, x, i), i=a..b):
end:

##PROBLEM 6 and PROBLEM 4##

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

#E is n*P'(1)
#V is n*(x*(x^-P'(1)*P(x))')'(1)
CheckCLT:=proc(P, x, n, c1, c2) local E, V:
    E:=subs(x=1, diff(P, x)):
    V:=n*subs(x=1, diff(x*diff(x^(-E)*P, x), x)):
    E:=n*E:
    return evalf([Pr(P, x, n, ceil(E+c1*sqrt(V)), floor(E+c2*sqrt(V))), Phi(c2) - Phi(c1)]):
end:

#CheckCLT((x+2*x^2+x^3)/4, x, 200, -2, 2); returned
#[0.9597692603, 0.9544997356].  These numbers are close to each other