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

read "C7.txt":

# 3.

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:

# 4.

NuAdjacentMates := proc(SP, n)
    parts := table():

    for k from 1 to n do
        for part from 1 to nops(SP) do
            if k in SP[part] then
                parts[k] := part:
                break:
            fi:
        od:
    od:

    add(ifelse(parts[i] = parts[i + 1], 1, 0), i=1..n-1):
end:

EstimateAvAdjacentMates := proc(n, K)
    add(NuAdjacentMates(RandSPn(n), n), j=1..K) / K:
end:

EstimateAvAdjacentMates(10, 1000);

# 1: 0
# 2: 0.5
# 3: 0.8
# 4: 1
# 5: 1.1....

# I don't see any *obvious* pattern yet, just by playing.

# (iii).

# Claim: The expectation is (n - 1) * bell(n - 1) / bell(n).

# expectation = sum(E[j is adjacent to j + 1], j).
# The number of partitions of [n] that j is adjacent to j + 1 is exactly bell(n
# - 1) (remove j, partition the remaining elements, put j into the part where j
# + 1 is). Thus, these smaller expectations are all bell(n - 1) / bell(n),
# except for j = n, which is 0. The total expectation is therefore

#   (n - 1) bell(n - 1) / bell(n).

with(combinat):

guess := n -> (n - 1) * bell(n - 1) / bell(n):

seq(evalf(guess(n)), n=1..10);

# Answers:

# 0., 0.5000000000, 0.8000000000, 1., 1.153846154, 1.280788177, 1.388825542, 1.482850242, 1.566179600,

# Looks right to me! Wow!

# Could we have guessed this? I don't see an *obvious* way, but I would like to
# know one.