# OK to post homework
# Robert Dougherty-Bliss
# 2021-02-14
# Assignment 9

read "C9.txt":

# These theorems are easy when the mean is 0, and f(x) / x^mean is the pgf of X
# - mean, which has mean 0.

SAc := proc(f, x, K)
    mu := subs(x = 1, diff(f, x)):
    F := f / x^mu:
    F := x * diff(x * diff(F, x), x):
    var := subs(x = 1, F):
    res := [mu, var]:

    for k from 3 to K do
        F := x * diff(F, x):
        res := [op(res), subs(x = 1, F) / var^(k / 2)]:
    od:

    res:
end:

X := RRV(100, 1000);
f := Pgf(X, x);
evalf(SAc(f, x, 8) = SA(X, 8));
evalb(evalf(SAc(f, x, 8) = SA(X, 8)));

# 2.

f := ((1 + x) / 2)^n;

moments := SAc(f, x, 15);
print(simplify(moments));

for k from 3 to 15 do
    print(limit(moments[k], n=infinity));
od:

# Assuming I haven't done anything wrong, the OEIS tells me that the nonzero
# moments scaled moments are probably double factorials of odd integers. This
# is probably just an annoying thing to verify.

# In fact, let's just look at the implied sum:

#   E[(X - n / 2)^k] = sum((j - n / 2)^k * binomial(n, j), j).

# This sum is symmetric in j, i.e., every term j < n / 2 can be paired up with
# its negative, a term j > n / 2. Thus when k is odd, these terms all cancel.
# The j = n / 2 term obviously vanishes, so this explains why the odd moments
# are 0.

# 3.

read "C7.txt":

RandSPn := proc(n)
    # What is the probability that a random set partition has exactly k parts?
    K := LD([seq(NuSPnk(n, k), k=1..n)]):

    RandSPnk(n, K):
end:

with(combinat):

meanTest := proc(n, K)
    local x, y, j, k:
    x := add(nops(RandSPn(n)), k=1..K) / K:
    y := add(j * NuSPnk(n, j), j=1..n) / bell(n):

    evalf([x, y, x - y]);
end:

print(partsTest(100, 1000));
print(partsTest(200, 1000));
print(partsTest(300, 1000));