# OK to post homework
# Robert Dougherty-Bliss
# 2021-02-04
# Assignment 4

read "C4.txt":

# 1. NuSnk(n, k) is some double sequence determined by the usual Pascal
# recurrence, which happens to satisfy the usual Pascal initial conditions.
# Fortunately for us, it is routine to verify that binomial(n, k) *also*
# satisfies these conditions!

simplify(convert(binomial(n, k) - binomial(n - 1, k) - binomial(n - 1, k - 1), factorial));

# 2.

# LD only cares about the proportion of the quantities, not the quantites
# themselves.

# In this case, binomial(n - 1, k) + binomial(n - 1, k - 1) = binomial(n, k),
# and binomial(n - 1, k) / binomial(n, k) = (n - k) / n, and
# binomial(n - 1, k - 1) / binomial(n, k) = k / n.
# So, with probability (n - k) / n we pick 1, and k / n we pick 2. But this is
# exactly what LD([n - k, k]) does.

# Obviously this is going to be faster becauase LD([A, B]) actually constructs
# a list of length A + B, and n <= binomial(n, k) for k >= 1.

RnkFast:=proc(n,k)  local c,S:
if k<0 or k>n then
 RETURN(FAIL):
fi:

if n=0 then
  if k=0 then
    RETURN({}):
   fi:
fi:

c:=LD([n - k, k]):

if c=1 then
  RETURN(Rnk(n-1,k)):
else
 RETURN(Rnk(n-1,k-1) union {n}):
fi:

end:

# 3.

EstimateAveMax := proc(n, k, K)
    evalf(add(max(RnkFast(n, k)), j=1..K) / K):
end:

# Claim: E[max(Rnk(n, k))] = k * (n + 1) / (k + 1).

# How many subsets have j as a maximum? Exactly binomial(j - 1, k - 1), the
# ones that use j and then k - 1 elements of {1, 2, ..., j - 1}. There are
# binomial(n, k) subsets, so the probability of having j as a maximum is
# binomial(j - 1, k - 1) / binomial(n, k).

# This means that the expectation is

#   sum(j * binomial(j - 1, k - 1), 1 <= j <= n) / binomial(n, k).

# In principle, this sum is utterly trivial. In practice, it's fun to do by
# hand. Write j * binomial(j - 1, k - 1) = k * binomial(j, k). (This is true
# for k >= 1, at least.) Then the sum is

#   k * sum(binomial(j, k), 1 <= j <= n)

# which is well-known via the discrete calculus to evaluate as

#       k (binomial(n + 1, k + 1) - binomial(1, k + 1))
#    =  k * binomial(n + 1, k + 1).

# So the full expectation is k * binomial(n + 1, k + 1) / binomial(n, k), which
# is k * (n + 1) / (k + 1).

# 4.

EstimateAveSum := proc(n, k, K)
    evalf(add(add(RnkFast(n, k)), j=1..K) / K):
end:

# Claim: E[sum(Rnk(n, k))] = k (n + 1) / 2.

# This expectation is exactly

#       sum(sum(S), S <= [n], |S| = k) / binomial(n, k).

# Every number in [n] shows up in exactly binomial(n - 1, k - 1) subsets, so
# the double sum reduces to

#   sum(binomial(n - 1, k - 1) * j, 1 <= j <= n) / binomial(n, k),

# and this is routine to evaluate and simplify.

# 5.

WordRecurrence := proc(a, b, initA, initB, fA, fB, comb)
    WR := (a, b) -> WordRecurrence(a, b, initA, initB, fA, fB, comb):

    if a < 0 or b < 0 then
        return FAIL:
    fi:

    if a = 0 then
        return initB(b)
    elif b = 0 then
        return initA(a)
    fi:

    start_1 := fA(WR(a - 1, b)):
    start_2 := fB(WR(a, b - 1)):

    comb(start_1, start_2):
end:

Wnk := proc(a, b)
    option remember:
    initA := a -> [[1 $ a]]:
    initB := b -> [[2 $ b]]:
    fA := ws -> [seq([1, op(w)], w in ws)]:
    fB := ws -> [seq([2, op(w)], w in ws)]:
    comb := (x, y) -> [op(x), op(y)]:
    WordRecurrence(a, b, initA, initB, fA, fB, comb):
end:

NuWnk := proc(a, b)
    option remember:
    initA := a -> 1:
    initB := b -> 1:
    fA := ws -> ws:
    fB := ws -> ws:
    comb := (x, y) -> x + y:
    WordRecurrence(a, b, initA, initB, fA, fB, comb):
end:

# What's the probability that a word chosen uniformly from Wnk(a, b) at random
# ends in a 1? Exactly NuWnk(a - 1, b) / NuWnk(a, b). So, with probability
# equal to that, end in a 1, and otherwise end in a 2.

# And we "know" that NuWnk(a, b) = binomial(a + b, a), so the probability is
# exactly a / (a + b). (Which feels obvious, I guess.)

RWnk := proc(a, b)
    if a = 0 then
        return [2 $ b]:
    elif b = 0 then
        return [1 $ a]
    fi:

    if rand(0..1.0)() < a / (a + b) then
        [1, op(RWnk(a - 1, b))]:
    else
        [2, op(RWnk(a, b - 1))]:
    fi:
end:

# Let me be annoying and prove that the probability that any particular word is
# chosen is indeed 1 / binomial(a + b, a). (Because, however intuitive it seems
# to be, I don't see how the proof is obvious!)

# Suppose that we can pick words of length < n uniformly at random. Grab any
# word of length n - what is the probability that we chose it?

# If it ends in 1, then it is

#       P(chose 1 with a, b) * P(chose word of length n - 1 with a - 1, b)
#    =  a / (a + b) / binomial(a - 1 + b, a - 1)
#    =  1 / binomial(a + b, a).

# A similar argument works if it ends in 2.

# 6 / 7. (I'm just going to to 7, which covers 6.)

# We could do a similar "generic" function for flexibility, but I'd rather be
# done in time for the Super Bowl tonight :)

WORDS := proc(L)
    if add(L) = 0 then
        return [[]]:
    end:

    res := []:

    for k from 1 to nops(L) do
        if L[k] > 0 then
            newL := L:
            newL[k] := newL[k] - 1:
            newWords := [seq([k, op(w)], w in WORDS(newL))]:
            res := [op(res), op(newWords)]:
        fi:
    od:

    res:
end:

NuWORDS := proc(L)
    option remember:
    if add(L) = 0 then
        return 1:
    end:

    total := 0:

    for k from 1 to nops(L) do
        if L[k] > 0 then
            newL := L:
            newL[k] := newL[k] - 1:
            total := total + NuWORDS(newL):
        fi:
    od:

    total:
end:

RandWORD := proc(L)
    if add(L) = 0 then
        return []
    fi:

    r := rand(0..1.0)():

    for k from 1 to nops(L) do
        if L[k] = 0 then
            next:
        fi:

        newL := L:
        newL[k] := newL[k] - 1:
        p := NuWORDS(newL) / NuWORDS(L):
        if r <= p then
            return [k, op(RandWORD(newL))]:
        fi:

        r := r - p:
    end:

    # We shouldn't ever get here!
    return FAIL:
end: