###################################################################### ##UrnSolitaire.txt: Save this file as UrnSolitaire2.txt # ## To use it, stay in the # ##same directory, get into Maple (by typing: maple ) # ##and then type: read UrnSolitaire2.txt # ##Then follow the instructions given there # ## # ##Written by Doron Zeilberger, Rutgers University , # #DoronZeil at gmail dot com # ###################################################################### #Created: Dec. 28, 2017 read `GuessHolo2.txt`: print(`Created: Dec. 28, 2017`): print(` This is UrnSolitaire.txt `): print(`It is one of the packages that accompany the article `): print(`How Many Rounds Should You Expect in Urn Solitaire?`): print(`by Shalosh B. Ekhad and Doron Zeilberger`): print(` 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): Ezra1:=proc() if args=NULL then print(` The supporting procedures are: Alpha, MakeNice, NuB, LC, OneMoveS, OneRound, OneRoundDet, OneRoundV, Psimple`): print(``): else ezra(args): fi: end: Ezra:=proc() if args=NULL then print(`The main procedures are: A, Asi, Arecs, ArecsS, Paper1, fx, fxS, fxrecs, G, Gdet, Gv, MC, MCdet, MCs, P, Pg, ProveHalf, Psimple `): print(` `): elif nops([args])=1 and op(1,[args])=A then print(`A(a,b): the expected number of rounds in Oakley-Perry Urn Solitaire with a White balls and`): print(`b Black Balls, until the urn only has one color, i.e. the last round is not counted. Try:`): print(`A(5,10);`): elif nops([args])=1 and op(1,[args])=Asi then print(`Asi(a,b): the expected number of rounds in Simple Urn Solitaire with a White balls and`): print(`b Black Balls. The last round (once all the balls are the same color) is not counted. Try:`): print(`Asi(5,10);`): elif nops([args])=1 and op(1,[args])=Arecs then print(`Arecs(m,n,M,N): gueses recurrence operators (that can be proved after the fact) for the function A(m,n), the expected`): print(`number of rounds in Urn Solitaire starting with m White balls and n Black balls. `): print(`MUST have the package GuessHolo2.txt in the same directory. The output is the pair [[ope1,ope2],InitialConditions]`): print(`Try: `): print(` Arecs(m,n,M,N); `): elif nops([args])=1 and op(1,[args])=ArecsS then print(`ArecsS(m,n,M,N): gueses recurrence operators (that can be proved after the fact) for the function Asi(m,n), the expected`): print(`number of rounds in Urn Solitaire starting with m White balls and n Black balls. `): print(`MUST have the package GuessHolo2.txt in the same directory.`): print(`Try: `): print(` Arecs(m,n,M,N); `): elif nops([args])=1 and op(1,[args])=Alpha then print(`Alpha(f,x,N): (borrowed from the Maple package HISTABRUT).`): print(`Given a probability generating function`): print(` f(x) (a polynomial, rational function and possibly beyond) `): print(` returns a list, of length N, whose `): print(` (i) First entry is the average `): print(` (ii): Second entry is the variance `): print(` for i=3...N, the i-th entry is the so-called alpha-coefficients `): print(`that is the i-th moment about the mean divided by the`): print(`variance to the power i/2 (For example, i=3 is the skewness`): print(`and i=4 is the Kurtosis)`): print(`If f(1) is not 1, than it first normalizes it by dividing`): print(`by f(1) (if it is not 0) .`): print(`For example, try:`): print(`Alpha(((1+x)/2)^100,x,4);`): print(``): elif nops([args])=1 and op(1,[args])=fx then print(`fx(a,b,x): the probability generating function, in x, for the r.v. "number of rounds" in the Oakley-Perry Urn Solitaire with a White balls and`): print(`b Black Balls. Try:`): print(`fx(5,10,x);`): elif nops([args])=1 and op(1,[args])=fxS then print(`fxS(a,b,x): the probability generating function, in x, for the r.v. "number of color changes" in Simple Urn Solitaire with a White balls and`): print(`b Black Balls. Try:`): print(`fxS(5,10,x);`): elif nops([args])=1 and op(1,[args])=fxrecs then print(`fxrecs(m,n,M,N,x): gueses recurrence operators (that can be proved after the fact) for the function fx(m,n,x), the expected`): print(`number of rounds in Urn Solitaire starting with m White balls and n Black balls. `): print(`MUST have the package GuessHolo2.txt in the same directory.`): print(`Try:`): print(`fxrecs(m,n,M,N,x);`): elif nops([args])=1 and op(1,[args])=G then print(` G(a,b): Simulating ONE game of Urn Solitaire with initially, a White balls and b Black balls. Terse version`): print(`it outputs the sequence of states after each round. Try:`): print(`G(10,20);`): elif nops([args])=1 and op(1,[args])=Gdet then print(`Gdet(m,n,W,B): simulates Oakley-Perry Urn Solitaire with m White Balls and n Black Balls `): print(`the White balls are called W[1], ..., W[m], the Black balls are called B[1], ..., W[n], `): print(`the output a list of lists of size 4. The first list is the order of the balls, parsed into rounds, that were taken out`): print(`and not replaced, in their order, the second list is the list of balls that were replaced, since`): print(`they finished a round (run of the same balls drawn). This followed by the sets of White Balls and Black Balls at the end`): print(`(one of them is empty, of course`): print(`Try:`): print(`Gdet(3,5,W,B);`): elif nops([args])=1 and op(1,[args])=Gv then print(` Gv(a,b): Verbose version of G(a,b( q.v.). `): print(`Simulating ONE game of Urn Solitaire with initially, a White balls and b Black balls. Terse version`): print(`it outputs the sequence of states after each round. Try:`): print(`Gv(10,20);`): elif nops([args])=1 and op(1,[args])=LC then print(`LC(a,b): in an urn with a White balls and b Black balls, randonly draws a ball. If it is White, the output is 1`): print(`if it is Black, the output is 2. Try:`): print(`LC(30,10);`): elif nops([args])=1 and op(1,[args])=MakeNice then print(`MakeNice(ope,N): making the recurrence operator without denominators. Try:`): print(`MakeNice(N^2/n+N/n^2+1,N);`): elif nops([args])=1 and op(1,[args])=MC then print(`MC(a,b,K,N): K games of G(a,b), (i.e. Oakley-Perry Urn Solitaire) returns the proportion of times White won, `): print(`followed by the list [average number of rounds, variance, ..., N-th moment (about the mean)].`): print(`Try: `): print(`MC(10,10,1000,4);`): elif nops([args])=1 and op(1,[args])=MCdet then print(`MCdet(a,b,K,N): K games of Gdet(a,b), (i.e. Oakley-Perry Urn Solitaire) returns the proportion of times White won, `): print(`followed by the list [average number of rounds, variance, ..., N-th moment (about the mean)].`): print(`Try: `): print(`MCdet(10,10,1000,4);`): elif nops([args])=1 and op(1,[args])=MCs then print(`MCs(a,b,K,N): K games of Gs(a,b), (i.e. SIMPLE Urn Solitaire) returns the proportion of times White won, `): print(`followed by the list [average number of rounds, variance, ..., N-th moment (about the mean)].`): print(`Try: `): print(`MCs(10,10,1000,4);`): elif nops([args])=1 and op(1,[args])=NuB then print(`NuB(L): given a walk, L, finds the number of blocks. Try:`): print(`NuB([[4,3],[4,2],[4,1],[3,1],[2,1],[1,1],[0,1]]);`): elif nops([args])=1 and op(1,[args])=Paper1 then print(`Paper1(m,n,K): Inputs symbols m and n , and pos. integer K, and outputs a paper about the recurrences satisfied by `): print(`the expected number of rounds in the Oakley-Perry Urn Solitaire, also for the diagonal, and gives the value for (K,K)`): print(`MUST have the package GuessHolo2.txt in the same directory.`): print(`Try: `): print(` Paper1(m,n,1000); `): elif nops([args])=1 and op(1,[args])=Paper2 then print(`Paper2(m,n,K): Inputs symbols m and n and a pos. integer K, outputs a paper about the recurrences satisfied by `): print(`the expected number of "rounds" (i.e. chunks) and gives the value for (K,K) in the SIMPLE Urn Solitaire.`): print(`It also gives a recurrence for the prob. gen. function.`): print(`MUST have the package GuessHolo2.txt in the same directory.`): print(`Try: `): print(` Paper2(m,n,1000); `): elif nops([args])=1 and op(1,[args])=OneRound then print(`OneRound(a,b): simulating one round of Oakley-Perry Urn Solitaire, terse version`): print(`It inputs a and b, the initial number of White and Black balls respectively in the urn, and outputs the`): print(`pair [a1,b1], that remain after one round (keep drawing balls until you get a different color, and then you put it back). Try:`): print(`OneRound(10,30);`): elif nops([args])=1 and op(1,[args])=OneRoundDet then print(`OneRoundDet(W,B): simulates One round of Oakley-Perry Urn Solitaire where the Urn contains a`): print(`set W of White Balls and a set B of Black balls.`): print(`It outputs the list of balls of the same color, followed by the first of the opposite color`): print(`that has to be put in the urn. Try:`): print(`OneRoundDet({w1,w2,w3,w4},{b1,b2,b3});`): elif nops([args])=1 and op(1,[args])=OneMoveS then print(`OneMoveS(a,b): simulating one move of SIMPLE Urn Solitaire, terse version`): print(`OneMoveS(10,30);`): elif nops([args])=1 and op(1,[args])=OneRoundV then print(`OneRoundV(a,b): verbose form of OneRound(a,b) (q.v.)`): print(`simulating one round of Oakley-Perry Urn Solitaire, terse version`): print(`It inputs a and b, the initial number of White and Black balls respectively in the urn, and outputs the`): print(`pair [a1,b1], that remain after one round (keep drawing balls until you get a different color, and then you put it back). Try:`): print(`OneRoundV(10,30);`): elif nops([args])=1 and op(1,[args])=P then print(`P(a,b): the probability of the last ball being White in Oakley-Perry Urn solitaire with initiall a White balls and b Black balls.`): print(`Using Dynamical programming. Of course it is always 1/2 (if a>0 and b>0). Try: `): print(`P(6,8);`): elif nops([args])=1 and op(1,[args])=Pg then print(`Pg(L,c);given an urn with nops(L) colors with L[i] balls of color i (1<=i<=nops(L)) and a color c (betweem 1 and nops(L)), find the probability that`): print(`the last ball is of colorc in Oakley-Perry Urn Solitaire. Try`): print(`Pg([3,5],1);`): elif nops([args])=1 and op(1,[args])=Psimple then print(`Psimplle(a,b): the probability of White winning Urn Solitaire without the rule that a different-colored ball is put back into the urn. Try:`): print(`Psimple(4,6);`): elif nops([args])=1 and op(1,[args])=ProveHalf then print(` ProveHalf: proves that P(a,b)=1/2 for a,b>0 . Try:`): print(`ProveHalf(); `): 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 #LC(a,b): in an urn with a White balls and b Black balls, randonly draws a ball. If it is White, the output is 1 #if it is Black, the output is 2. Try: #LC(30,10); LC:=proc(a,b) local ra: ra:=rand(1..a+b)(): if ra<=a then RETURN(1): else RETURN(2): fi: end: #OneRound(a,b): one round of Urn Solitaire OneRound:=proc(a,b) local c,i,j: c:=LC(a,b): if c=1 then for i from 1 to a-1 while LC(a-i,b)=1 do od: RETURN([a-i,b]): elif c=2 then for j from 1 to b-1 while LC(a,b-j)=2 do od: RETURN([a,b-j]): fi: end: #OneRoundV(a,b): one round of Urn Solitaire OneRoundV:=proc(a,b) local c,i,j: c:=LC(a,b): if c=1 then for i from 1 to a-1 while LC(a-i,b)=1 do od: print(`We drew`, i, `consecutive White balls, until we got a Black ball that we put back`): print(`there are currently`, a-i, `White balls and still `, b, `black balls `): RETURN([a-i,b]): elif c=2 then for j from 1 to b-1 while LC(a,b-j)=2 do od: print(`We drew`, j, `consecutive Black balls, until we got a White ball that we put back`): print(`there are currently`, b-j, `Black balls and still `, a, `white balls `): RETURN([a,b-j]): fi: end: #G(a,b): a random game G:=proc(a,b) local gu,mu: mu:=[a,b]: gu:=[mu]: while mu[1]>0 or mu[2]>0 do mu:=OneRound(op(mu)): gu:=[op(gu),mu]: od: [op(1..nops(gu)-1,gu)]: end: #Gv(a,b): a random game, verbose version Gv:=proc(a,b) local gu,mu: mu:=[a,b]: gu:=[mu]: print(`We start with`, a, `White balls and`, b, `Black balls in our urn, and we are playing Urn Solitaire`): print(`invented by B.E. Oakley and R.L. Perry, "A Sampling Process", The Mathematical Gazette, vol. 49, No. 367 (Feb. 1965), pp. 42-44`): print(`and polpularized in Peter Winkler's wonderful book "Mathematical Mind Benders", A.K. Peters, 2007, pp. 76 and 79`): print(``): #while mu<>FAIL and (mu<>FAIL and mu[1]>0) and (mu<>FAIL and mu[2]>0) do while mu[1]>0 or mu[2]>0 do mu:=OneRoundV(op(mu)): if mu=FAIL then print(`Game is over`): RETURN(gu): else gu:=[op(gu),mu]: fi: od: print(``): print(`----------------------------`): print(``): print(`Now the Game is over, and it took`, nops(gu)-2, `color changes`): [op(1..nops(gu)-1,gu)]: end: #MC(a,b,K,N): K games of G(a,b), returns the proportion of times White won, #followed by the list [average number of rounds, variance, ..., N-th moment (about the mean)]. #Try: #MC(10,10,1000,4); MC:=proc(a,b,K,N) local c,gu,mu,x,i: c:=0: mu:=x: for i from 1 to K do gu:=G(a,b): mu:=mu+x^(nops(gu)-1): if gu[-1][1]<>0 and gu[-1][2]=0 then c:=c+1: fi: od: evalf(c/K),evalf(Alpha(mu,x,N)): end: P:=proc(a,b) local k: option remember: if b=0 then RETURN(1): fi: if a=0 then RETURN(0): fi: add(binomial(a,k)/binomial(a+b,k)*b/(a+b-k)*P(a-k,b),k=1..a)+ add(binomial(b,k)/binomial(a+b,k)*a/(a+b-k)*P(a,b-k),k=1..b): end: #ProveHalf: proves that P(a,b)=1/2 for a,b>0 ProveHalf:=proc() local k,a,b: option remember: evalb( normal(convert( (1/2)*sum(binomial(a,k)/binomial(a+b,k)*b/(a+b-k),k=1..a-1)+ (1/2)*sum(binomial(b,k)/binomial(a+b,k)*a/(a+b-k),k=1..b-1)+1/binomial(a+b,a),factorial) -1/2)=0): end: #fx(a,b,x): the probability generating function, in x, for the r.v. "number of rounds" in Urn Solitaire with a White balls and #b Black Balls. Try: #fx(5,10,x); fx:=proc(a,b,x) local k: option remember: if b=0 then RETURN(1): fi: if a=0 then RETURN(1): fi: expand( x*( add(binomial(a,k)/binomial(a+b,k)*b/(a+b-k)*fx(a-k,b,x),k=1..a)+ add(binomial(b,k)/binomial(a+b,k)*a/(a+b-k)*fx(a,b-k,x),k=1..b))): end: #A(a,b): the expected number of rounds (until there is a clear winner) it takes to play Urn Solitaire with a White balls and b Black balls. Try: #A(5,5); A:=proc(a,b) local k: option remember: if b=0 or a=0 then RETURN(0): else 1+add(binomial(a,k)/binomial(a+b,k)*b/(a+b-k)*A(a-k,b),k=1..a)+ add(binomial(b,k)/binomial(a+b,k)*a/(a+b-k)*A(a,b-k),k=1..b): fi: end: #Arecs(m,n,M,N): gueses recurrence operators (that can be proved after the fact) for the function A(m,n), the expected #number of rounds in Urn Solitaire starting with m White balls and n Black balls. #MUST have the package GuessHolo2.txt in the same directory. #Try: #Arecs(m,n,M,N); Arecs:=proc(m,n,M,N) local lu,i,j: lu:=GH2wic((i,j)->A(i,j),m,n,M,N,10,10): if not CheckOpe((i,j)->A(i,j),lu[1][1],M,N,m,n,60) then RETURN(FAIL): elif not CheckOpe((i,j)->A(i,j),lu[1][2],M,N,m,n,60) then RETURN(FAIL): else RETURN(lu): fi: end: #Paper1(m,n,K): verbose form of Arecs(m,n,M,N) (q.v.), also the diagonal recurrence and tells you the terms of A(n,n) from n=K-9 to n=K #Try: #Paper1(m,n,K); Paper1:=proc(m,n,K) local INI,lu,M,N,E,ope1,ope2,k,i,t0,ope,gu,i1,j1: t0:=time(): lu:=Arecs(m,n,M,N): ope1:=MakeNice(lu[1][1],M): ope2:=MakeNice(lu[1][2],N): INI:=lu[2]: print(`Linear Recurrence Equations Satisfied by the Expected Number of Rounds in Oakley-Perry Urn Solitaire`): print(``): print(`By Shalosh B. Ekhad `): print(``): print(`We start with`, m, `White balls and`, n, `Black balls in our urn, and we are playing Urn Solitaire`): print(`invented by B.E. Oakley and R.L. Perry, "A Sampling Process", The Mathematical Gazette, vol. 49, No. 367 (Feb. 1965), pp. 42-44`): print(`and polpularized in Peter Winkler's wonderful book "Mathematical Mind Benders", A.K. Peters, 2007, pp. 76 and 79`): print(``): print(`In this game, you draw a ball at random, and keep drawing balls if they are of the same color as the initial ball`): print(`as soon as you get a ball of the other color, you put it back, and start a new round. Let E(m,n) be the expected`): print(`number of rounds.`): print(``): print(`Theorem 1:`): print(``): print(`E(m,n) satisfies the following linear recurrences in m`): print(``): print(add(coeff(ope1,M,i)*E(m+i,n),i=0..degree(ope1,M))=0): print(``): print(`and in Maple input format`): print(``): lprint(add(coeff(ope1,M,i)*E(m+i,n),i=0..degree(ope1,M))=0): print(``): print(`E(m,n) also satisfies (by symmetry) the following linear recurrences in n`): print(``): print(add(coeff(ope2,N,i)*E(m,n+i),i=0..degree(ope2,N))=0): print(``): print(`and in Maple input format`): print(``): lprint(add(coeff(ope2,N,i)*E(m,n+i),i=0..degree(ope2,N))=0): print(``): print(`subject to the initial conditions`): print(``): print(seq(seq(E(i1,j1)=INI[i1][j1],j1=1..degree(ope2,N)),i1=1..degree(ope1,M))): print(``): print(`Proof `): print(`E(m,n) satisfies the following recurrence`): print(``): print(E(m,n)= 1+Sum(binomial(m,k)/binomial(m+n,k)*n/(m+n-k)*E(m-k,n),k=1..m)+Sum(binomial(n,k)/binomial(m+n,k)*m/(m+n-k)*E(m,n-k),k=1..n)): print(``): print(``): print(`subject to the boundary conditions E(m,0)=0 and E(0,n)=0. `): print(`Verifying this is a routine exercise in the Holonomic Ansatz (e.g. using Christoph Koutschan's package) that we omit`): print(`since we verified it for many cases, and since everything is finite, this is proof enough for us.`): print(``): print(`Corollary:`): print(``): ope:=findrec([seq(A(i,i),i=1..40)],7,3,n,N): ope:=MakeNice(ope,N): print(`The diagonal sequence, E(n,n), satisfies the linear recurrence equation with polynomial coefficies`): print(``): print(add(coeff(ope,N,i)*E(n+i,n+i),i=0..3)=0): print(``): print(``): print(`and in Maple input format`): print(``): lprint(add(coeff(ope,N,i)*E(n+i,n+i),i=0..3)=0): print(``): print(`subject to the initial conditions`, seq(E(i,i)=A(i,i),i=1..3)): print(``): print(`Proof: Routine computations in the holonomic ansatz`): print(``): print(`Using this recurrence, we can fompute many terms`): gu:=SeqFromRec(ope,n,N,[seq(A(i,i),i=1..3)],K): print(``): print(`for example, E(i,i)/i, for i from`, K-9, ` to `,K, ` are ` ): print(evalf([seq(gu[i]/i,i=K-9..K)])): print(``): lu:=ValueFromGH2(lu,m,n,M,N,2*K,3*K): print(`Finally, to illustrate the efficiency of the above-holonomic representation, just for fun, the expected number of rounds when `): print(``): print(`there are`, 2*K, `White balls and`, 3*K, ` Black balls is `): print(``): print(lu): print(``): print(`and in floating-point it is`, evalf(lu)): print(``): print(`Note that using the two dimensional, defining recurrence for E(m,n) would have taken for ever. `): print(`-------------------------------------`): print(``): print(`This ends this article, that took`, time()-t0, `seconds to generate. `): end: #fxrecs(m,n,M,N,x): gueses recurrence operators (that can be proved after the fact) for the function fx(m,n,x), the expected #number of rounds in Urn Solitaire starting with m White balls and n Black balls. #MUST have the package GuessHolo2.txt in the same directory. #Try: #fxrecs(m,n,M,N,x); fxrecs:=proc(m,n,M,N,x) local lu,i,j: lu:=GH2((i,j)->fx(i,j,x),m,n,M,N,15,15): if not CheckOpe((i,j)->fx(i,j,x),lu[1],M,N,m,n,40) then RETURN(FAIL): elif not CheckOpe((i,j)->fx(i,j,x),lu[2],M,N,m,n,40) then RETURN(FAIL): else RETURN(lu): fi: end: #Psimple(a,b): the probability of White winning Urn Solitaire without the rule that a different-colored ball is put back into the urn #try: #Psimple(4,6); Psimple:=proc(a,b) option remember: if a=0 then RETURN(0): elif b=0 then RETURN(1): else a/(a+b)*Psimple(a-1,b)+b/(a+b)*Psimple(a,b-1): fi: end: #Psimple1(a,b): Like Psimple(a,b) but done in chunks Psimple1:=proc(a,b) local k: option remember: if b=0 then RETURN(1): fi: if a=0 then RETURN(0): fi: add(binomial(a,k)/binomial(a+b,k)*b/(a+b-k)*Psimple1(a-k,b-1),k=1..a-1) + add(binomial(b,k)/binomial(a+b,k)*a/(a+b-k)*Psimple1(a-1,b-k),k=1..b-1) +1/binomial(a+b,b): end: #Asi(a,b): The expected number of color changes in simple Urn Solitaire Asi:=proc(a,b) local k: option remember: if b=0 then RETURN(0): fi: if a=0 then RETURN(0): fi: 1+add(binomial(a,k)/binomial(a+b,k)*b/(a+b-k)*Asi(a-k,b-1),k=1..a-1)+add(binomial(b,k)/binomial(a+b,k)*a/(a+b-k)*Asi(a-1,b-k),k=1..b-1): end: #fxS(a,b,x): The prob. gen. function for the r.v. "number of color changes" in simple Urn Solitaire fxS:=proc(a,b,x) local k: option remember: if b=0 then RETURN(1): fi: if a=0 then RETURN(1): fi: expand(x*(add(binomial(a,k)/binomial(a+b,k)*b/(a+b-k)*fxS(a-k,b-1,x),k=1..a-1)+add(binomial(b,k)/binomial(a+b,k)*a/(a+b-k)*fxS(a-1,b-k,x),k=1..b-1)) +x*(1/binomial(a+b,b)+1/binomial(a+b,a))): end: #OneRoundS(a,b): one round of Simple Urn Solitaire OneRoundS:=proc(a,b) local c,i,j: c:=LC(a,b): if c=1 then for i from 1 to a-1 while LC(a-i,b)=1 do od: RETURN([a-i,b-1]): elif c=2 then for j from 1 to b-1 while LC(a,b-j)=2 do od: RETURN([a-1,b-j]): fi: end: #OneMoveS(a,b): one move in the Simple Urn Solitaire. Try: OneMoveS(10,7); OneMoveS:=proc(a,b) local c: c:=LC(a,b): if c=1 then RETURN([a-1,b]): else RETURN([a,b-1]): fi: end: #Gs(a,b): a random game in Simple Urn Solitaire Gs:=proc(a,b) local gu,mu: mu:=[a,b]: gu:=[mu]: while mu[1]>0 or mu[2]>0 do mu:=OneMoveS(op(mu)): gu:=[op(gu),mu]: od: gu: end: #NuLs(v):The number of Ls in a list v NuLs:=proc(v) local i,v1: if v=[] or nops(convert(v,set))=1 then RETURN(0): fi: for i from 1 to nops(v) while v[i]=v[1] do od: NuLs([op(i+1..nops(v),v)])+1: end: #NuB(L): given a walk, L, finds the number of runs. Try: #NuB([[4,3],[4,2],[4,1],[3,1],[2,1],[1,1],[0,1]]); NuB:=proc(L) local i,lu,c: lu:=[]: for i from 1 to nops(L)-1 do if L[i][1]=L[i+1][1] then lu:=[op(lu),1]: else lu:=[op(lu),2]: fi: od: NuLs(lu): end: #MCs(a,b,K,N): K games of Gs(a,b), returns the proportion of times White won, #followed by the list [average number of rounds, variance, ..., N-th moment (about the mean)]. #Try: #MCs(10,10,1000,4); MCs:=proc(a,b,K,N) local c,gu,mu,x,i: c:=0: mu:=x: for i from 1 to K do gu:=Gs(a,b): mu:=mu+x^(NuB(gu)): if gu[nops(gu)-1][1]=1 then c:=c+1: fi: od: evalf(c/K),evalf(Alpha(mu,x,N)): end: #ArecsS(m,n,M,N): gueses recurrence operators (that can be proved after the fact) for the function A(m,n), the expected #number of rounds in Urn Solitaire starting with m White balls and n Black balls. #MUST have the package GuessHolo2.txt in the same directory. #Try: #ArecsS(m,n,M,N); ArecsS:=proc(m,n,M,N) local lu,i,j: lu:=GH2((i,j)->Asi(i,j),m,n,M,N,10,10): if not CheckOpe((i,j)->A(i,j),lu[1],M,N,m,n,60) then RETURN(FAIL): elif not CheckOpe((i,j)->A(i,j),lu[2],M,N,m,n,60) then RETURN(FAIL): else RETURN(lu): fi: end: #Paper2(m,n,K): Paper abot the simple Urn Solitaire #Try: #Paper2(m,n,K); Paper2:=proc(m,n,K) local N,E,i,t0,ope1,ope2,ope3,gu,F,x,V,mu: t0:=time(): print(`Linear Recurrence Equations Satisfied by the Expected Number of Rounds in Simple Urn Solitaire`): print(``): print(`By Shalosh B. Ekhad `): print(``): print(`We start with`, m, `White balls and`, n, `Black balls in our urn, and we are playing Simple Urn Solitaire`): print(``): print(`In this game, you draw a ball at random, and keep drawing balls if they are of the same color as the initial ball`): print(`as soon as you get a ball of the other color, the round ends, and a new round starts`): print(`Let E(m,n) be the expected number of rounds.`): print(``): print(`Theorem 1:`): print(``): ope1:=findrec([seq(Asi(i,i),i=1..40)],3,2,n,N): ope1:=MakeNice(ope1,N): print(`The diagonal sequence, E(n,n), satisfies the linear recurrence equation with polynomial coefficies`): print(``): print(add(coeff(ope1,N,i)*E(n+i,n+i),i=0..3)=0): print(``): print(``): print(`and in Maple input format`): print(``): lprint(add(coeff(ope1,N,i)*E(n+i,n+i),i=0..3)=0): print(``): print(`subject to the initial conditions`, seq(E(i,i)=Asi(i,i),i=1..3)): print(``): print(`Proof: Routine computations in the holonomic ansatz`): print(``): print(`Using this recurrence, we can compute many terms`): gu:=SeqFromRec(ope1,n,N,[seq(Asi(i,i),i=1..2)],K): print(``): print(`for example, E(i,i)/i, for i from`, K-9, ` to `,K, ` are ` ): print(evalf([seq(gu[i]/i,i=K-9..K)])): print(``): print(`We also have:`): print(``): print(`Theorem 2: Let F(n,n,x) be the probability generating function of the r.v. "number of rounds" in Simple Urn Solitaire`): print(``): gu:=[seq(fxS(i,i,x),i=1..60)]: ope2:=findrec(gu,2,3,n,N): ope2:=MakeNice(ope2,N): print(``): print(` It satisfies the linear recurrence equation with polynomial coefficies`): print(``): print(add(coeff(ope2,N,i)*F(n+i,n+i,x),i=0..degree(ope2,N))=0): print(``): print(`and in Maple input format`): print(``): lprint(add(coeff(ope2,N,i)*F(n+i,n+i,x),i=0..degree(ope2,N))=0): print(``): print(`Subject to the initial conditions`, seq(F(i,i,x)=fxS(i,i,x),i=1..3)): mu:= [seq( subs(x=1,diff(diff(x*gu[i],x),x))-Asi(i,i)^2,i=1..nops(gu))]: ope3:=Findrec(mu,n,N,10): ope3:=MakeNice(ope3,N): if ope3<>FAIL then print(`Theorem 3: The sequence of variances V(n,n) satisfies the linear recurrence equation`): print(``): print(add(coeff(ope3,N,i)*V(n+i,n+i),i=0..degree(ope3,N))=0): print(``): print(`and in Maple input format`): print(``): lprint(add(coeff(ope3,N,i)*V(n+i,n+i),i=0..degree(ope3,N))=0): print(``): print(`Subject to the initial conditions`, seq(V(i,i)=mu[i],i=1..degree(ope3,N))): fi: print(`Using this recurrence, we can compute many terms`): gu:=SeqFromRec(ope3,n,N,[op(1..degree(ope3,N),mu)],K): print(``): print(`for example, V(i,i)/i, for i from`, K-9, ` to `,K, ` are ` ): print(evalf([seq(gu[i]/i,i=K-9..K)])): print(`-------------------------------------`): print(``): print(`This ends this article, that took`, time()-t0, `seconds to generate. `): [ [ope1, [Asi(1,1),Asi(2,2)]] , [ope2,[seq(fxS(i,i,x),i=1..3)] ], [ope3, [op(1..degree(ope3,N), mu)]] ] : end: #OneRoundDet(W,B): simulates One round of Oakley-Perry Urn Solitaire where the Urn contains a #set W of White Balls and a set B of Black balls. #It outputs the list of balls of the same color, followed by the first of the opposite color #that has to be put in the urn. Try: #OneRoundDet({w1,w2,w3,w4},{b1,b2,b3}); OneRoundDet:=proc(W,B) local Ur,kadur,kadur1, gu1,lavan: Ur:=W union B: kadur:=Ur[rand(1..nops(Ur))()]: gu1:=[kadur]: Ur:=Ur minus {kadur}: if member(kadur,W) then lavan:=1: else lavan:=0: fi: if lavan=1 then while Ur<>{} do kadur1:=Ur[rand(1..nops(Ur))()]: if member(kadur1,B) then RETURN(gu1,kadur1): else Ur:=Ur minus {kadur1}: gu1:=[op(gu1),kadur1]: fi: od: fi: if lavan=0 then while Ur<>{} do kadur1:=Ur[rand(1..nops(Ur))()]: if member(kadur1,W) then RETURN(gu1,kadur1): else gu1:=[op(gu1),kadur1]: Ur:=Ur minus {kadur1}: fi: od: fi: end: #Gdet(m,n,W,B): simulates Oakley-Perry Urn Solitaire with m White Balls and n Black Balls #the White balls are called W[1], ..., W[m], the Black balls are called B[1], ..., W[n], #the output a list of lists of size 4. The first list is the order of the balls, parsed into rounds, that were taken out #and not replaced, in their order, the second list is the list of balls that were replaced, since #they finished a round (run of the same balls drawn). This followed by the sets of White Balls and Black Balls at the end #(one of them is empty, of course #Try: #Gdet(3,5,W,B); Gdet:=proc(m,n,W,B) local WhB,BlB,i,gu1,gu2,mu: WhB:={seq(W[i],i=1..m)}: BlB:={seq(B[i],i=1..n)}: gu1:=[]: gu2:=[]: while WhB<>{} and BlB<>{} do mu:=OneRoundDet(WhB,BlB): WhB:=WhB minus {op(mu[1])}: BlB:=BlB minus {op(mu[1])}: gu1:=[op(gu1),mu[1]]: gu2:=[op(gu2),mu[2]]: od: [gu1,gu2,WhB,BlB]: end: #MCdet(a,b,K,N): K games of Gdet(a,b), returns the proportion of times White won, #followed by the list [average number of rounds, variance, ..., N-th moment (about the mean)]. #Try: #MCdet(10,10,1000,4); MCdet:=proc(a,b,K,N) local c,gu,mu,x,i,W,B: c:=0: mu:=x: for i from 1 to K do gu:=Gdet(a,b,W,B): mu:=mu+x^(nops(gu[1])): if gu[3]<>{} then c:=c+1: fi: od: evalf(c/K),evalf(Alpha(mu,x,N)): end: #Gdet1(m,n,W,B): simulates Oakley-Perry Urn Solitaire with m White Balls and n Black Balls #the White balls are called W[1], ..., W[m], the Black balls are called B[1], ..., W[n], #the output a list of lists of size 4. The first list is the order of the balls, parsed into rounds, that were taken out #and not replaced, in their order, the second list is the list of balls that were replaced, since #they finished a round (run of the same balls drawn). This followed by the sets of White Balls and Black Balls at the end #(one of them is empty, of course #Try: #Gdet1(3,5,W,B); Gdet1:=proc(m,n,W,B) local WhB,BlB,i,gu,mu,pi: WhB:={seq(W[i],i=1..m)}: BlB:={seq(B[i],i=1..n)}: gu:=[]: while WhB<>{} and BlB<>{} do mu:=OneRoundDet(WhB,BlB): WhB:=WhB minus {op(mu[1])}: BlB:=BlB minus {op(mu[1])}: gu:=[op(gu),op(mu[1]),mu[2]]: od: if WhB<>{} then pi:=randperm(nops(WhB)): gu:=[op(gu), seq(WhB[pi[i]],i=1..nops(pi))]: elif WhB<>{} then pi:=randperm(nops(WhB)): gu:=[op(gu), seq(WhB[pi[i]],i=1..nops(pi))]: fi: gu: end: #MakeNice(ope,N): making the recurrence operator without denominators MakeNice:=proc(ope,N) local ope1,i: if ope=FAIL then RETURN(FAIL): fi: ope1:=numer(normal(ope)): add(factor(coeff(ope1,N,i))*N^i,i=0..degree(ope1,N)): end: #Pg(L,c);given an urn with nops(L) colors with L[i] balls of color i (1<=i<=nops(L)) and a color c (betweem 1 and nops(L)), find the probability that #the last ball is of colorc in Oakley-Perry Urn Solitaire. Try #Pg([3,5],1); Pg:=proc(L,c) local k,i,gu,x,lu: option remember: lu:=add(x[L[i]],i=1..nops(L)): if coeff(lu,x[0],1)=nops(L)-1 then for i from 1 to nops(L) while L[i]=0 do od: if i=c then RETURN(1): else RETURN(0): fi: fi: gu:=0: for i from 1 to nops(L) do gu:=gu+add(binomial(L[i],k)/binomial(convert(L,`+`),k)*(convert(L,`+`)-L[i])/(convert(L,`+`)-k)*Pg([op(1..i-1,L),L[i]-k,op(i+1..nops(L),L)],c), k=1..L[i]): od: gu: end: