Digits:=500:
interface(displayprecision=10):

####################################################################################################
#
#GeneralizedShepp.txt: Save this file as GeneralizedShepp.txt
#In order to use it, point Maple to the directory this file is in and type: read "GeneralizedShepp.txt":
#You may then follow the instrustions on screen.
#
#Written by Richard Voepel, Rutgers University
#rzv2@math.rutgers.edu
#
####################################################################################################

print(`This is GeneralizedShepp.txt`):
print(`It is a small Maple package that accompanies the following article:`):
print(`EXPERIMENTAL RESULTS AND GENERALIZATIONS FOR THE STOPPING PROBLEM OF SHEPP’S URN`):
print(`by Richard Voepel`):
print(``):
print(`Please report any bugs to Richard Voepel at: rzv2@math.rutgers.edu`):
print(`The most current version of this package and the associated paper can be found at:`):
print(`http://sites.math.rutgers.edu/~rzv2/Shepp`):
print(``):
print(`For a list of procedures, type the following: Help();`):

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

Help:=proc()

if nargs=0 then
print(`Contains primary procedures: RiskFuncGen, PGFTableGen, MmntTableGen`):
print(`For information about a primary procedure P, type Help(P).`):
print(`For auxiliary procedures, type HelpAux().`):

elif nargs=1 and args[1]=RiskFuncGen then
print(`RiskFuncGen(d,q,x): Returns a functional on PGFs that encodes the probability of losing d or more dollars`):
print(`being less than q. For example, try typing the following: F:=RiskFuncGen(1,1/8,x):`):

elif nargs=1 and args[1]=PGFTableGen then
print(`PGFTableGen(M,P,F,x): Returns a table of PGFs for urn random variables U(m,p), for all m+p<=M+P,`):
print(`based on the strategy given by a single functional F on PGFs. The PGFs are always Laurent polynomials`):
print(`of the indeterminate x. The boolean conditional for the strategy is of the form F(g)>0. However,`):
print(`the base cases for the urns are allways U(0,p)=x^p and U(m,0)=1 by convention.`):
print(`For example, try typing the following: F:=RiskFuncGen(1,1/8,x): P:=PGFTableGen(20,20,F,x):`):

elif nargs=1 and args[1]=MmntTableGen then
print(`MmntTableGen(PGFTable,M,P,r,x): Returns a table M, where the entry M[m,p] is a list of the first r standardized`):
print(`moments of the PGF located in PGFTable[m,p], starting from the 0th moment, for all m+p<=M+P. The indeterminate x`):
print(`must be the same indeterminate used in PGFTable, and PGFTable must be defined for all m+p<=M+P. The value of r`):
print(`must be greater than or equal to 3.`):
print(`For example, try typing the following: F:=RiskFuncGen(1,1/8,x): P:=PGFTableGen(20,20,F,x): M:=MmntTableGen(P,20,20,6,x):`):

else
print(`There is no help for`, args):
fi:

end:

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

HelpAux:=proc()

if nargs=0 then
print(`Contains auxiliary procedures: CPF, FactorialMoment, CentralizedMoment, StandardizedMoment`):
print(`For information about an auxiliary procedure P, type HelpAux(P).`):

elif nargs=1 and args[1]=CPF then
print(`CPF(n): The Cantor Pairing Function, which maps naturals {0,1,...} to pairs of naturals.`):
print(`For example, try typing the following: CPF(17);`):

elif nargs=1 and args[1]=FactorialMoment then
print(`FactorialMoment(f,r,x): Given a probability generating function f in the indeterminate x, returns`):
print(`the first r (starting from 0) factorial moments.`):
print(`For example, try typing the following: FactorialMoment(1/2+1/2*x,6,x);`):

elif nargs=1 and args[1]=CentralizedMoment then
print(`CentralizedMoment(l): Takes a list l of factorial moments, and returns a list of centralized moments.`):
print(`For example, try typing the following: CentralizedMoment(FactorialMoment(1/2+1/2*x,6,x));`):

elif nargs=1 and args[1]=StandardizedMoment then
print(`StandardizedMoment(l): Take a list of l centralized moments, and if the length of l is at least 3,`):
print(`returns a list of standardized moments. If the variance is too close to 0, the original argument`):
print(`will be returned instead.`):
print(`For example, try typing the following: StandardizedMoment(CentralizedMoment(FactorialMoment(1/2+1/2*x,6,x)));`):

else
print(`There is no help for`, args):
fi:

end:

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

#CPF(n), the Cantor Pairing Function, which maps naturals {0,1,...} to pairs of naturals.
CPF:=proc(n) local w,t,m,p:
w:=floor((sqrt(8*n+1)-1)/2):
t:=(w^2+w)/2:
m:=n-t:
p:=w-m:
RETURN([m,p]):
end:

#RiskFuncGen(q,d,x) returns a special form of functional on PGFs that encodes the probability of losing d or more dollars being less than q.
RiskFuncGen:=proc(d,q,x) local i: 
RETURN(f->q-subs(x=1,add(coeff(f,x,i)*x^i,i=ldegree(f)..-d))):
end:


#PGFTableGen(M,P,F,x) produces a table of probability generating functions for the Shepp Urns U(m,p), for all m+p<=M+P, based on the strategy given by a single functional F.
#The probability generating functions are functions of the indeterminate x, and are always Laurent polynomials.
#The boolean conditional for the strategy is of the form F(g)>0 for the given functional F on probability generating functions.
#Regardless of the conditional, the probability generating function for urns U(0,p) is x^p, and for U(m,0) is 1, to reflect that no strategy that is remotely rational would
#pass on an urn with only plus balls or draw from an urn with only negative balls.
PGFTableGen:=proc(M,P,F,x) local N,i,m,p,U:
N:=M+P+1:
for i from 0 to (N^2+N)/2-1 do
m:=CPF(i)[1]:
p:=CPF(i)[2]:
if evalb(p=0)then
U[m,p]:=1:
elif evalb(m=0)then
U[m,p]:=x^p:
else
U[m,p]:=expand(simplify(m/(p+m)*x^(-1)*U[m-1,p]+p/(p+m)*x*U[m,p-1])):
fi:
if evalb(F(U[m,p])<=0) then
U[m,p]:=1:
fi:
od:
RETURN(U):
end:


#MmntTableGen(PGFTable,M,P,r,x) produces a table of the first r (including 0) standardized moments corresponding to an input table of PGFs.
#The table of PGFs has entries [m,p] for all m+p<=M+P, and each PGF is in terms of x.
MmntTableGen:=proc(PGFTable,M,P,r,x) local T,N,i,m,p:
if r<3 then
RETURN(`ERROR: r>=3 required.`)
fi:
N:=M+P+1:
for i from 0 to (N^2+N)/2-1 do
m:=CPF(i)[1]:
p:=CPF(i)[2]:
T[m,p]:=StandardizedMoment(CentralizedMoment(FactorialMoment(PGFTable[m,p],r,x))):
od:
RETURN(T):
end:


FactorialMoment:=proc(f,r,x) local g,h,mu,list:
g:=diff(f,x):
mu:=subs(x=1,g):
h:=expand(simplify(f/x^(mu))):
list:=[subs(x=1,h)]:
while nops(list)<r do
list:=[op(list),subs(x=1,diff(h,x$nops(list)))]
od:
RETURN(list):
end:


CentralizedMoment:=proc(l) local list,m,k:
list:=[1]:
while nops(list)<nops(l) do
m:=add(Stirling2(nops(list),k)*l[k+1],k=1..nops(list)):
list:=[op(list),m]:
od:
RETURN(list):
end:


StandardizedMoment:=proc(l) local k,list:
if evalb(nops(l)<3) then
RETURN(`ERROR: Length >=3 required.`)
fi:
list:=[0$nops(l)]:
list[1]:=l[1]:
if evalb(abs(l[3])<10^(floor(-Digits/4))) then
RETURN(l):
fi:
for k from 2 to nops(l) do
list[k]:=l[k]/sqrt(l[3])^(k-1):
od:
RETURN(list):
end: