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

#Read packages
with(`combinat`):

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

#Help procedure
Help:=proc():
    print(` MHg(p) , MHgs(p) , MHgsExperiment(pset, N:=1000) `):
    print(` SDsim(N) , BDg(k, N, r) , BDgsim(k, N, r, M)`):
    print(` BDgf(k, N, r) , IsDerangement(pi) , HowLikelyIsDer(n, K) `):
end:

##PROBLEM 2##

#MHg(p): p is a rational number between 0 and 1 (inclusive),
#and your probability of guessing the correct door
#(by our convention, door #1) is p.
#Simulates ONE instance of the Monty Hall game show
#where there is a car behind one door (#1)
#and the other doors have goats.
#implementing the switching strategy
MHg:=proc(p, verbose:=true) local door1, UnpickedDoors, MontyDoors, MontyDoor, ND, printv, i, roll:
    printv:=proc(str)
        if verbose then
            print(op(str)):
        fi:
    end:
    roll:=rand(1..denom(p))():
    if roll <= numer(p) then
        door1:=1:
    else
        door1:=rand(2..3)():
    fi:
    printv([`You picked door`, door1]):
    UnpickedDoors:={1, 2, 3} minus {door1}:
    MontyDoors:=UnpickedDoors minus {1}:
    MontyDoor:=MontyDoors[1]:
    printv([`Monty shows you that door`, MontyDoor, `has a goat`]):
    printv([`and he kindly offers you to switch your`]):
    printv([`choice, and you agree`]):
    printv([`Your new choice`]):
    
    ND:={1, 2, 3} minus {door1, MontyDoor}:
    ND:=ND[1]:
    
    printv([`You decided to change your guessed door to door`, ND]):
    if ND = 1 then
        printv([`Congratulations! You won a Porche!`]):
        return true:
    else
        printv([`Monty tricked you; have fun with the goat`]):
        return false:
    fi:
end:

##PROBLEM 3##

#MHgs(p): Silent version of MHg
MHgs:=proc(p)
    return MHg(p, false):
end:

#Experiment with MHgs for various p
#N=number of iterations
MHgsExperiment:=proc(pset, N:=1000) local p, ret, count, n:
    ret:=[]:
    for p in pset do
        count:=0:
        for n from 1 to N do
            if MHgs(p) then
                count:=count + 1:
            fi:
        od:
        ret:=[op(ret), [p, evalf(count/N)]]:
    od:
end:

#Output was [[0, 1.], [1/10, 0.8860000000], [1/5, 0.7820000000], [3/10, 0.6860000000], [2/5, 0.5960000000], [1/2, 0.4960000000], [3/5, 0.4250000000], [7/10, 0.2970000000], [4/5, 0.1940000000], [9/10, 0.1040000000], [1, 0.]], so confirmed empirically

##PROBLEM 4##

#SDsim(N): simulates (kind of) SD(1/N,1/(N-1), 0), i.e. in a
#population of N people (whose names are 1,2, ...N) there is
#exactly one sick person (but no one knows who he or she is),
#each equally likely to be the sick one, and he or she would
#definitely be diagonsed sick (since FN=0), and out of the N-1
#healthy people, exactly one would be diagnosed (falsely) sick.
#So in every run one person would be real sick, and one person
#would be healthy-but-declared-sick. The ouput is a pair of
#(distint) integers [rs,fs].
SDsim:=proc(N) local st, i, rs:
    st:={seq(i, i=1..N)}:
    rs:=st[rand(1..nops(st))()]:
    st:=st minus {rs}:
    return [rs, st[rand(1..nops(st))()]]:
end:

#By running it many times, I have convinced myself that each
#person would be roughly an equal number of times really sick
#and falsely sick.

##PROBLEM 5##

#BDg(k, N, r): returns true if you can't find two persons,
#out of the k people that are <=r apart.
#k people, N days
BDg:=proc(k, N, r) local i, B, BDiffs:
    B:=[seq(rand(1..N)(), i=1..k)]:
    B:=sort(B):
    BDiffs:={seq(B[i]-B[i-1], i=2..nops(B)), B[1]-B[-1]+N}:
    return evalb(min(BDiffs) > r):
end:

#Run BDg M times, return proportion of true
BDgsim:=proc(k, N, r, M) local count, i:
    count:=0:
    for i from 1 to M do
        if BDg(k, N, r) then
            count:=count + 1:
        fi:
    od:
    return evalf(count/M):
end:

##PROBLEM 6##

#The number of ways to have no collisions (where a collision is
#two birthdays within r days of each other) is
#N*(N-2r-1)*(N-4r-2)*...*(N-2*(k-1)*r-(k-1))
#since each birthday creates a "window of exclusion" around it
#There are N^k total configurations, so the probability of
#no collisions is
#product(N-(2*r+1)*i, i=0..k-1)/N^k
#
#I don't know a way to remove the product and express this
#in terms of factorials

#Program implementing this formula
BDgf:=proc(k, N, r) local i:
    return evalf(mul(N-(2*r+1)*i, i=0..k-1)/N^k)
end:

#This formula checks out experimentally

##PROBLEM 7##

#IsDerangement(pi): inputs a permutation pi, and outputs true
#if the permutation is a derangement (i.e. has no fixed points)
#and false otherwise.
#For example IsDerangement([2,3,4,1]) should be true but
#IsDerangement([2,4,3,1]) should be false.
IsDerangement:=proc(pi) local i:
    for i from 1 to nops(pi) do
        if pi[i] = i then
            return false:
        fi:
    od:
    return true:
end:

#HowLikelyIsDer(n, K): inputs positive integers n and K, and by
#using combinat[randperm](n) and IsDerangement(pi), K times,
#counts the ratio of those that turned out to be derangements to
#the total K, in floating-point.
HowLikelyIsDer:=proc(n, K) local i, count:
    count:=0:
    for i from 1 to K do
        if IsDerangement(randperm(n)) then
            count:=count + 1:
        fi:
    od:
    return evalf(count/K):
end:

#HowLikelyIsDer(100,2000) outputted 0.3635000000 

#The number of permutations of length n with
#at least k fixed points is
#(n choose k)*(n-k)!=n!/k!
#since we can choose the k fixed points and then permute the rest
#arbitrarily. By Inclusion-Exclusion, the number of derangements
#on n elements is equal to
#sum((-1)^k*n!/k!, k=0..n)
#Since there are n! permutations, the probability that pi is
#a derangement is this divided by n!, namely
#sum((-1)^k/k!, k=0..n)
#This approaches 1/e as n approaches infinity.