###################################################################### ##PosPileGames.txt: Save this file as PosPileGames.txt # ## To use it, stay in the # ##same directory, get into Maple (by typing: maple ) # ##and then type: read PosPileGames.txt # ##Then follow the instructions given there # ## # ##Written by Doron Zeilberger, Rutgers University , # #DoronZeil at gmail dot com # ###################################################################### #Created: July 2019 print(`Created: July , 2019`): print(` This is PosPileGames.txt `): print(`It is one of the Maple packages that accompany the article `): print(` "A Multi-Computational Exploration of Some Games of Pure Chance" `): print(`by Thotsaporn Aek Thanatipanonda 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 Simulation procedures type ezraSi();, for help with`): print(`a specific procedure, type ezra(procedure_name); .`): print(``): print(`---------------------------------------`): print(`---------------------------------------`): print(`For a list of the STORY procedures type ezraSt();, 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): ezraSi:=proc() if args=NULL then print(` The Simulations procedures are: OneGame, ManyGames `): print(``): else ezra(args): fi: end: ezraSt:=proc() if args=NULL then print(` The STORY procedures are: MomsAMv, MomsFnkV, SSqSeqRecV `): print(``): else ezra(args): fi: end: ezra1:=proc() if args=NULL then print(` The supporting procedures are: Findrec, GR, Moms, SMomsAM ,RtoCoe, SSq `): print(``): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(`The main procedures are: GFmom, GGF, MomsAM, MomsFnk, ProveAN, SSqSeq , SSqSeqRec `): print(` `): elif nops([args])=1 and op(1,[args])=Findrec then print(`Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with`): print(`poly coeffs. `): print(`of maximum DEGREE+ORDER<=MaxC. Try:`): print(`Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);`): elif nops([args])=1 and op(1,[args])=GFmom then print(`GFmom(G,x,K): inputs a list G describing a game where a loaded die with k:=nops(G) faces labeled, say`): print(`with POSITIVE dots a1, a2, .., ak with probabilities p1,..., pk, where G is written as a list of k pairs`): print(`[[a1,p1], ..., [ak,pk]], It also inputs variables x. It outputsa list of length K, whose r-th entry is the generating functions in x`): print(` for the r-th (straight) moment of the duration of the game, where `): print(` the game ends the FIRST time the score of the four-year-old is n or above. Try: `): print(` GFmom([[1,1/2],[2,1/2]],x,4); `): print(` GFmom([[1,p],[2,1-p]],x,4); `): print(` GFmom([seq([i,1/6],i=1..6)],x,4); `): elif nops([args])=1 and op(1,[args])=GGF then print(`GGF(G,x,t): inputs a list G describing a game where a loaded die with k:=nops(G) faces labeled, say`): print(`with POSITIVE dots a1, a2, .., ak with probabilities p1,..., pk, where G is written as a list of k pairs`): print(`[[a1,p1], ..., [ak,pk]], It also inputs variables x and t. It outputs the generating function in x`): print(`whose coefficient of x^n is the prob. generating function, in t, of the duration of the game, where`): print(`the game ends the FIRST time the score of the four-year-old is n or above. Try:`): print(`GGF([[1,1/2],[2,1/2]],x,t);`): print(`GGF([[1,p],[2,1-p]],x,t);`): print(`GGF([seq([i,1/6],i=1..6)],x,t);`): elif nops([args])=1 and op(1,[args])=GR then print(`GR(L,x): fits a rational function R(x), to the list L such that R(i)=L[i]`): elif nops([args])=1 and op(1,[args])=ManyGames then print(`ManyGames(G,n1,K,k), plays OneGame(G,n1) K times and returns the statistical information. Try:`): print(`ManyGames([[1,1/2],[2,1/2]],100,1000,4); `): elif nops([args])=1 and op(1,[args])=Moms then print(`Moms(G,n,K): inputs a list G describing a game where a loaded die with k:=nops(G) faces labeled, say`): print(`with POSITIVE dots a1, a2, .., ak with probabilities p1,..., pk, where G is written as a list of k pairs`): print(` [[a1,p1], ..., [ak,pk]], It also inputs variables n. It outputs a list of length K, whose r-th entry is the polynomial `): print(`expression for the r-th (straight) moment of the duration of the game up to exponentially small error, where `): print(`the game ends the FIRST time the score of the four-year-old is n or above, starting at 0. Try:`): print(`The first entry is the expectation, the secod entry is the second moment, etc.`): print(` Moms([[1,1/2],[2,1/2]],n,4); `): print(`Moms([[1,p],[2,1-p]],n,4);`): print(`Moms([seq([i,1/6],i=1..6)],n,4);`): elif nops([args])=1 and op(1,[args])=MomsAM then print(`MomsAM(G,n,K): inputs a list G describing a game where a loaded die with k:=nops(G) faces labeled, say`): print(`with POSITIVE dots a1, a2, .., ak with probabilities p1,..., pk, where G is written as a list of k pairs`): print(` [[a1,p1], ..., [ak,pk]], It also inputs variables n. It outputs a list of length K, whose r-th entry is the polynomial `): print(`expression for the r-th moment ABOUT THE MEAN, of the duration of the game up to exponentially small error, where `): print(`the game starts at 0 and ends the FIRST time the score of the four-year-old is n or above.`): print(`The first entry is the expectation, the secod entry is the second moment, etc. Try:`): print(` MomsAM([[1,1/2],[2,1/2]],n,4); `): print(`MomsAM([[1,p],[2,1-p]],n,4);`): print(`MomsAM([seq([i,1/6],i=1..6)],n,4);`): elif nops([args])=1 and op(1,[args])=MomsAMv then print(`MomsAMv(G,n,K): Story version of MomsAM(G,n,k) (q.v.). Try:`): print(`MomsAMv([[1,1/2],[2,1/2]],n,4);`): print(`MomsAMv([[1,p],[2,1-p]],n,4);`): print(`MomsAMv([seq([i,1/6],i=1..6)],n,4);`): elif nops([args])=1 and op(1,[args])=MomsFnk then print(`MomsFnk(n,L,k,K) : inputs symbols n and k, a positive integer L (for the number of data points to use) and a positive integer K`): print(`outputs explicit polynomial expressions in n, with coefficients that are rational functions in k for`): print(`expressions (mod exponentially small value) for `): print(` (i) average (ii) variance, followed... i-th comement about the mean for i from 3 to K. Try: `): print(` MomsFnk(n,20,k,2); `): elif nops([args])=1 and op(1,[args])=MomsFnkV then print(`MomsFnkV(n,L,k,K): Verbose form MomsFnk(n,L,k,K) (q.v.) . Try:`): print(`MomsFnkV(n,20,k,2);`): elif nops([args])=1 and op(1,[args])=ProveAN then print(`ProveAN(G,K): proves asympotic normality up to order K of the r.v. "duration of game with die G"`): print(`ProveAN([[1,1/2],[2,1/2]],6);`): print(` ProveAN([[1,p],[2,1-p]],8); `): print(`ProveAN([seq([i,1/6],i=1..6)],8);`): elif nops([args])=1 and op(1,[args])=OneGame then print(`OneGame(G,n1): `): print(`Given a prob. distrubution of steps (all forward) and a positive integer n1, simulates one game with the goal of getting to >=n1 . `): print(`OneGame([seq([i,1/6],i=1..6)],100);`): elif nops([args])=1 and op(1,[args])=RtoCoe then print(`RtoCoe(R,x,n): inputs a rational function whose denominator is divisible by a power of (1-x), and outputs the`): print(`polynomial expression that approximates (with exponentially small error) the coefficient of x^n. It is assumed that`): print(`the error is exponentially small. Try:`): print(`RtoCoe(1/(1-x)^2/(2-x),x,n);`): elif nops([args])=1 and op(1,[args])=SMomsAM then print(`SMomsAM(G,n,K): inputs a list G describing a game where a loaded die with k:=nops(G) faces labeled, say`): print(`with POSITIVE dots a1, a2, .., ak with probabilities p1,..., pk, where G is written as a list of k pairs`): print(` [[a1,p1], ..., [ak,pk]], It also inputs variables n. It outputs a list of length K, whose r-th entry is the polynomial `): print(`expression for the r-th SCALED moment ABOUT THE MEAN, of the duration of the game up to exponentially small error, where `): print(`the game starts at 0 and ends the FIRST time the score of the four-year-old is n or above.`): print(`The first entry is the expectation, the secod entry is the second moment, etc. Try:`): print(` SMomsAM([[1,1/2],[2,1/2]],n,4); `): print(`SMomsAM([[1,p],[2,1-p]],n,4);`): print(`SMomsAM([seq([i,1/6],i=1..6)],n,4);`): elif nops([args])=1 and op(1,[args])=SSq then print(`SSq(f,t): Inputs a polynomial f of the variable f outputs the sum of the squares of its coefficients.`): print(`Try:`): print(`SSq(1/2+1/2*t,t);`): elif nops([args])=1 and op(1,[args])=SSqSeq then print(`SSqSeq(G,N): inputs a list G describing a game where a loaded die with k:=nops(G) faces labeled, say`): print(`with POSITIVE dots a1, a2, .., ak with probabilities p1,..., pk, where G is written as a list of k pairs`): print(` [[a1,p1], ..., [ak,pk]], and a positive integer N, outputs the first N entries of the `): print(` sequence: "sum of the squares of the prob. gen. for the duration `): print(` the duration of the game, where `): print(` the game ends the FIRST time the score of the four-year-old is n or above. Try: `): print(` SSqSeq([[1,1/2],[2,1/2]],100); `): print(` SSqSeq([[1,p],[2,1-p]],30); `): print(` SSqSeq([seq([i,1/6],i=1..6)],100); `): elif nops([args])=1 and op(1,[args])=SSqSeqRec then print(`SSqSeqRec(G,n,N,K): inputs a list G describing a game where a loaded die with k:=nops(G) faces labeled, say`): print(`with POSITIVE dots a1, a2, .., ak with probabilities p1,..., pk, where G is written as a list of k pairs`): print(`[[a1,p1], ..., [ak,pk]], symbols n and N (the shift operator in n), and a positive integer K,`): print(`outputs a linear recurrence operator with polynomial coefficients of complexity<=K `): print(` of complexity <=K, annihilating `): print(` sequence: "sum of the squares of the prob. gen. for `): print(` the duration of the game, where `): print(` the game ends the FIRST time the score of the four-year-old is n or above. `): print(` If none exists, then it returns FAIL, and you are welcome to try again with a larger K. `): print(` It also returns the list of initial conditions. The output is a list of length 2, [ope,INI], where INI is a list of numbers`): print(`of the lenght of the order of ope (i.e. the degree in N of ope. Try:`): print(` SSqSeqRec([[1,1/2],[2,1/2]],n,N,8);`): print(`SSqSeqRec([[1,p],[2,1-p]],n,N,8);`): print(` SSqSeqRec([[1,1/3],[2,1/3],[3,1/3]],n,N,16);`): elif nops([args])=1 and op(1,[args])=SSqSeqRecV then print(`SSqSeqRecV(G,n,K1,K2): Verbose form of SSqSeqRec(G,n,N,K1) (q.v.), also gives the exact value of the sequence for n=K2. Try:`): print(`SSqSeqRecV([[1,1/2],[2,1/2]],n,8,2000);`): print(`SSqSeqRecV([[1,p],[2,1-p]],n,8,20);`): print(` SSqSeqRecV([[1,1/3],[2,1/3],[3,1/3]],n,16,2000);`): else print(`There is no ezra for`,args): fi: end: ### From Findrec #findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating #the sequence f of degree DEGREE and order ORDER #For example, try: findrec([seq(i,i=1..10)],0,2,n,N); findrecVerbose:=proc(f,DEGREE,ORDER,n,N) local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a: if (1+DEGREE)*(1+ORDER)+3+ORDER>nops(f) then ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE): fi: ope:=0: var:={}: for i from 0 to ORDER do for j from 0 to DEGREE do ope:=ope+a[i,j]*n^j*N^i: var:=var union {a[i,j]}: od: od: eq:={}: for n0 from 1 to (1+DEGREE)*(1+ORDER)+2 do eq1:=0: for i from 0 to ORDER do eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f): od: eq:= eq union {eq1}: od: var1:=solve(eq,var): kv:={}: for i from 1 to nops(var1) do mu:=op(i,var1): if op(1,mu)=op(2,mu) then kv:= kv union {op(1,mu)}: fi: od: ope:=subs(var1,ope): if ope=0 then RETURN(FAIL): fi: ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}: if nops(ope)>1 then print(`There is some slack, there are `, nops(ope)): print(ope): RETURN(Yafe(ope[1],N)[2]): elif nops(ope)=1 then RETURN(Yafe(ope[1],N)[2]): else RETURN(FAIL): fi: end: #findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating #the sequence f of degree DEGREE and order ORDER #For example, try: findrec([seq(i,i=1..10)],0,2,n,N); findrec:=proc(f,DEGREE,ORDER,n,N) local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a: option remember: #if not findrecEx(f,DEGREE,ORDER,ithprime(20)) then # RETURN(FAIL): #fi: #if not findrecEx(f,DEGREE,ORDER,ithprime(40)) then # RETURN(FAIL): #fi: #if not findrecEx(f,DEGREE,ORDER,ithprime(80)) then # RETURN(FAIL): #fi: if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE): fi: ope:=0: var:={}: for i from 0 to ORDER do for j from 0 to DEGREE do ope:=ope+a[i,j]*n^j*N^i: var:=var union {a[i,j]}: od: od: eq:={}: for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do eq1:=0: for i from 0 to ORDER do eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f): od: eq:= eq union {eq1}: od: var1:=solve(eq,var): kv:={}: for i from 1 to nops(var1) do mu:=op(i,var1): if op(1,mu)=op(2,mu) then kv:= kv union {op(1,mu)}: fi: od: ope:=subs(var1,ope): if ope=0 then RETURN(FAIL): fi: ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}: if nops(ope)>1 then RETURN(Yafe(ope[1],N)[2]): elif nops(ope)=1 then RETURN(Yafe(ope[1],N)[2]): else RETURN(FAIL): fi: end: Yafe:=proc(ope,N) local i,ope1,coe1,L: if ope=0 then RETURN(1,0): fi: ope1:=expand(ope): L:=degree(ope1,N): coe1:=coeff(ope1,N,L): ope1:=normal(ope1/coe1): ope1:=normal(ope1): ope1:= convert( [seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`): factor(coe1),ope1: end: #Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with #poly coffs. #of maximum DEGREE+ORDER<=MaxC #e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2); Findrec:=proc(f,n,N,MaxC) local DEGREE, ORDER,ope,L: for L from 0 to MaxC do for ORDER from 0 to L do DEGREE:=L-ORDER: if (2+DEGREE)*(1+ORDER)+4>=nops(f) then print(`Insufficient data for degree`, DEGREE, `and order `,ORDER): RETURN(FAIL): fi: ope:=findrec([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE,ORDER,n,N): if ope<>FAIL then RETURN(ope): fi: od: od: FAIL: end: #SeqFromRec(ope,n,N,Ini,K): Given the first L-1 #terms of the sequence Ini=[f(1), ..., f(L-1)] #satisfied by the recurrence ope(n,N)f(n)=0 #extends it to the first K values SeqFromRec:=proc(ope,n,N,Ini,K) local ope1,gu,L,n1,j1: ope1:=Yafe(ope,N)[2]: L:=degree(ope1,N): if nops(Ini)<>L then ERROR(`Ini should be of length`, L): fi: ope1:=expand(subs(n=n-L,ope1)/N^L): gu:=Ini: for n1 from nops(Ini)+1 to K do gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)), j1=1..L)]: od: gu: end: ### End From Findrec GR1:=proc(L,x,d1,d2) local eq,var, gu,a,b,i: gu:=(x^d1+add(a[i]*x^i,i=0..d1-1))/add(b[i]*x^i,i=0..d2): var:={seq(a[i],i=0..d1-1),seq(b[i],i=0..d2)}: if nops(L)-nops(var)<=0 then print(`Make L bigger`): RETURN(FAIL): fi: eq:={seq(subs(x=i,gu)-L[i],i=1..nops(L))}: var:=solve(eq,var): if var=NULL then FAIL: else subs(var,gu): fi: end: #GR(L,x): fits a rational function R(x), to the list L such that R(i)=L[i] GR:=proc(L,x) local d1,d2,gu: for d1 from 0 to trunc((nops(L)-2)/2) do for d2 from 0 to trunc((nops(L)-2)/2)-d1 do gu:=GR1(L,x,d1,d2): if gu<>FAIL then RETURN(factor(gu)): fi: od: od: FAIL: end: #GGF(G,x,t): inputs a list G describing a game where a loaded die with k:=nops(G) faces labeled, say #with POSITIVE dots a1, a2, .., ak with probabilities p1,..., pk, where G is written as a list of k pairs #[[a1,p1], ..., [ak,pk]], It also inputs variables x and t. It outputs the generating function in x #whose coefficient of x^n is the prob. generating function, in t, of the duration of the game, where #the game ends the FIRST time the score of the four-year-old is n or above. Try: #GGF([[1,1/2],[2,1/2]],x,t); #GGF([[1,p],[2,1-p]],x,t); #GGF([seq([i,1/6],i=1..6)],x,t); GGF:=proc(G,x,t) local B,T,i,j: B:=1-t*add(G[i][2]*x^G[i][1],i=1..nops(G)): T:=1+ add(G[i][2]*t*add(x^j,j=1..G[i][1]-1),i=1..nops(G)): T/B: end: #GFmom(G,x,K): inputs a list G describing a game where a loaded die with k:=nops(G) faces labeled, say #with POSITIVE dots a1, a2, .., ak with probabilities p1,..., pk, where G is written as a list of k pairs #[[a1,p1], ..., [ak,pk]], It also inputs variables x. It outputsa list of length K, whose r-th entry is the generating functions in x #for the r-th (straight) moment of the duration of the game, where #the game ends the FIRST time the score of the four-year-old is n or above. Try: #GFmom([[1,1/2],[2,1/2]],x,4); #GFmom([[1,p],[2,1-p]],x,4); #GFmom([seq([i,1/6],i=1..6)],x,4); GFmom:=proc(G,x,K) local t,gu,gu1,lu,i: gu:=GGF(G,x,t): gu1:=gu: lu:=[]: for i from 1 to K do gu1:=t*diff(gu1,t): lu:=[op(lu), factor(normal(subs(t=1,gu1)))]: od: lu: end: #RtoCoe(R,x,n): inputs a rational function whose denominator is divisible by a power of (1-x), and outputs the #polynomial expression that approximates (with exponentially small error) the coefficient of x^n. It is assumed that #the error is exponentially small. Try: #RtoCoe(1/(1-x)^2/(2-x),x,n); RtoCoe:=proc(R,x,n) local R1,gu,R1a,k,lu,i: R1:=convert(factor(R),parfrac,x): gu:=0: for i from 1 to nops(R1) do R1a:=op(i,R1): if subs(x=1,denom(R1a))=0 then k:=degree(denom(R1a),x): lu:=normal(R1a*(1-x)^k): gu:= expand(gu+lu*expand(binomial(n+k-1,k-1))): fi: od: gu: end: #Moms(G,n,K): inputs a list G describing a game where a loaded die with k:=nops(G) faces labeled, say #with POSITIVE dots a1, a2, .., ak with probabilities p1,..., pk, where G is written as a list of k pairs #[[a1,p1], ..., [ak,pk]], It also inputs variables n. It outputs a list of length K, whose r-th entry is the polynomial #expression for the r-th (straight) moment of the duration of the game up to exponentially small error, where #the game ends the FIRST time the score of the four-year-old is n or above. Try: #The first entry is the expectation, the secod entry is the second moment, etc. #Moms([[1,1/2],[2,1/2]],n,4); #Moms([[1,p],[2,1-p]],n,4); #Moms([seq([i,1/6],i=1..6)],n,4); Moms:=proc(G,n,K) local x,gu,lu,gu1,i: gu:=GFmom(G,x,K): lu:=[]: for i from 1 to K do gu1:=RtoCoe(gu[i],x,n): if gu1=FAIL then RETURN(lu): else lu:=[op(lu),gu1]: fi: od: lu: end: #MomsAM(G,n,K): inputs a list G describing a game where a loaded die with k:=nops(G) faces labeled, say #with POSITIVE dots a1, a2, .., ak with probabilities p1,..., pk, where G is written as a list of k pairs #[[a1,p1], ..., [ak,pk]], It also inputs variables n. It outputsa list of length K, whose r-th entry is the polynomial #expression for the r-th moment ABOUT THE MEAN of the duration of the game up to exponentially small error, where #the game ends the FIRST time the score of the four-year-old is n or above. Try: #The first entry is the expectation, the secod entry is the second moment, etc. #MomsAM([[1,1/2],[2,1/2]],n,4); #MomsAM([[1,p],[2,1-p]],n,4); #MomsAM([seq([i,1/6],i=1..6)],n,4); MomsAM:=proc(G,n,K) local av,lu,gu,i,j,gu1: option remember: lu:=Moms(G,n,K): if nops(lu)1 then print(`The variance`, gu[2], ` is too big, no concentrarion about the mean`): RETURN(false): fi: for i from 3 to K do if i mod 2=1 then if degree(gu[i],n)>=i/2 then print(`The `,i, `-th moment about the mean`, gu[i], ` is too big`): RETURN(false): fi: else if degree(gu[i],n)<>i/2 then print(`The `,i, `-th moment about the mean`, gu[i], ` is not what is should be`): RETURN(false): fi: lu:=coeff(gu[i],n,i/2)/coeff(gu[2],n,1)^(i/2): if lu<>mul(2*j-1,j=1..i/2) then print(`The scaled `, i, `-th moment`, lu, ` is not `, mul(2*j-1,j=1..i/2)): RETURN(false): fi: fi: od: true: end: #MomsFnk(n,L,k,K): inputs symbols n and k, a positive integer L (for the number of data points to use) and a positive integer K #outputs explicit polynomial expressions in n, with coefficients that are rational functions in k for #expressions (mod exponentially small value) for ##(i) average (ii) variance, followed... i-th comement about the mean for i from 3 to K. Try: #MomsFnk(n,20,k,2); MomsFnk:=proc(n,L,k,K) local i,j,gu,k1,i1,gu0,FA,A,B,vu,d: option remember: FA:=[]: gu:=[seq(MomsAM([seq([i,1/k1],i=1..k1)],n,K) ,k1=1..L)]: gu0:=[seq(gu[i][1],i=1..L)]: A:=GR([seq(coeff(gu0[i],n,1),i=1..L)],k): if A=FAIL then RETURN(FAIL): fi: B:=GR([seq(coeff(gu0[i],n,0),i=1..L)],k): if B=FAIL then RETURN(FAIL): fi: FA:=[A*n+B]: for i from 2 to K do gu0:=[seq(gu[j][i],j=1..L)]: vu:=0: for d from 0 to trunc(i/2) do A:=GR([seq(coeff(gu0[i1],n,d),i1=1..L)],k): if A=FAIL then RETURN(FA): else vu:=vu+A*n^d: fi: od: FA:=[op(FA),vu]: od: FA: end: #SSq(f,t): Inputs a polynomial f of the variable f outputs the sum of the squares of its coefficients. #Try: #SSq(1/2+1/2*t,t); SSq:=proc(f,t) local i: expand(add(coeff(f,t,i)^2,i=ldegree(f,t)..degree(f,t))): end: #SSqSeq(G,N): inputs a list G describing a game where a loaded die with k:=nops(G) faces labeled, say #with POSITIVE dots a1, a2, .., ak with probabilities p1,..., pk, where G is written as a list of k pairs #[[a1,p1], ..., [ak,pk]], and a positive integer N, outputs the first N entries of the #sequence: "sum of the squares of the prob. gen. for the duration #the duration of the game, where #the game ends the FIRST time the score of the four-year-old is n or above. Try: #SSqSeq([[1,1/2],[2,1/2]],100); #SSqSeq([[1,p],[2,1-p]],30); #SSqSeq([seq([i,1/6],i=1..6)],100) SSqSeq:=proc(G,N) local gu,x,t,mu,i: gu:=taylor(GGF(G,x,t),x=0,N+1): mu:=[seq(expand(coeff(gu,x,i)),i=1..N)]: [seq(SSq(mu[i],t),i=1..N)]: end: #MomsAMv(G,n,K): Story version of MomsAM(G,n,k) (q.v.). Try: #MomsAMv([[1,1/2],[2,1/2]],n,4); #MomsAMv([[1,p],[2,1-p]],n,4); #MomsAMv([seq([i,1/6],i=1..6)],n,4); MomsAMv:=proc(G,n,K) local gu,i,a,b,t0: t0:=time(): gu:=MomsAM(G,n,K): if gu=FAIL then RETURN(FAIL): fi: print(``): print(`Statistical Analysis of the Number of Donors a Beggar Needs Before He (or She) can go home`): print(``): print(`By Shalosh B. Ekhad `): print(``): print(`In a certain country there are only`, nops(G), `types of coins, of denominations`, seq(G[i][1],i=1..nops(G)) ): print(`A person passing a beggar throws`): print(``): for i from 1 to nops(G) do print(`A coin of`, G[i][1], `cents , with probability`, G[i][2]): od: print(`Each person's donation is independent of the others.`): print(``): print(`The beggar decides to go home as soon as the amount in his hat exceeds n cents. Let X(n) be the random variable`): print(`"Number of donors until he goes home".`): print(``): print(`Since the expected size of a single donation is`): print(``): a:=add(G[i][1]*G[i][2],i=1..nops(G)): print(a): print(``): print(`The expectation of X(n) is roughly`, n/a): print(``): print(`More precisely, up to exponentially small error, the Expectation of X(n) is`): print(``): print(gu[1]): print(``): print(`and in Maple notation it is`): print(``): lprint(gu[1]): print(``): print(`The variance of X(n), up to exponentially small error is`): print(``): print(gu[2]): print(``): print(`and in Maple notation it is`): print(``): lprint(gu[2]): print(``): b:=coeff(gu[2],n,1): for i from 3 to nops(gu) do print(`The `, i, `-th moment about the mean of of X(n), up to exponentially small error is`): print(``): print(gu[i]): print(``): print(`and in Maple notation it is`): print(``): lprint(gu[i]): print(``): print(`The limit of the scaled `, i , `-th moment as n goes to infinity is`): print(coeff(gu[i],n,trunc((i+1)/2))/b^(i/2)): print(``): od: print(`-----------------------`): print(``): print(`This ends this article that took`, time()-t0 , `seconds to produce. `): print(``): end: #MomsFnkV(n,L,k,K): Verbose form MomsFnk(n,L,k,K) (q.v.) . Try: #MomsFnkV(n,20,k,2); MomsFnkV:=proc(n,L,k,K) local gu,a,i,b,t0: t0:=time(): gu:=MomsFnk(n,L,k,K): print(``): print(`Statistical Analysis of the Number of Spins needed to play a Snakes And Ladders Game (w/o Snakes and Ladders)`): print(``): print(`With a fair k-faced die, terminating as soon as the total meets or exceeds n `): print(``): print(`By Shalosh B. Ekhad `): print(``): print(`Suppose that you keep tossing a fair k-faced die marked with 1 through k dots, and at each step`): print(`you advance according to the number of dots shown by the die. You keep playing until the first time`): print(`your total meets or exceeds n`): print(``): print(`Let X(n) be the random variable, number of rolls until finishing. Of course it depens on k. `): print(``): print(`Since the expected size of a single donation is`): print(``): a:=(k+1)/2: print(a): print(``): print(`The expectation of X(n) is roughly`, n/a): print(``): print(`More precisely, up to exponentially small error, the Expectation of X(n) is`): print(``): print(gu[1]): print(``): print(`and in Maple notation it is`): print(``): lprint(gu[1]): print(``): print(`The variance of X(n), up to exponentially small error is`): print(``): print(gu[2]): print(``): print(`and in Maple notation it is`): print(``): lprint(gu[2]): print(``): b:=coeff(gu[2],n,1): for i from 3 to nops(gu) do print(`The `, i, `-th moment about the mean of of X(n), up to exponentially small error is`): print(``): print(gu[i]): print(``): print(`and in Maple notation it is`): print(``): lprint(gu[i]): print(``): print(`The limit of the scaled `, i , `-th moment as n goes to infinity is`): print(coeff(gu[i],n,trunc((i+1)/2))/b^(i/2)): print(``): od: print(`-----------------------`): print(``): print(`This ends this article that took`, time()-t0 , `seconds to produce. `): print(``): end: #SSqSeqRec(G,n,N,K): inputs a list G describing a ladder game where a loaded die with k:=nops(G) faces labeled, say #with POSITIVE dots a1, a2, .., ak with probabilities p1,..., pk, where G is written as a list of k pairs #[[a1,p1], ..., [ak,pk]], symbols n and N (the shift operator in n), and a positive integer K, #outputs a linear recurrence operator with polynomial coefficients of complexity<=K #of complexity <=K, annihilating #sequence: "sum of the squares of the prob. gen. for #the duration of the game, where #the game ends the FIRST time the score of the four-year-old is n or above. #If none exists, then it returns FAIL, and you are welcome to try again with a larger K. #It also returns the list of initial conditions. #SSqSeqRec([[1,1/2],[2,1/2]],n,N,8); #SSqSeqRec([[1,p],[2,1-p]],n,N,8); SSqSeqRec:=proc(G,n,N,K) local gu,Nu,i,ope: option remember: Nu:=max(seq((2+K-i)*(1+i)+10,i=0..K)): gu:=SSqSeq(G,Nu): ope:=Findrec(gu,n,N,K): if ope=FAIL then FAIL: else [ope,[op(1..degree(ope,N),gu)]]: fi: end: #SSqSeqRecV(G,n,K1,K2): Verbose form of SSqSeqRec(G,n,N,K1) (q.v.). with the value of the K2-th term just for kick. Try: #SSqSeqRecV([[1,1/2],[2,1/2]],n,N,8,1000); SSqSeqRecV:=proc(G,n,K1,K2) local Ope,ope,N,lu,a,d,t0,i,ope1,ka: t0:=time(): Ope:=SSqSeqRec(G,n,N,K1): if Ope=FAIL then print(`Make `, K, `bigger (if you wish), or forget about it. `): RETURN(FAIL): fi: ope:=Ope[1]: d:=degree(ope,N): lu:=SSqSeq(G,d): ope1:=expand(subs(n=n-d,ope)/N^d): ope1:=ope1/coeff(ope1,N,0): ope1:=add(factor(normal(coeff(ope1,N,-i)))*N^(-i),i=0..d): ope1:=1-ope1: print(``): print(`On the Probabilities of the First Player Winning in a Two-Player Pile game where the winner is the first to accumulate n dollars`): print(``): print(`By Shalosh B. Ekhad `): print(``): print(`Consider the following simple 2-Player Pile game`): print(``): print(`Players take turn tossing a die with`, nops(G), `faces marked as follows, with the given probabilities `): print(``): for i from 1 to nops(G) do print(`A face marked with `, G[i][1], `dots , with probability`, G[i][2]): od: print(`and they get the number of chips indicated by the die `): print(`The player to be the FIRST to accumulate n chips is declared the winner. `): print(``): print(`Theorem: The probability of the FIRST player winning is`): print(``): print(1/2+ 1/2*a(n)): print(``): print(`where a(n) is a squence of numbers satisfying the linear recurrence`): print(``): print(a(n)=add(coeff(ope1,N,-i)*a(n-i),i=1..d)): print(``): print(`Subject to the initial conditions`): print(``): print(seq(a(i)=lu[i],i=1..d)): print(``): print(`and in Maple notation`): print(``): lprint(a(n)=add(coeff(ope1,N,-i)*a(n-i),i=1..d)): print(``): print(``): lprint(seq(a(i)=lu[i],i=1..d)): print(``): if type(G[1][2],numeric) then ka:=SeqFromRec(Ope[1],n,N,Ope[2],K2)[K2]: print(`Just for kicks`, a(K2), `equals `, ka): print(``): print(`Hence the probability of the first player winning if the goal is to reach`, K2, `chips is `): print(``): print((1+ka)/2): print(``): print(`that, in floating point is`): print(evalf((1+ka)/2)): print(``): fi: print(`-----------------------`): print(``): print(`This ends this article that took`, time()-t0 , `seconds to produce. `): print(``): end: ##start simulation procedures #OneRoll(L): inputs a list of positive rational mumbers that add up to 1 and outputs # a random outcome from 1 to nops(L) according to this prob. dist. #Try #OneRoll([2/3,1/4,1/12]); OneRoll:=proc(L) local i,M,L1,j,ra: if convert(L,`+`)<>1 then RETURN(FAIL): fi: if min(op(L))<0 then RETURN(FAIL): fi: M:=lcm(seq(denom(L[i]),i=1..nops(L))): L1:=[seq(M*L[i],i=1..nops(L))]: L1:=[seq(add(L1[j],j=1..i),i=1..nops(L))]: L1: M:=L1[nops(L1)]: ra:=rand(1..M)(): for j from 1 to nops(L1) while ra>L1[j] do od: j: end: #OneGame(G,n1): #Given a prob. distrubution of steps (all forward) and a positive integer n1, simulates one game with the goal of getting to >=n1 . #OneGame([seq([i,1/6],i=1..6)],100): OneGame:=proc(G,n1) local n,j,co,lu,L,i1: j:=0: co:=0: L:=[seq(G[i1][2],i1=1..nops(G))]: while j