###################################################################### ##SHEPP.txt: Save this file as SHEPP.txt # ## To use it, stay in the # ##same directory, get into Maple (by typing: maple ) # ##and then type: read SHEPP # ##Then follow the instructions given there # ## # ##Written by Doron Zeilberger, Rutgers University , # #zeilberg at math dot rutgers dot edu # ###################################################################### #Created: Aug.-Sept. 2012 print(`Created: Aug.-Sept. 2012`): print(` This is SHEPP.txt `): print(`It is one of the packages that accompany the article `): print(`and also available from Zeilberger's website`): print(``): print(`Please report bugs to zeilberg at math dot rutgers dot edu`): print(``): print(`The most current version of this package and paper`): print(` are available from`): print(`http://www.math.rutgers.edu/~zeilberg/ .`): print(`For a list of the procedures type ezra();, for help with`): print(`a specific procedure, type ezra(procedure_name); .`): print(``): with(combinat): ezra1:=proc() if args=NULL then print(` The supporting procedures are: Moms, NorMoms`): print(``): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(`The main procedures are: B, beta, beta1, betaN, Larry1, Larry2,`): print(` Seqb, SeqB, V, Vra, Vrax, Vx `): print(` `): elif nops([args])=1 and op(1,[args])=B then print(`B(m,p): V(m,p0*binomial(n+p,p), always an integer thanks to `): print(`William C. Boyce. Try: B(4,4);`): elif nops([args])=1 and op(1,[args])=beta then print(`beta(p): If there are p plus balls, the largest m for which`): print(` V(m,p)>0, try beta(5); `): elif nops([args])=1 and op(1,[args])=beta1 then print(`beta1(p): The smallest m for which V(m,p)=0 follwed`): print(`by true if stop and go are equally good`): elif nops([args])=1 and op(1,[args])=betaN then print(`betaN(p):evalf((beta(p)-p)/sqrt(2*p)):`): elif nops([args])=1 and op(1,[args])=Larry1 then print(`Larry1(N): the list whose p-th item is the smallest m`): print(`for which the expected gain in the Shepp Urn game is 0`): print(`followed by the list of those for which it is OK to go`): print(`or stop. Larry Shepp conjectured that the only value is m=2 p=1`): elif nops([args])=1 and op(1,[args])=Larry2 then print(`Larry2(N): the list of pairs [[m,p],ExpectedLoss]`): print(`ranked according to the loss in continuing`): print(`where m is the smallest m for which [m,p]`): print(`for which the expected gain in the Shepp Urn game is 0`): print(`followed by the expected (non=positive) gain if you stupidly`): print(`keep playing`): elif nops([args])=1 and op(1,[args])=Moms then print(`Moms(f,x,R): all the R-th first moments (about the mean)`): elif nops([args])=1 and op(1,[args])=NorMoms then print(`NorMoms(f,x,R): all the Normalized R-th first moments (about the mean)`): elif nops([args])=1 and op(1,[args])=Seqb then print(`Seqb(N): the list of length N whose p-th entry`): print(`gives you beta(p) followed the list of`): print(`of B(p,p) for p from 1 to N`): print(`followed by the list of evalf of (beta(p)-p)/sqrt(2*p)`): print(`followed by the list of evalf of B(p,p)/binomial(2*p,p)/sqrt(p)`): print(`Try Seqb(100);`): elif nops([args])=1 and op(1,[args])=SeqB then print(`SeqB(N): the list of length N whose n-th entry`): print(`gives you the list of length n+1 whose m-th entry is`): print(`B(m,n-m); `): elif nops([args])=1 and op(1,[args])=V then print(`V(m,p): The expected gain in the Shepp Urn game if`): print(`right now there are m minus balls and p plus balls in the urn`): print(`Try V(3,3);`): elif nops([args])=1 and op(1,[args])=Vra then print(`Vra(m,p,m0,CutOff): The expected gain in`): print(`a risk-averse Shepp Urn game. If`): print(`right now there are m minus balls and p plus balls in the urn`): print(`and the play is risk averse, at any stage the player stops`): print(`unless (i) the expected gain is positive AND`): print(`(ii) the probability of losing >= m0 dollars`): print(`is more than CutOff. For example try:`): print(`Try Vra(11,10,1,1/5);`): elif nops([args])=1 and op(1,[args])=Vrax then print(`Vrax(m,p,m0,CutOff,x): the prob. generating function in x,`): print(`for risk-averse Shepp Urn game, if`): print(`right now there are m minus balls and p plus balls in the urn`): print(`and the play is risk averse, at any stage the player stops`): print(`unless (i) the expected gain is positive AND`): print(`(ii) the probability of losing >= m0 dollars`): print(`is more than CutOff. For example try:`): print(`Try Vrax(11,10,1,1/5,x);`): elif nops([args])=1 and op(1,[args])=Vx then print(`Vx(m,p,x): the prob. generating function for the money gained`): print(`following the maximize expectation strategy,`): print(`Try Vx(10,10,x);`): else print(`There is no ezra for`,args): fi: end: #V(m,p): The expected gain in the Shepp Urn game if #right now there are m minus balls and p plus balls in the urn #Try V(3,3); V:=proc(m,p) option remember: if m=0 then RETURN(p): elif p=0 then RETURN(0): fi: max(0,p/(m+p)*V(m,p-1)+m/(p+m)*V(m-1,p)+(p-m)/(p+m)): end: #beta(p): If there are p plus balls the largest m for which #V(m,p)>0, try beta(5); beta:=proc(p) local m: for m from p while V(m,p)>0 do od: m-1: end: #betaB(p): If there are p plus balls the largest m for which #B(m,p)>0, try beta(5); betaB:=proc(p) local m: for m from p while V(m,p)>0 do od: m-1: end: #betaN(p):evalf((beta(p)-p)/sqrt(2*p)): betaN:=proc(p): evalf((beta(p)-p)/sqrt(2*p)):end: Vf:=proc(m,p) option remember: if m=0 then RETURN(p): elif p=0 then RETURN(0): fi: max(0,evalf(p/(m+p)*Vf(m,p-1)+m/(p+m)*Vf(m-1,p)+(p-m)/(p+m))): end: betaf:=proc(p) local m: for m from p while Vf(m,p)>0 do od: m-1: end: betaNf:=proc(p): evalf((beta(p)-p)/sqrt(2*p)):end: #Moms(f,x,R): all the R-th first moments (about the mean) Moms:=proc(f,x,R) local f1,av,gu,i: av:=subs(x=1,diff(f,x)): f1:=expand(f/x^av): gu:=[av]: f1:=x*diff(f1,x): for i from 1 to R-1 do f1:=x*diff(f1,x): gu:=[op(gu),subs(x=1,f1)]: od: gu: end: #NorMoms(f,x,R): all the Normalized R-th first moments (about the mean) NorMoms:=proc(f,x,R) local f1,av,gu,i: av:=subs(x=1,diff(f,x)): f1:=expand(f/x^av): gu:=[av]: f1:=x*diff(f1,x): for i from 1 to R-1 do f1:=x*diff(f1,x): gu:=[op(gu),subs(x=1,f1)]: od: evalf([gu[1],sqrt(gu[2]),seq(gu[i]/gu[2]^(i/2),i=3..nops(gu))]): end: #beta1(p): The smallest m for which V(m,p)=0 follwed #by true if stop and go are equally good beta1:=proc(p) local m: for m from p while V(m,p)>0 do od: if p/(m+p)*V(m,p-1)+m/(p+m)*V(m-1,p)+(p-m)/(p+m)=0 then RETURN([m,true]): else RETURN([m,false]): fi: end: #Larry1(N): the list whose p-th item is the smallest m #for which the expected gain in the Shepp Urn game is 0 #followed by the list of those for which it is OK to go #or stop. Larry Shepp conjectured that the only value is m=2 p=1 Larry1:=proc(N) local gu,mu,p,ka: gu:=[]: mu:=[]: for p from 1 to N do ka:=beta1(p): gu:=[op(gu),ka[1]]: if ka[2] then mu:=[op(mu),ka[1]]: fi: od: gu,mu: end: #beta2(p): The smallest m for which V(m,p)=0 followed #by how much you lose if you do continue beta2:=proc(p) local m: for m from p while V(m,p)>0 do od: [m,p/(m+p)*V(m,p-1)+m/(p+m)*V(m-1,p)+(p-m)/(p+m)]: end: #Larry2(N): the list of pairs [[m,p],ExpectedLoss] #ranked according to the loss in continuing #where m is the smallest m for which [m,p] #for which the expected gain in the Shepp Urn game is 0 #followed by the expected (non=positive) gain if you stupidly #keep playing Larry2:=proc(N) local p,T,gu1,gu: gu:=[seq(beta2(p),p=1..N)]: for p from 1 to nops(gu) do T[-gu[p][2]]:=p: od: gu1:=[seq(-gu[p][2],p=1..nops(gu))]: gu1:=sort(gu1): [seq( [T[gu1[p]],beta2(T[gu1[p]]) ],p=1..nops(gu1))]: end: #Vx(m,p,x): the prob. generating function for the money gained #following the maximize expectation strategy, #Try Vx(10,10,x); Vx:=proc(m,p,x) option remember: if m=0 then RETURN(x^p): fi: if p=0 then RETURN(1): fi: if V(m,p)=0 then RETURN(1): else RETURN(expand(p/(m+p)*x*Vx(m,p-1,x)+m/(p+m)*(1/x)*Vx(m-1,p,x))): fi: end: #Vrax(m,p,m0,CutOff,x): the prob. generating function in x, #for risk-averse Shepp Urn game, if #right now there are m minus balls and p plus balls in the urn #and the play is risk averse, at any stage the player stops #unless (i) the expected gain is positive AND #(ii) the probability of losing >= m0 dollars #is more than CutOff. For example try: #Try Vrax(11,10,1,1/5,x); Vrax:=proc(m,p,m0,CutOff,x) local gu,i: option remember: if m=0 then RETURN(x^p): elif p=0 then RETURN(1): fi: gu:=expand(p/(m+p)*x*Vrax(m,p-1,m0,CutOff,x)+ m/(p+m)*(1/x)*Vrax(m-1,p,m0,CutOff,x)): if subs(x=1,gu)>0 and add(coeff(gu,x,i),i=ldegree(gu,x)..-m0)= m0 dollars #is more than CutOff. For example try: #Try Vra(11,10,1,1/5); Vra:=proc(m,p,m0,CutOff) local x: subs(x=1,diff(Vrax(m,p,m0,CutOff,x),x)): end: #B(m,p): V(m,p0*binomial(n+p,p), always an integer thanks to #William C. Boyce. Try: B(4,4); B:=proc(m,p) option remember: if m<0 or p<0 or (m=0 and p=0) then 0: else max(0,binomial(m+p-1,m)-binomial(m+p-1,p)+B(m-1,p)+B(m,p-1)): fi: end: #SeqB(N): the list of length N whose n-th entry #gives you the list of length n+1 whose m-th entry is #B(m,n-m); SeqB:=proc(N) local gu,p,m,lu1,lu,mu: mu:=[1,1]: gu:=[mu]: for p from 2 to N do lu1:=p: lu:=[]: for m from 1 while lu1>0 do lu:=[op(lu),lu1]: if m+1<=nops(mu) then lu1:=max(0,binomial(m+p-1,m)-binomial(m+p-1,p)+mu[m+1]+lu1): else lu1:=max(0,binomial(m+p-1,m)-binomial(m+p-1,p)+lu1): fi: od: gu:=[op(gu),lu]: mu:=lu: od: gu: end: #Seqb(N): the list of length N whose p-th entry #gives you beta(p) followed the list of #of B(p,p) for p from 1 to N #Try Seqb(100); Seqb:=proc(N) local gu1,gu2,p,m,lu1,lu,mu,ru1,ru2: mu:=[1,1]: gu1:=[1]: gu2:=[1]: for p from 2 to N do lu1:=p: lu:=[]: for m from 1 while lu1>0 do lu:=[op(lu),lu1]: if m+1<=nops(mu) then lu1:=max(0,binomial(m+p-1,m)-binomial(m+p-1,p)+mu[m+1]+lu1): else lu1:=max(0,binomial(m+p-1,m)-binomial(m+p-1,p)+lu1): fi: od: gu1:=[op(gu1),nops(lu)-1]: gu2:=[op(gu2),lu[p+1]]: mu:=lu: od: gu1,gu2: ru1:=[seq(evalf((gu1[p]-p)/sqrt(2*p)),p=1..N)]: ru2:=[seq(evalf(gu2[p]/binomial(2*p,p)/sqrt(p)),p=1..N)]: [gu1,gu2,ru1,ru2]: end: