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

#Read packages
with(`combinat`):

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

#Help procedure
Help:=proc():
    print(` PGL(p, N, x) , IDlimit(f, x, k) `):
end:

##PROBLEM 2##

#PGL(p, N, x): inputs a numeric or symbolic probability p, a
#numeric positive integer N, and a variable x, and outputs the
#list of length N-1 whose n-th entry is the probability generating
#function for the random variable: "duration of game" started with
#n dollars, ending with 0 or N, and with probability of a win in
#each round p.
PGL:=proc(p, N, x) local E, eq, var, i, sol:
    var:={seq(E[i],i=0..N)}:
    eq:={E[0]=x, E[N]=x, seq(E[i]=x*(p*E[i+1]+(1-p)*E[i-1]), i=1..N-1)}:
    sol:=solve(eq, var):
    return subs(sol, [seq(E[i],i=1..N-1)]):
end:

##PROBLEM 3##

#IDlimit(f, x, k): Given the probability generating function f
#in the variable x, for some r.v. under some probability
#distribution outputs the LIST [Exp, Var, std. moms]
#NOTE: the Exp is the true expectation, not automatically zero
#Never evaluates at x=1, just takes limits
IDlimit:=proc(f, x, k) local i, g, av, M:
    g:=f/limit(f, x=1):
    av:=limit(diff(g, x), x=1):
    M:=[av]:
    g:=g/x^av:
    g:=x*diff(g, x):
    for i from 2 to k do
        g:=x*diff(g, x):
        M:=[op(M), limit(g, x=1)]:
    od:
    return [M[1], M[2], seq(M[i]/M[2]^(i/2), i=3..k)]:
end:

##PROBLEM 4##

#I ran
#IDlimit(PGFL(1/2,n,2*n,x),x,4);
#It outputs
#[n^2, 2/3*n^4-2/3*n^2, (4/15*n^2-4/3*n^4+16/15*n^6)/(2/3*n^4-2/3*n^2)^(3/2), (-104/15*n^6+412/105*n^8+8/105*n^2+44/15*n^4)/(2/3*n^4-2/3*n^2)^2]
#These are the expectation, variance, and standardized third and
#fourth moments
#The third moment goes to 4*sqrt(6)/5 as n goes to infinity
#The fourth moment goes to 309/35 as n goes to infinity
#Conclusion: Not asymptotically normal