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

#Read packages
with(`combinat`):

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

#Help procedure
Help:=proc():
    print(` Alphas(f, x, k) , Momslc(n, p, k) , SymDie(f, n, k) `):
    print(` CP(P1, L1, P2, L2) , FindDice(minface, maxface) `):
end:

##AUXILIARY PROCEDURES (for multiple problems)##

#Alphas(f, x, k): the list [av, var, standardized k-th moment]
#for p.g.f. f
Alphas:=proc(f, x, k) local M, i:
    M:=MomNs(f, x, k):
    return [M[1], M[2], seq(M[i]/M[2]^(i/2), i=3..k)]:
end:

##PROBLEM 1##

#Rigorously derives expressions in n and p for the expected value,
#variance, and the third through sixth standardized moments for the
#random variable "total amount gained" after n tosses of a loaded
#coin with the probability of Heads is p and the probability of
#Tails is 1-p, and you win a dollar if you get Heads and lose a
#dollar if you get Tails.
Momslc:=proc(n, p, k) local x:
    return Alphas((p*x+(1-p)/x)^n, x, k):
end:
#simplify(Momslc(6)) returns
#[n*(2*p-1),
#4*n*p-4*n*p^2,
#-(2*p-1)/(-n*p*(p-1))^(1/2),
#(3*n*p^2-6*p^2-3*n*p+6*p-1)/n/p/(p-1),
#-(-24*p^3+20*n*p^3+36*p^2-30*n*p^2-14*p+10*n*p+1)/(-n*p*(p-1))^(1/2)/n/p/(p-1),
#(120*p^4-130*n*p^4+15*n^2*p^4+260*n*p^3-30*n^2*p^3-240*p^3-155*n*p^2+150*p^2+15*n^2*p^2+25*n*p-30*p+1)/n^2/p^2/(p-1)^2]

#L:=simplify(Momslc(6))
#limit(L[3], n=infinity) returns 0
#limit(L[4], n=infinity) returns 3
#limit(L[5], n=infinity) returns 0
#limit(L[6], n=infinity) returns 15

##PROBLEM 2##

#SymDie(f, n, k): inputs a specific positive integer, f, a symbol
#(or number) n, and a positive integer k larger than 2, and outputs
#the expectation, variance, and third through k-th moment of the
#random variable "total number of dots" if you roll a FAIR f-faced
#die n times.
#if standardized is true, return standardized moments
#otherwise, return straight moments
SymDie:=proc(f, n, k, standardized:=false) local momfun, x, i:
    if standardized then
        momfun:=Alphas:
    else
        momfun:=MomNs:
    fi:
    return momfun(add(x^i, i=1..f)^n, x, k):
end:

#seq(SymDie(f, n, 2), f=1..10); outputs
#[n, 0], [3/2*n, 1/4*n], [2*n, 2/3*n], [5/2*n, 5/4*n], [3*n, 2*n],
#[7/2*n, 35/12*n], [4*n, 4*n], [9/2*n, 21/4*n], [5*n, 20/3*n],
#[11/2*n, 33/4*n]
#Conjecture: The expectation is (f+1)*n/2
#Conjecture: The variance is (f-1)*(f+1)*n/12 (obtained by observing
#denominators of 2*f and then OEIS-ing the numerators and obtaining
#the tetrahedral numbers)
#[seq([seq(limit([seq(SymDie(f, n, 10, true), f=1..10)][i][j], n=infinity), j=3..10)], i=1..nops(L))]; outputs
#[[0, 3, 0, 15, 0, 105, 0, 945], [0, 3, 0, 15, 0, 105, 0, 945],
#[0, 3, 0, 15, 0, 105, 0, 945], [0, 3, 0, 15, 0, 105, 0, 945],
#[0, 3, 0, 15, 0, 105, 0, 945], [0, 3, 0, 15, 0, 105, 0, 945],
#[0, 3, 0, 15, 0, 105, 0, 945], [0, 3, 0, 15, 0, 105, 0, 945],
#[0, 3, 0, 15, 0, 105, 0, 945]]
#This proves that these random variables (at least for f up to 10)
#are asymptotically normal (at least up to the 10th moment)

##PROBLEM 3##

#CP(P1, L1, P2, L2): inputs
#1. a finite probability distribution P1, of length n1, say,
#and a r.v. L1, of length n1,
#2. a finite probability distribution P1, of length n2, say,
#and a r.v. L2, of length n2,
#and outputs the pair [P, L] where P is a list of size n1*n2 where
#the probability distribution defined on pairs [i, j]
#(1 <= i <= n1, 1 ≤ j ≤ n2, ordered in lexicographic order),
#and L is the r.v. defined by
#L[i,j]:=L1[i]+L2[j]
CP:=proc(P1, L1, P2, L2) local i, j:
    return [[seq(seq(P1[i]*P2[j], j=1..nops(P2)), i=1..nops(P1))],
        [seq(seq(L1[i]+L2[j], j=1..nops(L2)), i=1..nops(L1))]]:
end:

#I verified experimentally that
#PGF(op(CP(P1,L1,P2,L2)),x)=PGF(P1,L1,x)*PGF(P2,L2,x)

##PROBLEM 4##

#The left side equals sum(sum(P1[i]*P2[j]*x^(L1[i]+L2[j]), j=1..nops(L2)), i=1..nops(L1))
#The right side equals sum(P1[i]*x^(L1[i]), i=1..nops(L1))*sum(P2[i]*x^(L2[i]), i=1..nops(L2))
#These are equal by the definition of polynomial multiplication

##PROBLEM 5##

#The standard die used in Backgammon and other games of chance has
#6 faces marked with 1,2,3,4,5,6. We are supposed to design a pair
#of two fair dice with faces showing
#a1,a2,a3,a4,a5,a6 (repeats are allowed)
#and
#b1,b2,b3,b4,b5,b6 (repeats are allowed)
#such that if you throw them, you would have the same probability
#of getting the sum to be k as if you threw the usual 2 dice, for
#every k from 2 to 12. In other words such that
#PGF(op(CP([(1/6)$6], [a1,a2,a3,a4,a5,a6], [(1/6)$6], [b1,b2,b3,b4,b5,b6])), x) = PGF(op(CP([(1/6)$6], [1,2,3,4,5,6], [(1/6)$6], [1,2,3,4,5,6])), x)

#I'm not very clever, so let's have Maple find these dice for us
FindDice:=proc(minface, maxface) local a, b, ret, i, j, pg:
    a:=[minface$6]:
    b:=[minface$6]:
    ret:={}:
    pg:=expand(PGF(op(CP([(1/6)$6], [1,2,3,4,5,6], [(1/6)$6], [1,2,3,4,5,6])), x)):
    while true do
        #Check the dice
        if evalb(expand(PGF(op(CP([(1/6)$6], [a[1],a[2],a[3],a[4],a[5],a[6]], [(1/6)$6], [b[1],b[2],b[3],b[4],b[5],b[6]])), x)) = pg) then
            ret:=ret union {[a, b]}:
        fi:
        #Update a
        i:=nops(a):
        #We must have at least one 1 on a, since 2 must be a sum
        while i > 1 do
            if a[i] < maxface then
                a[i]:=a[i] + 1:
                for j from i+1 to nops(a) do
                    a[j]:=a[i]:
                od:
                break:
            else
                a[i]:=minface:
            fi:
            i:=i-1:
        od:
        if i <= 1 then
            #Update b
            i:=nops(b):
            #We must have at least one 1 on b,
            #since 2 must be a sum
            while i > 1 do
                #We don't need both dice to exceed 6, since we
                #can't have the total exceed 12
                if b[i] < 6 then
                    b[i]:=b[i] + 1:
                    for j from i+1 to nops(a) do
                        b[j]:=b[i]:
                    od:
                    break:
                else
                    b[i]:=minface:
                fi:
                i:=i-1:
            od:
            if i <= 1 then
                break:
            fi:
        fi:
    od:
    return ret:
end:

#FindDice(1, 8) returned
#{[[1, 2, 3, 4, 5, 6], [1, 2, 3, 4, 5, 6]],
#[[1, 3, 4, 5, 6, 8], [1, 2, 2, 3, 3, 4]]}
#The first element of this set is the standard dice, which clearly
#have the same probability distribution as standard dice.
#The second element, a=[1,3,4,5,6,8] and b=[1,2,2,3,3,4] are
#another pair of dice with combined probability distribution
#the same as standard dice.

#AFTER I found these dice, I read about them on Wikipedia, which
#made me realize that factoring the generating function would have
#been much easier (though this computation took less than a minute)
#As I said, I'm not very clever.