#OK to post homework
#Karnaa Mistry, 10/18/20, Assignment HW #12

with(combinat):

# 1. 

#CPg(L): Given list of lists of 'databases', outputs CP(L1,L2,...,Lk)
CPg:=proc(L) local i,Lx,Ly:

Lx:=L[1]:
for i from 2 to nops(L) do
  Ly:=CP(Lx,L[i]):
  Lx:=Ly:
od:

RETURN(Lx):

end:




# 2.

#AveGF(f,x): AveClever, though with already known weight-enumerator
AveGF:=proc(f,x):
subs(x=1,diff(f,x))/subs(x=1,f):
end:


#kthMomentGF(f,x,k): kthMomentClever, but with known weight-enumerator
kthMomentGF:=proc(f,x,k) local mu,f1,i,c:
mu:=AveGF(f,x):
f1:=f/x^mu:
c:=0:

for i from 0 to degree(f,x) do
  c:=c+coeff(f,x,i):
od:

for i from 1 to k do
  f1:=expand(x*diff(f1,x)):
od:

subs(x=1,f1)/c:
end:




# 3.

#ScaledMomentGF(f,x,k): inputs a given weight-enumerator, and finds the scaled moment about the mean
ScaledMomentGF:=proc(f,x,k) local mk,m2,c:

mk:=kthMomentGF(f,x,k):
m2:=kthMomentGF(f,x,2):

if m2 = 0 then RETURN("Error: 2nd moment is 0"): fi:

c:=mk/(m2)^(k/2):

RETURN(c):

end:




# 4.

# First, let us note that odd k values give us 0 from ScaledMomentGF.  So, let's start inspection at k=2.
# For k=2, ScaledMomentGF((1+x)/2)^n,x,2) with n=1..10 gives us all 1's.  So, f_2 = 1.
# For k=4, with n=1..5 we have 1, 4/2, 7/3, 10/4, 13/5, which gives us a change of 1 in the denominators,
#    and a change of 3 in the numerators, and so we have f_4 = (3n-2)/n
# For k=6, it gets harder, and we have 1, 16/4, 61/9, 136/8...  The denominator is just n^2, though the
#    numerator is much trickier.  1,16,61,136 is not even in the OEIS.  I wasn't sure how to obtain an
#    explicit expression, but did some research and found a package in Maple that lets you input data
#    points, and it uses Polynomial Interpolation to find an appropriate function.  Since the first changes
#    for the numerator are 15,45,75,etc. and the second changes are equally 30, we know that this will be
#    some type of polynomial.  So, using with(Student[NumericalAnalysis]): we have something like
#    expand(PolynomialInterpolation(xy,independentvar=x,method=newton)); to give us the numerator of
#    15n^2 - 30n + 16.  So, f_6 = (15n^2 - 30n + 16) / n^2
# For k=8, we do something similar, and see the denominators are n^3, and use Maple to find the numerator
#    function, 105n^3 - 420n^2 + 588n - 272.  So, f_8 = (105n^3 - 420n^2 + 588n - 272) / n^3
# For k=10, I was unable to obtain a consistent function using the above methods.

# Calculating limits, we get a sequence like (1,0,3,0,15,0,105 ...), and so this is in the OEIS, with
# number A123023, "a(n)=(n-1)*a(n-2),a(0)=1,a(1)=0"
# a(n) is also commented as the n-th moment of the standard normal distribution