#!/usr/local/bin/maple
# -*- maplev -*-

# Nathaniel Shar
# HW 14
# Experimental Mathematics

# It is okay to link to this assignment on the course webpage.

# From c14.txt:

################################################################################
#March 11, 2013, Larry Shepp Urn
Help:=proc(): print(` Larry(m,p), LarryE(m,p) `):
print(`LarryTale(N),beta(p) , LarryPGF(m,p,x)`):
print(`PrDist(m,p) `):
print(` Stat(P,x,K) `):
end:

#LarryE(m,p): input non.neg. integers m and p
#representing a box (urn) with m (-1)dollar bills
#and p 1-one dollar bills, outputs the
#expected gain of the game (where at any given time)
#you pick at random one or the other dollars bills
#Pr. of picking a dollar bill is p/(m+p) and
#Pr. of picking a -1 dollar bill m/(m+p)
#And you quit whenever you want,
LarryE:=proc(m,p)
option remember:
max(Larry(m,p),0):
end:

#Larry(m,p): input non.neg. integers m and p
#representing a box (urn) with m (-1)dollar bills
#and p 1-one dollar bills, outputs the
#expected gain of the game (where at any given time)
#you pick at random one or the other dollars bills
#Pr. of picking a dollar bill is p/(m+p) and
#Pr. of picking a -1 dollar bill m/(m+p)
#And you quit whenever you want
Larry:=proc(m,p)
option remember:

if m=0 then
 RETURN(p):
fi:

if p=0 then
 RETURN(-1):
fi:

p/(m+p)*(1+LarryE(m,p-1))+ m/(m+p)*(-1+LarryE(m-1,p)):



end:

#LarryTable(N):
LarryTable:=proc(N) local m,p:
for p from N by -1 to 0 do
 lprint(seq(evalf(LarryE(m,p),3),m=0..N)):
od:




end:


#beta(p): The largest m for which LarryE(m,p) is >0
#beat(p)-p is asymtotically sqrt(p) times a (explicit) Constant
beta:=proc(p) local m:

for m from p while LarryE(m,p)>0 do od:
m-1:

end:


#LarrtPGF(m,p,x): The generalized "polynomial"
#whose exponets are numbers (rational) and
#where the coeff. of x^a is the prob. of exactly getting
#a dollars when you end the game (but following the
#stupid max. exp. strategy)
LarryPGF:=proc(m,p,x)
option remember:

if LarryE(m,p)=0 then
  RETURN(1):
fi:

if m=0 then
 RETURN(x^p):
fi:



expand(p/(m+p)*x*LarryPGF(m,p-1,x)+
m/(m+p)/x*LarryPGF(m-1,p,x)):



end:


#PrDist(m,p): the list of pairs [i,p(i)] where
#the prob. of getting i is p(i)
PrDist:=proc(m,p) local olivia,x,i:

olivia:=LarryPGF(m,p,x):

[seq([i, coeff(olivia,x,i)],
i=ldegree(olivia,x)..degree(olivia,x))]:

end:

#Stat(P,x,K): inputs a prob. g.f. P in x and a pos. integer K
#outputs a list of length K where the first entry if the
#average, the second the variance, and the r-th is the s-called
#alpha coefficiient m[r]/m[2]^(r/2) where m[r] is the
#r-th moment about the mean
#try:
#evalf(Stat(((1+x)/2)^1000,x,10));
#Taken from homework problem, not done in class
Stat:=proc(P,x,K) local av,va,L,P1,r:

P1:=P/subs(x=1,P):

av:=subs(x=1,diff(P1,x)):
P1:=expand(P1/x^av):
P1:=expand(x*diff(P1,x)):
P1:=expand(x*diff(P1,x)):
va:=subs(x=1,P1):

L:=[av,va]:

for r from 3 to K do
P1:=expand(x*diff(P1,x)):
L:=[op(L),subs(x=1,P1)/va^(r/2)]:
od:

L:
end:

################################################################################

#############
# Problem 1 #
#############

SheppEst := proc(n):
    return (beta(n)-n)/sqrt(2*n):
end:

# The values go up and down, but averaging the last 100 that we can
# get before maple refuses to recurse further, we get 0.83, so that's
# my estimate.

#############
# Problem 2 #
#############


PrWinBreakEvenLose := proc(m,p) local GF,x:
    GF := LarryPGF(m,p,x):
    return [add(coeff(GF, x, i), i=1..degree(GF)), coeff(GF, x, 0), add(coeff(GF, x, i), i=ldegree(GF)..-1)]:
end:

# It appears that the probability of losing when m = beta(p) is small
# and roughly decreases as p increases. You end up with a high
# probability of gaining 1 and a small probability of losing, but you
# may lose a large amount.

#############
# Problem 3 #
#############

VerifyTheorem3point1a := proc(N) local m,p,s:
    for p from 0 to N do:
        for m from 0 to beta(p) do:
            if LarryE(m, p+1) < LarryE(m,p) + 1/(m+p+1) then:
                return false:
            fi:
        od:
    od:
    return true:
end:

VerifyTheorem3point1b := proc(N) local m,p,s:
    for p from 0 to N do:
        for m from 0 to N do:
            if LarryE(m, p+1) > LarryE(m,p) + 1 then:
                return false:
            fi:
        od:
    od:
    return true:
end:

VerifyTheorem3point2a := proc(N) local m,p,s:
    for p from 0 to N do:
        for m from 0 to beta(p)-1 do:
            if LarryE(m, p) < LarryE(m+1,p) + 1/(m+p+1) then:
                return false:
            fi:
        od:
    od:
    return true:
end:

VerifyTheorem3point2b := proc(N) local m,p,s:
    for p from 0 to N do:
        for m from 0 to N do:
            if LarryE(m, p) > LarryE(m+1,p) + 1 then:
                return false:
            fi:
        od:
    od:
    return true:
end:

VerifyTheorem3point8 := proc(N) local m,p,s:
    for p from 0 to N do:
        for m from 0 to beta(p) do:
            if LarryE(m+1, p+1) <= LarryE(m,p) then:
                return false:
            fi:
        od:
    od:
    return true:
end:

# They are all true.

#############
# Problem 4 #
#############

# LarryE(2,n) = p*(p-1)/(p+1)
# LarryE(3,n) = (p-1)*(p^3+4*p^2-12)/((p+1)*(p+2)*(p+3))