######################################################################
##Urn.txt: Save this file as  Urn.txt                                #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read Urn.txt<Enter>                                #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#DoronZeil at gmail dot com                                          #
######################################################################
 
#Created: Dec. 29, 2017
 
print(`Created:  Dec. 29, 2017`):
print(` This is Urn.txt `):
print(`It is one of the packages that accompany the article `):
print(`How Rounds Should You Expect in Urn Solitaire?`):
print(`by Shalosh B. Ekhad and Doron Zeilberger`):

print(`and also available from Zeilberger's website`):
print(``):
print(`Please report bugs to DoronZeil at gmail dot com `):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://sites.math.rutgers.edu/~zeilberg/  .`):

print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type Ezra1();, for help with`):
 print(`a specific procedure, type Ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type Ezra();, for help with`):
 print(`a specific procedure, type Ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

with(combinat):

read `GuessHolo2.txt`:

Ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: CC, Wt, Wt1 `):
 print(``):

else
Ezra(args):
fi:

end:

Ezra:=proc()

if args=NULL then
 print(`The main procedures are:  GWalks, MC, PrW, PtWx, PrWG, RandG, Walks `):
 print(` `):


elif nops([args])=1 and op(1,[args])=CC then
print(`CC(L, WhB, BlB): inputs a list of elements in WhB union BlB, outputs the number of color changes`):
print(`Try: CC([w1,b1,w2,w3,b2],{w1,w2,w3},{b1,b2});`):

elif nops([args])=1 and op(1,[args])=GWalks then
print(`GWalks(a,b,c,d): the set of walks from [a,b] and [c,d] corresponding to the number of white balls and black balls`):
print(`in an urn that starts with a White balls and b Black balls such that c/d>=a/b. The walks from [a,b] to [c,d]`):
print(`with c<=a d<=b with unit negative steps, staying at c/d>=a/b. Try:`):
print(`GWalks(5,4,2,1);`):

elif nops([args])=1 and op(1,[args])=MC then
print(`MC(m,n,L,K): simulates K Simple Urn Solitaire games, outputs the ratio of those that ended in white, and`):
print(`the avearge number of color changes, followed by the variance, and alpha coefficients up to the L-th. Try:`):
print(`MC(4,7,4,100);`):

elif nops([args])=1 and op(1,[args])=PrW then
print(`PrW(a,b,c,d): prob. that a walk from [a,b] makes to [c,d]. Try:`):
print(`PrW(4,3,2,1);`):

elif nops([args])=1 and op(1,[args])=PrWx then
print(`PrWx(a,b,c,d,x): prob. generating function for the set of  walks from [a,b] to [c,d] for the r.v. number of color changes. Try:`):
print(`PrWx(4,3,2,1,x);`):

elif nops([args])=1 and op(1,[args])=PrWG then
print(`PrWG(a,b,c,d): prob. that a walk from [a,b] makes to [c,d] and it always lines in c/d>=a/b. Try:`):
print(`PrWG(4,3,2,1);`):

elif nops([args])=1 and op(1,[args])=RandG then
print(`RandG(m,n,W,B): Given a urn with m White balls labelled W[1], ..., W[m], and n Black balls, labelled B[1], ..., B[n],`):
print(` outputs a random way of picking them in order. Try: `):
print(` RandG(4,6,W,B); `):

elif nops([args])=1 and op(1,[args])=Walks then
print(`Walks(a,b,c,d): the set of walks from [a,b] and [c,d] corresponding to the number of white balls and black balls`):
print(`in an urn that starts with a White balls and b Black balls. The walks from [a,b] to [c,d]`):
print(`with c<=a d<=b with unit negative steps. Try:`):
print(`Walks(5,4,2,1);`):

elif nops([args])=1 and op(1,[args])=Wt then
print(`Wt(L): the weight of a walk L. Try:`):
print(`Wt([[3,2],[3,1],[2,1],[1,1]]);`):

elif nops([args])=1 and op(1,[args])=Wt1 then
print(`Wt1(pt1,pt2): the weight of a step from pt1 to pt2`):

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

##From Histabrut
#AveAndMoms(f,x,N): Given a probability generating function
#f(x) (a polynomial, rational function and possibly beyond)
#returns a list whose first entry is the average 
#(alias expectation, alias mean)
#followed by the variance, the third moment (about the mean) and
#so on, until the N-th moment (about the mean).
#If f(1) is not 1, than it first normalizes it by dividing
#by f(1) (if it is not 0) .
#For example, try:
#AveAndMoms(((1+x)/2)^100,x,4);

AveAndMoms:=proc(f,x,N) local mu,F,memu1,gu,i:
mu:=simplify(subs(x=1,f)):

if mu=0 then
print(f, `is neither a prob. generating function nor can it be made so`):
RETURN(FAIL):
fi:

F:=f/mu:


memu1:=simplify(subs(x=1,diff(F,x))):

gu:=[memu1]:

F:=F/x^memu1:

F:=x*diff(F,x):
for i from 2 to N do
 F:=x*diff(F,x):
 gu:=[op(gu),simplify(subs(x=1,F))]:
od:

gu:

end:

#Alpha(f,x,N): Given a probability generating function
#f(x) (a polynomial, rational function and possibly beyond)
#returns a list, of length N, whose 
#(i) First entry is the average
#(ii): Second entry is the variance
#for i=3...N, the i-th entry is the so-called alpha-coefficients
#that is the i-th moment about the mean divided by the
#variance to the power i/2 (For example, i=3 is the skewness
#and i=4 is the Kurtosis)
#If f(1) is not 1, than it first normalizes it by dividing
#by f(1) (if it is not 0) .
#For example, try:
#Alpha(((1+x)/2)^100,x,4);
Alpha:=proc(f,x,N) local gu,i:
 gu:=AveAndMoms(f,x,N):
 if gu=FAIL then
  RETURN(gu):
 fi:

if gu[2]=0 then
print(`The variance is 0`):
  RETURN(FAIL):
fi:

[gu[1],gu[2],seq(gu[i]/gu[2]^(i/2),i=3..N)]:

end:

##End From Histabrut

#RandG(m,n,W,B): Given a urn with m White balls labelled W[1], ..., W[m], and n Black balls, labelled B[1], ..., B[n],
#outputs a random way of picking them in order. Try:
#RandG(4,6,W,B);
RandG:=proc(m,n,W,B) local gu,S,i,gu1:
 S:={seq(W[i],i=1..m),seq(B[i],i=1..n)}:
gu:=[]:

while S<>{} do
 gu1:=S[rand(1..nops(S))()]:
 gu:=[op(gu),gu1]:
 S:=S minus {gu1}:
od:

gu:

end:

#CC(L, WhB, BlB): inputs a list of elements in WhB union BlB, outputs the number of color changes
#Try: CC([w1,b1,w2,w3,b2],{w1,w2,w3},{b1,b2});
CC:=proc(L, WhB, BlB) local i,c:

c:=0:
if convert(L,set) minus (WhB union BlB)<>{} then
   RETURN(FAIL):
fi:

for i from 1 to nops(L)-1 do
 if (member(L[i],WhB) and  member(L[i+1],BlB)) or (member(L[i],BlB) and  member(L[i+1],WhB)) then
   c:=c+1:
 fi:
od:
c:
end:





#MC(m,n,L,K): simulates K Simple Urn Solitaire games, outputs the ratio of those that ended in white, and
#the avearge number of color changes, followed by the variance, and alpha coefficients up to the L-th. Try
#MC(4,7,4,100);
MC:=proc(m,n,L,K) local i,g,x,c,W,B,WhB,BlB,mu:

 WhB:={seq(W[i],i=1..m)}:
 BlB:={seq(B[i],i=1..n)}:
c:=0:

mu:=0:

for i from 1 to K do
g:=RandG(m,n,W,B):
 if member(g[-1],WhB) then
   c:=c+1:
 fi:
mu:=mu+x^CC(g,WhB,BlB):
od:

evalf([c/K,Alpha(mu,x,L)]):

end:


#Walks(a,b,c,d): the set of walks from [a,b] and [c,d] corresponding to the number of white balls and black balls
#in an urn that starts with a White balls and b Black balls. The walks from [a,b] to [c,d]
#with c<=a d<=b with unit negative steps. Try:
#Walks(5,4,2,1);
Walks:=proc(a,b,c,d) local gu,w,gu1:
option remember:

if c>a or c<0 or d>b or d<0 then
 RETURN({}):
fi:

if c=a and d=b then
 RETURN({[[a,b]]}):
fi:


gu1:=Walks(a,b,c+1,d):
gu:={seq([op(w),[c,d]],w in gu1)}:
gu1:=Walks(a,b,c,d+1):
gu:=gu union {seq([op(w),[c,d]],w in gu1)}:
gu:
end:

#GWalks(a,b,c,d): the set of walks from [a,b] and [c,d] corresponding to the number of white balls and black balls
#in an urn that starts with a White balls and b Black balls and such that the ratio of white balls to black balls
# is >=a/b
# The walks from [a,b] to [c,d]
#with c<=a d<=b with unit negative steps staying in the region bx-ay>=0. Try:
#Walks(5,4,2,1);
GWalks:=proc(a,b,c,d) local gu,w,gu1:
option remember:

if c>a or c<0 or d>b or d<0 or b*c-a*d<0 then
 RETURN({}):
fi:

if c=a and d=b then
 RETURN({[[a,b]]}):
fi:


gu1:=GWalks(a,b,c+1,d):
gu:={seq([op(w),[c,d]],w in gu1)}:
gu1:=GWalks(a,b,c,d+1):
gu:=gu union {seq([op(w),[c,d]],w in gu1)}:
gu:
end:


#Wt1(pt1,pt2): the weight of a step from pt1 to pt2
Wt1:=proc(pt1,pt2) local a,b:
a:=pt1[1]:
b:=pt1[2]:

if pt2=[a-1,b] then
 a/(a+b):
elif pt2=[a,b-1] then
 b/(a+b):
else
FAIL:
fi:
end:

#Wt(L): the weight of a walk L
Wt:=proc(L) local i:
mul(Wt1(L[i],L[i+1]),i=1..nops(L)-1):
end:



#PrW(a,b,c,d): prob. that a walk from [a,b] makes to [c,d]. Try:
#PrW(4,3,2,1);
PrW:=proc(a,b,c,d) option remember:

if c>a or c<0 or d>b or d<0 then
 RETURN(0):
fi:

if c=a and d=b then
  RETURN(1):
fi:

(c+1)/(c+1+d)*PrW(a,b,c+1,d)+(d+1)/(c+1+d)*PrW(a,b,c,d+1):

end:

#PrWG(a,b,c,d): prob. that a walk from [a,b] makes to [c,d] and it is good. Try:
#PrWG(4,3,2,1);
PrWG:=proc(a,b,c,d) option remember:

if c>a or c<0 or d>b or d<0 or b*c-a*d<0 then
 RETURN(0):
fi:

if c=a and d=b then
  RETURN(1):
fi:

(c+1)/(c+1+d)*PrWG(a,b,c+1,d)+(d+1)/(c+1+d)*PrWG(a,b,c,d+1):

end:


#PrWW(a,b,c,d): prob. that a walk from [a,b] makes to [c,d] and ends with a White Ball. Try:
#PrWW(4,3,2,1);
PrWW:=proc(a,b,c,d) option remember:

if c>a or c<0 or d>b or d<0 or c=a and d=b then
 RETURN(0):
fi:

if c=a-1 and d=b then
  RETURN(a/(a+b)):
fi:

(c+1)/(c+1+d)*(PrWW(a,b,c+1,d)+PrWB(a,b,c+1,d)):

end:

#PrWB(a,b,c,d): prob. that a walk from [a,b] makes to [c,d] and ends with a Black Ball. Try:
#PrWB(4,3,2,1);
PrWB:=proc(a,b,c,d) option remember:

if c>a or c<0 or d>b or d<0 or c=a and d=b then
 RETURN(0):
fi:

if c=a and d=b-1 then
  RETURN(b/(a+b)):
fi:

(d+1)/(c+1+d)*(PrWW(a,b,c,d+1)+PrWB(a,b,c,d+1)):

end:

#PrWWx(a,b,c,d,x): prob. generating function of the r.v.
#Number of color changes for  walks from [a,b] makes to [c,d] and ends with a White Ball. Try:
#PrWW(4,3,2,1);
PrWWx:=proc(a,b,c,d,x) option remember:

if c>a or c<0 or d>b or d<0 or c=a and d=b then
 RETURN(0):
fi:

if c=a-1 and d=b then
  RETURN(a/(a+b)):
fi:

expand((c+1)/(c+1+d)*(PrWWx(a,b,c+1,d,x)+x*PrWBx(a,b,c+1,d,x))):

end:

#PrWBx(a,b,c,d,x): prob. generating function for the set of  walks from [a,b] makes to [c,d] and ends with a Black Ball. Try:
#PrWBx(4,3,2,1,x);
PrWBx:=proc(a,b,c,d,x) option remember:

if c>a or c<0 or d>b or d<0 or c=a and d=b then
 RETURN(0):
fi:

if c=a and d=b-1 then
  RETURN(b/(a+b)):
fi:

expand((d+1)/(c+1+d)*(x*PrWWx(a,b,c,d+1,x)+PrWBx(a,b,c,d+1,x))):

end:

#PrWx(a,b,c,d,x): prob. generating function for the set of  walks from [a,b] makes to [c,d]. Try:
#PrWx(4,3,2,1,x);
PrWx:=proc(a,b,c,d,x) option remember:
expand(PrWWx(a,b,c,d,x)+PrWBx(a,b,c,d,x)):
end: