A story of balls, randomness and PDEs

TL;DR

This paper introduces a novel PDE-based method to solve complex probability problems involving balls and random removals, providing explicit formulas and demonstrating the approach's versatility over traditional inductive solutions.

Contribution

It presents a systematic PDE approach to derive probability distributions in a complex ball removal problem, offering a new algorithmic and analytical framework.

Findings

Derived explicit probability formulas for remaining balls.

Solved the problem using a PDE approach with boundary conditions.

Validated the method by reproducing known results.

Abstract

Several differential equations usually appearing in mathematical physics are solved through a power series expansion, which reduces in solving difference equations. In this paper a probability problem is presented whose solution follows a completely reversed but systematic approach. Hence, this work is about illustrating how complex probability problems could be tackled with the more powerful techniques of a better studied and well understood field, that of differential equations. The problem is defined as follows: Inside a box containing r red and w white balls random removals occur. The balls are removed successively according to the three following rules. Rule I: If a white ball is chosen it is immediately discarded. If a red ball is chosen, it is placed back into the box and a new ball is randomly chosen. The second ball is then removed irrespective of the color. Rule II: Once one ball is removed, the game continues from Rule I. Rule III: The game ends once all the red balls are removed. The question posed is the determination of the probability that k white balls remain where k = 0, 1, 2, ..., w. Ending the game once all the white balls are removed, a second question is the determination of the probability that k red balls remain where k = 0, 1, 2, ..., r. While inductive solutions are possible, the current approach demonstrates a different and algorithmic route. In particular, the law of total probability yields a recursion that is transformed into a linear inhomogeneous 2D PDE, with suitable boundary conditions. The PDE solutions, which are found analytically, provide the generating functionals of the required probabilities as a function of r, w and k. Using the functionals, the probability formulas for any r, w and k are finally obtained in a closed form. Reproducing existing results of the literature this method is quite generic and adaptable to a large class of problems.

Equations

$$\boxed{\begin{array}{rl}{p}_{III}^{(k)}(r,w)& =k!\frac{r!(r+w+1)}{(r+k+1)!}\frac{rw!(r+w-k-1)!}{(r+w)!(w-k)!}\text{mbox},RuleIII,\\ p_{IV}^{(k)}(r,w)& =(2k+1)\frac{r!(r+w+k)!}{(r+k+1)!(w+r)!}\frac{wr!(r+w-k-1)!}{(r+w)!(r-k)!}\text{mbox},RuleIV.\end{array}}$$

$$\boxed{\begin{array}{r}\sum _{k=0}^w{p}_{III}^{(k)}(r,w)=1,\sum _{k=0}^r{p}_{IV}^{(k)}(r,w)=1,\end{array}}$$

$$\displaystyle\delta_{j}^{i}=\mbox{$1$ if $i=j$ $(i,\,j\,=\,\mathbb{Z})$ and $0$ otherwise}.$$

$$p_{III}^{(k)}(r,w)\equiv$$

$$\equiv$$

$$p_{III}^{(k)}(0,w)\equiv$$

$$p_{III}^{(0)}(0,0)\equiv$$

$$p_{III}^{(k)}(r\ge 0,0)\equiv$$

$$p_{III}^{(k)}(r,w<k)\equiv$$

$$p_{IV}^{(k)}(r,w)\equiv$$

$$\equiv$$

$$p_{IV}^{(k)}(r,0)\equiv$$

$$p_{IV}^{(0)}(0,0)\equiv$$

$$p_{IV}^{(k)}(0,w\ge 0)\equiv$$

$$p_{IV}^{(k)}(r<k,w)\equiv$$

$$p^{(k)}(r,w)=\frac{r}{r+w}{p}^{(k)}(r-1,w)+\frac{w}{r+w}{p}^{(k)}(r,w-1).$$

$$p^{(k)}(r,w)=\frac{r!w!}{(r+w)!}{f}^{(k)}(r,w)$$

$$f^{(k)}(r,w)=f^{(k)}(r-1,w)+f^{(k)}(r,w-1).$$

$$\sum _{r=0}^{\infty }\sum _{w=0}^{\infty }{f}_{r,w}^{(k)}{z}^{-r}{u}^{-w}-\sum _{w=0}^{\infty }{f}_{0,w}^{(k)}{u}^{-w}-\sum _{r=0}^{\infty }{f}_{r,0}^{(k)}{z}^{-r}+f_{0,0}=$$

$$\frac{1}{z}\left(\sum _{r=0}^{\infty }\sum _{w=0}^{\infty }{f}_{r,w}^{(k)}{z}^{-r}{u}^{-w}-\sum _{r=0}^{\infty }{f}_{r,0}^{(k)}{z}^{-r}\right)+\frac{1}{u}\left(\sum _{r=0}^{\infty }\sum _{w=0}^{\infty }{f}_{r,w}^{(k)}{z}^{-r}{u}^{-w}-\sum _{w=0}^{\infty }{f}_{0,w}^{(k)}{u}^{-w}\right)$$

$$Y^{(k)}(z,u)=\sum _{r=0}^{\infty }\sum _{w=0}^{\infty }{f}_{r,w}^{(k)}{z}^{-r}{u}^{-w}$$

$$Y^{(k)}(z,u)=\frac{1}{zu-z-u}\left(u(z-1)\sum _{r=0}^{\infty }{f}_{r,0}^{(k)}{z}^{-r}+z(u-1)\sum _{w=0}^{\infty }{f}_{0,w}^{(k)}{u}^{-w}-zu{f}_{0,0}^{(k)}\right).$$

$$\displaystyle Y^{(k>0)}(z,u)=\frac{z(u-1)}{zu-z-u}\times u^{-k}\mbox{ (because $f^{(k>0)}_{r,0}=0\,\,\forall\,r$ and $f^{(k)}_{0,w}=\delta^{k}_{w}$).}$$

$$\frac{z(u-1)}{zu-z-u}$$

$$=\sum _{r=0}^{\infty }\sum _{w=0}^{\infty }\frac{(r+w)!}{r!w!}\frac{1}{z^r}\frac{1}{u^w}-\sum _{r=0}^{\infty }\sum _{w=1}^{\infty }\frac{(r+w-1)!}{r!(w-1)!}\frac{1}{z^r}\frac{1}{u^w}$$

$$=\sum _{r=0}^{\infty }\sum _{w=0}^{\infty }\left(\frac{(r+w)!}{r!w!}-\frac{(r+w-1)!}{r!(w-1)!}\right)\frac{1}{z^r}\frac{1}{u^w}\left(\frac{\text{mbox}(i)assuming(-1)!=\Gamma (0)=\infty ,}{\text{mbox}(ii)considering(\text{ref}pa).}\right)$$

$$=\sum _{r=0}^{\infty }\sum _{w=0}^{\infty }\frac{r(r+w-1)!}{r!w!}\frac{1}{z^r}\frac{1}{u^w}$$

$$Y^{(k)}(z,u)=\sum _{r=0}^{\infty }\sum _{w=k}^{\infty }\frac{r(r+w-k-1)!}{r!(w-k)!}\frac{1}{z^r}\frac{1}{u^w}=\sum _{r=0}^{\infty }\sum _{w=0}^{\infty }\frac{r(r+w-k-1)!}{r!(w-k)!}\frac{1}{z^r}\frac{1}{u^w}$$

$$p^{(k)}(r,w)=\frac{r(r+w-k-1)!w!}{(r+w)!(w-k)!},\forall r,w\ge 0,k>0.$$

$$\displaystyle Y^{(k=0)}(z,u)=\frac{z(u-1)}{zu-z-u}\mbox{ (because $f^{(k=0)}_{r,0}=1\,\,\forall\,r$ and $f^{(k=0)}_{0,w}=\delta^{0}_{w}$).}$$

$$\boxed{p^{(k)}(r,w)=\frac{r(r+w-k-1)!w!}{(r+w)!(w-k)!}},\forall k,r,w\ge 0,k\le w,$$

$$\sum _{k=0}^w{p}^{(k)}(r,w)=1.$$

$$A=$$

$$B=$$

$$r$$

$$p^{(k)}(r,w)$$

$$=p(B)\frac{1}{1+k}=\frac{\left(\genfrac{}{}{0ex}{}{r}{r-1}\right)\left(\genfrac{}{}{0ex}{}{w}{w-k}\right)}{\left(\genfrac{}{}{0ex}{}{r+w}{r+w-k-1}\right)}\frac{1}{1+k}=\frac{r(r+w-k-1)!w!}{(r+w)!(w-k)!},$$

$$G(r,w;z)\equiv \sum _{w=0}^k{p}^{(k)}(r,w)z^k=\frac{r}{r+w}{}_2{F}_1\left(1,-w,1-r-w;z\right).$$

$$p^{(k)}(r,w)=\frac{r^2}{(r+w{)}^2}{p}^{(k)}(r-1,w)+\frac{w^2+2rw}{(r+w{)}^2}{p}^{(k)}(r,w-1)$$

Full text

the three following rules. Rule I: If a white ball is chosen it is immediately discarded from the box. If a red ball is chosen, it is placed back into the box and a new ball is randomly chosen. The second ball is then removed irrespective of the color. Rule II: Once one ball is removed, the game continues from Rule I. Rule III: The game ends once all the red balls are removed. The first question posed is the determination of the probability that $k$ white balls remain where $k=0,\,1,\,2,\,...,\,w$. Ending the game once all the white balls are removed, a second question is the determination of the probability that k red balls remain where $k=0,\,1,\,2,\,...,\,r$. While inductive solutions are possible, the current approach demonstrates a different and algorithmic route. In particular, the law of total probability yields a recursive equation that is transformed into a linear PDE in two dimensions with inhomogeneous source terms, and suitable boundary conditions that depend on $k$. The PDE solutions, which are found analytically, provide up to a known rescaling factor, the (two) generating functionals of the required probabilities as a function of $r$, $w$ and $k$. Using then the derived functionals, the required probability formulas for any $r$, $w$ and $k$ are finally obtained in a closed form; the probability distributions turn out to be linear combinations of hypergeometric functions of type ${}_{3}F_{2}$. Reproducing existing results of the academic literature, which are special but less involved cases of the current completely generic solution, this method is quite generic and adaptable to a large class of problems.

1. 1 Introduction

2. 2 A toy model

3. 2.1 Solution through $z$-transformations

4. 2.2 A combinatorial derivation

5. 3 Setting up the problem: From the law of total probability to the required PDE

6. 3.1 The PDE of the generating functional with rescaling

7. 3.2 The PDE of the generating functional without rescaling

8. 4 Reproducing known results from the literature

9. 5 Deriving the generic generating functional for any number of remaining white balls $k$, Rule III

10. 6 The general probability formula for Rule III

11. 6.1 The case $k=0$

12. 6.2 The case $k=1$

13. 6.3 The case $k=2$

14. 6.4 The general solution

15. 6.4.1 Maximal probability

16. 6.4.2 Limiting cases

17. 7 The general probability formula for Rule IV

18. 7.1 The PDE, the boundary conditions and the solution, Rule IV

19. 7.2 The probabilities for the first few cases, Rule IV

20. 7.2.1 The case $k=0$

21. 7.2.2 The case $k=1$

22. 7.2.3 The case $k=2$

23. 7.3 The general solution

24. 7.3.1 Maximal probability

25. 7.3.2 Limiting cases

26. 8 Discussion

27. A Gamma function expansions near poles

28. B Alternative boundary conditions, Problem 1, Rule III

29. C A second derivation of (44)

30. D Finite summation identities involving a product of binomials and a rational function

31. E Proving the normalization of the probability formula, Rule III

32. F Solution to Problem 1, Rule III: a combinatorial-inductive derivation

33. G Problem 4: A challenging problem from the literature

34. H Python simulation code and some results

1 Introduction

During interviews for quantitative researchers in investment banks, hedge funds, asset management firms, fin-tech companies etc, the candidates are often requested to solve various probability problems. The current problem is motivated by such interviews and is a more complicated version of these kind of questions.

Problem 1

Inside a box containing $r$ red and $w$ white balls random removals occur. The balls are removed one by one according to the three following rules. Rule I: If a white ball is chosen it is immediately discarded from the box. If a red ball is chosen, it is placed back into the box and a new ball is randomly chosen. The second ball is then removed irrespective of the color. Rule II: Once a ball is removed, the game continues from Rule I. Rule III: The game ends once all the red balls are removed. The question posed is the determination of the probability that $k$ white balls remain where $k=0,\,1,\,2,\,...,\,w$. Changing rule Rule III into Rule IV: The game ends once all the white balls are removed. The second question posed is the determination of the probability that $k$ red balls remain where $k=0,\,1,\,2,\,...,\,r$.

This paper in not only about providing the solution to Problem 1, as, according to appendix F, a combinatorial-inductive solution111We are particularly grateful to D. Christofides for providing us with the solution through such an approach. is (also) possible. This paper is rather about illustrating how difficult problems involving difference equations usually appearing in probability questions, could be tackled systematically with the more powerful techniques of a better studied and well understood field, that of differential equations. As it is known, differential equations appearing in mathematical physics, such as Bessel functions, Hypergeometric functions, Hermite polynomials, to name a few, are solved perturbatively using the Frobenius method. During the course of the solution, recursive relations that connect higher order coefficients with lower order ones appear; that is is difference equations.

In this problem, we work in the reversed order following an algorithmic recipe explained below and summarized in five steps in section 8. That is starting from a probability problem, which is reduced to a linear difference equation in two variables, we deduce an inhomogeneous linear partial differential equation (PDE) with suitable boundary conditions that we solve analytically. The solution, up to a known constant, provides the probability generating functional of the problem at hand. Expanding then the functional into a Laurent series, we explicitly obtain the probability solutions in closed form, which, as we will compute through our step by step systematic approach, are given by equations (75) and (91) for the Rule III and the Rule IV respectively, and are reproduced below. In particular, the probability $p^{(k)}(r,w)$, where $k$ is the number of the white or red remaining balls, and $r$ and $w$ are the red and the white initial balls inside the box respectively, is given by

[TABLE]

Satisfying the normalization conditions

[TABLE]

the two probability formulas consist of the main result of this paper.

The proposed approach combines several beautiful branches of mathematics: probability theory, differential equations, special functions and peripherally, through analytic continuation, residues and Laurent expansions of functions, complex analysis. Most of the steps between equations are either done explicitly or explained in detail. Therefore, being self contained, this paper can be served as an introductory set of lecture notes for advanced undergraduate and graduate students with mathematical background that are interested in applications of the aforementioned fields of mathematics.

For a quick but yet a detailed exposition of the reader to the main ideas of this work, it suffices to restrict to sections 2.1, 3.1 and 4.

This work is organized as follows.

In section 2, using $z$-transformations in two dimensions, we solve an easier problem, Problem 2, that involves algebraic but not differential equations. This serves as an introduction to Problem 1 and to the subsequent sections, and clarifies the formalism and the notation.

In section 3 we return to the original problem where using the law of total probability and $z$-transformations, we eventually obtain the required PDE (of Problem 1). The boundary conditions are discussed in the subsequent sections.

Section 4 deals with Problem 3, which is a special case of Problem 1, and whose result is known in the community [1]. Hence, we cross-check our method against known cases.

Returning to Problem 1, Rule III, section 5 discusses the boundary conditions and derives the generic solution of the PDE for any remaining white balls $k=0,\,1,\,2\,...,\,w$.

Using the solutions of the PDE from section 5, the first part of section 6 derives the explicit probability formulas for the cases $k=0,\,1,\,2$ (Rule III). Given these solutions and by observing the pattern (there exists an even more systematic approach), the remaining section derives the generic probabilities $p_{III}^{(k)}(r,w)$ as a function of the initial balls $r$ and $w$ and of the remaining white balls $k$, which consists the one of the two main results of this paper. The probabilities $p_{III}^{(k)}(r,w)$ are subsequently being investigated.

Armed with the experience from earlier sections and in particular from section 6, section 7 outlines and eventually provides the second main result of this paper, namely the formula of $p_{IV}^{(k)}(r,w)$ for Problem 1, Rule IV. This section is much shorter than sections 5 and 6.

The last section, section 8, concludes by discussing the advantages of our approach.

Appendices A, C, D and E are reserved for longer calculations and for proving several useful identities. Appendix B discusses alternative boundary conditions for Problem 1, Rule III. Appendix F provides the solution to Problem 1, Rule III through a combinatorial-inducive approach. On the other hand, appendix G documents a different and hard question posed in [2], which serves as another reference for these kind of problems. Lastly, appendix H provides a python code that simulates Problem 1 and provides the results for a few cases contrasting them against the two analytically obtained formulas.

Notation:

We use lower case latin letters for the remaining balls inside the box. In particular, we use the letter $r$ for the red balls and the letter $w$ for the white balls that the box contains at the beginning of the games. The letter $k$ is used for the remaining white or red balls (for Rule III and for Rule IV respectively) after the game ends. We use other lower case latin letters except from $k$, $r$ and $w$ in order to sum over other indices that may appear along the way. Finally, during the $z-$transformations, the letters $z$ and $u$ are associated with the red and the white balls respectively.

We define the Kronecker delta by

[TABLE]

The following definitions apply for Problem 1, Rule III (and for Problem 2).

[TABLE]

Equation (2c) is a special case of (2d) but we make it explicit. What it should be emphasized is that equations (2b)-(2e) should be consistent with equation (2a); this can be checked explicitly once the formula for $p_{III}^{(k)}(r,w)$ is obtained. In particular, as we will see in section 5, either equation (2b) or equation (2d), under particular modifications that we will derive in what follows, should be used as the boundary conditions for the continuum version of the problem.

The motivation behind the specification of the ordering of the limits in equation (2a) is because we would like to write compact expressions for the final probability formulas by giving a meaning to subtle expressions such as $w/(w+r)$ or $r/(w+r)$ when $w=r=0$. Another example would be $(r+w-3)!/(r-1)!$ for $r=0$ and $w=2$ and so forth. The particular ordering of the limits implies that for the base case $r=0$ and $w=0$ (the game starts without any red balls, i.e. the game ends before it begins), the probability $p_{III}^{(k)}(r=0,w=0)$ is, by our limit convention, equal to one when $k=0$ and zero otherwise. To rephrase, the ordering of the limits essentially defines equations (2b) and (2d) and hence, according to the discussion of the previous paragraph, such an ordering defines the boundary conditions of the problem, which are necessary for the uniqueness of the solution. We will see explicitly the necessity of the orderings of the limits and the applications of (2b)-(2e) in what follows.

For Problem 1, Rule IV, the roles of $r$ and $w$ are interchanged. In particular, we have

[TABLE]

Analogous comments apply here as the comments below equation (2).

2 A toy model

The following problem is given.

Problem 2

Inside a box containing $r$ red and $w$ white balls random removals occur. The balls are removed randomly one by one. The game ends once all the red balls are removed. The question posed is the determination of the probability that $k$ white balls remain where $k=0,\,1,\,2,\,...,w$.

2.1 Solution through $z$-transformations

Using the law of total probability by conditioning on the first removal, the following recursive relation is obtained

[TABLE]

We now make the educated ansatz

[TABLE]

obtaining the much simpler equation for $f^{(k)}$

[TABLE]

This ansatz is motivated because we want to factor out $\frac{r!w!}{(r+w)!}$, which provides the probability to remain with $k=w$ white balls; that is if we keep removing only red balls (equivalently only white balls).

In the next step, we multiply both sides of (6) by $z^{-r}$ and $u^{-w}$ and sum over $r$ and $w$ from one to infinity. Shifting some indices around and adding and subtracting suitable terms, we eventually obtain

[TABLE]

where we simplified the notation on $f^{(k)}$ setting $f^{(k)}(r,w)=f^{(k)}_{r,w}$.

Defining the generating functional $Y^{(k)}$ by

[TABLE]

equation (7) yields

[TABLE]

Case $k>0$

In this case, most of the terms in the right hand side of (9) are zero. In particular, $f^{(k>0)}_{r,0}=0\,\forall\,r\geq 0$ because the probability to remain with a positive number of white balls $k$ starting without any white balls at all must be zero. Another way to see it is through equation (2d). Also in the view of (2b) all $f^{(k)}_{0,w}$’s are zero except from the term $f^{(k)}_{0,k}$, which is equal to one. Hence,

[TABLE]

The $f^{(k)}_{r,w}$’s are then obtained through a Laurent expansion of (10) in inverse powers of $z$ and $u$. In fact, the expansion can be derived from the expansion of $z(u-1)/(zu-z-u)$ followed by a shifting of the $w$ index of $u$ by $k$ due to the overall $u^{-k}$ factor in (10). Hence,

[TABLE]

where in the third line we have extended the summation from $w=0$ because the denominator of the second term in the bracket vanishes assuming the factorial is analytically continued to a Gamma function such that $n!=\Gamma(n+1)$. Also the term $w=r=0$ in the sum of the last equality should be understood, in the view of (2b) and (2c), as taking the limit $w\rightarrow 0$ first followed by $r\rightarrow 0$ yielding a unit coefficient. Generally, the previous arguments employ the ordering of the limits defined in (2a) (see also the discussion below (2)).

Combining now (11) with (10) we obtain

[TABLE]

where the in the second equality we have again extended the sum from $w=0$ in the view of the fact that $\Gamma(-n)=\infty,$ when $n\,=\,0,\,1,\,2,\,...\,.$ From (12) and (8) we identify the coefficients $f^{(k)}_{r,w}$. Using then (5) we obtain the desired probability formula

[TABLE]

Case $k=0$

In this case, the non-zero terms in the right hand side of (9) are the $f^{(k=0)}_{r,0}\,\forall\,r\geq 0$ including the $f^{(k=0)}_{0,0}$ term. In fact, all of these terms are equal to one in the view of (2d); they represent the probability to remain with $k=0$ white balls when we start without any white balls at all and this should happen with probability equal to one. The fact that $f^{(k=0)}_{r,0}=\delta_{0}^{0}=1$ has already taken into account the rescaling pre-factor of (5), which also reduces to one when $w=0$.

Resumming then the right hand side of (9) and simplifying yields

[TABLE]

Comparing the right hand side of (14) with that of (10) we realize that the probabilities for $k=0$ can be obtained from (13) by taking the limit $k\to 0$.

To conclude, the general solution is given by

[TABLE]

and it satisfies the normalization condition

[TABLE]

One may check that formula (15) satisfies (4) and has the right behavior in several limiting cases. Indeed, $p^{(k)}(r,w<k)=0$ where we use the fact that the factorial of negative integer has a first order pole (tends to infinity). Also, in the view of (2a), the probability formula yields $p^{(k)}(0,w)=\delta_{w}^{k}$.

2.2 A combinatorial derivation

In section 2.1 we followed a rather complicated route in order to illustrate part of the formalism of Problem 1. However, Problem 2 may also be solved using simple combinatorics, which we present below for comparison.

We begin by defining the events $A$ and $B$ where

[TABLE]

The following two key facts are emphasized: (i) The complement of the event $A$, $A^{\intercal}$, and the event to remain with $k$ white balls starting with $r$ reds and $w$ whites, are mutually exclusive. (ii) The event $B$ does not impose any constrains on the possible orderings where the $r-1$ red and the $w-k$ white balls could be chosen.

With these two key facts at hand we apply the law of total probability by conditioning on the events $A$ and $A^{\intercal}$. Defining $P(A)$ and $P(B)$ as the probabilities for the occurrence of the events $A$ and $B$ respectively we obtain

[TABLE]

which is identical to (15) and this completes the derivation.

It is interesting to note that the probability of the event $B$, according to the second equality of the second line of (18), essentially follows a type ${}_{2}F_{1}$ hypergeometric distribution222We remind the reader that “In probability theory and statistics, the hypergeometric distribution is a discrete probability distribution that describes the probability of $k$ successes (random draws for which the object drawn has a specified feature) in $n$ draws, without replacement, from a finite population of size $N$ that contains exactly $K$ objects with that feature, wherein each draw is either a success or a failure. If a random variable $X$ follows the hypergeometric distribution then its probability mass function is given by $p_{X}(k)=\frac{\binom{K}{k}\binom{N-K}{n-k}}{\binom{N}{n}}$” (Wikipedia).. In problem $1$, as we will see further below in sections 5, 6 and 7, the hypergeometric function is a key object.

For later comparison with the results of Problem 1, we provide the probability generating functional of Problem 2 for fixed $r$ and $w$, which is given by

[TABLE]

It is not a coincidence that equation (19) involves a type ${}_{2}F_{1}$ hypergeometric function.

3 Setting up the problem: From the law of total probability to the required PDE

We now return to Problem 1. Working as in previous section, we will derive a PDE rather than an algebraic equation for the analogous $Y^{{(k)}}$ we saw earlier.

Using the law of total probability by conditioning on the first removal, we obtain

[TABLE]

where the weighting probability of $p^{(k)}(r-1,w)$ is complement to the one of $p^{(k)}(r,w-1)$, and can be understood as $w/(w+r)+rw/(r+w)^{2}$: the probability to remove a white ball equals the probability to either chose a white ball directly or to choose a red followed by a white ball.

3.1 The PDE of the generating functional with rescaling

In this case too, we make an educated ansatz on (20) as follows

[TABLE]

obtaining the much simpler equation333Note added after the completion of this work: According to (F4), the ansatz $p^{(k)}(r,w)=\frac{(r!)^{2}k!w!}{((r+w)!)^{2}}f^{(k)}(r,w)$ yields an equation of the form (F8), which is an even simpler equation for $f^{(k)}(r,w)$ than equation (22). for $f^{(k)}$

[TABLE]

This ansatz, as in Problem 2, is motivated because the expression $(r!w!)^{2}/((r+w)!)^{2}$ provides the probability to remain with $k=w$ white balls; that is if we keep removing only red balls. Working then as in Problem 1 by multiplying both sides of (22) by $z^{-r}$ and $u^{-w}$ and summing over $r$ and $w$ from one to infinity, after some algebra, we eventually obtain

[TABLE]

Defining $Y^{(k)}(z,u)$ as in (8) and making the substitutions

[TABLE]

on the $r$ and the $w$ factors that multiply the quantity $f^{(k)}_{r,w}z^{-r}u^{-w}$, equation (23) yields

[TABLE]

In what follows, we proceed using standard techniques of PDEs. In particular, we proceed using the fact that a first order PDE in two dimensions is equivalent to a system of two first order ODEs. In order to see that we define and solve the first ODE

[TABLE]

where $x_{2}$ is the integration constant. The next step is to solve for $x_{2}=x_{2}(z,u)$ and make the change of variables

[TABLE]

where we note that as $z,\,u\to\infty$ we have that $x_{1}\to\infty$ and $x_{2}\to-\sqrt{z}\to-\infty$. In the final step, we apply the transformations (27) and the differential equation (25) becomes

[TABLE]

where, as expected, the PDE reduces to a first order ODE. This is the second equation in the equivalent set of the two first order ODEs mentioned earlier.

In the following sections, we will study (28) in three different versions depending on the right hand side. The first version, section 4, corresponds to known results in the literature (Problem 3), the second version, sections 5 and 6, corresponds to the original problem, Problem 1, Rule III, while the third version corresponds to section 7, Problem 1, Rule IV.

3.2 The PDE of the generating functional without rescaling

For completeness, we present the PDE that corresponds to the initial recursion (20) rather to (22), and, which provides the actual probability generating functional. Working as before and defining

[TABLE]

the PDE satisfied by $\tilde{Y}^{(k)}(z,u)$ is given by

[TABLE]

We now have two cases to deal with. We start with Rule III for which conditions (2b) and (2d) yield $p_{III}^{(k)}(0,w)=\delta_{w}^{k}$ and $p_{III}^{(k)}(r,0)=\delta_{0}^{k}$ respectively. Given this and using the fact that

[TABLE]

we conclude that for the Rule III case, equation (30) reduces to

[TABLE]

The second case is about Rule IV. Using an analogous approach as in the case of Rule III, it is eventually found that in the case of Rule IV, equation (30) reduces to

[TABLE]

where $D_{z,u}$ is the (same) differential operator defined in (32b).

We observe that the resulting PDEs are $2^{nd}$ order linear PDEs with mixed derivatives, which are generally harder equations to solve.

4 Reproducing known results from the literature

The following known problem [1] is given.

Problem 3

Inside a box containing $r$ red and $w$ white balls random removals occur. The balls are removed one by one according to the three following rules. Rule I: If a white ball is chosen it is immediately discarded from the box. If a red ball is chosen, it is placed back into the box and a new ball is randomly chosen. The second ball is then removed irrespective of the color. Rule II: Once a ball is removed, the game continues from Rule I. Rule III: The game ends once all balls, except from one ball that remains, are removed. The question posed is the determination of the probability that the remaining ball is white.

Solution

All the steps that lead to (28) are identical for this problem as well and hence, (28) is the starting point. Given Rule III, which in particular, it implies that even at the event where all reds are removed and where $k>1$ white balls remain, we are still allowed to keep removing white balls (with probability one) until we end up with a single white ball. Thus, we now drop the superscript $k=1$ from $Y^{(k)}$ and $f^{(k)}$ as, according to Rule III, we will keep removing balls until one remains. Therefore, $f_{0,0}=0$ and $f_{0,w}=1,\,\,\forall\,w>0$ 444In particular, equation (2b) ceases to apply. and hence, the right hand side of (28) becomes $-1/u=-(1+x_{1}+x_{2}\sqrt{x_{1}})/x_{1}$ (see (26) and (27)). Hence, (28) simplifies to

[TABLE]

where the boundary condition implies that when no white balls exist ($u\to\infty$), which is equivalent to $x_{2}\to-\frac{1}{\sqrt{x_{1}}}-\sqrt{x_{1}}$ (see (27b)), the probabilities and hence, the generating functional to remain with one white ball at the end, should vanish.

The general solution to (34) is given by

[TABLE]

where $C(x_{2})$ is an arbitrary integration constant to be specified by the boundary condition, which in particular, it implies

[TABLE]

where the simplification from the second to the third line uses the following

[TABLE]

The last simplification is very convenient when expanding in Laurent series.

Combining now (35) with (36) we arrive at the desired generating functional

[TABLE]

where $x_{1}$ and $x_{2}$, according to (27), are functions of $z$ and $u$. It is straightforward to check that (38) satisfies the set of equations in (34).

The final step is to use (38) and (27) and Laurent expand in inverse powers of $z$ and of $u$ noting from (27) that as $z,\,u\to\infty$ we have $x_{1}\to\infty$ and $\,x_{2}\to-\infty$. In order to expand we make use of the following identities

[TABLE]

where we use Gamma functions instead of factorials for reasons that will become evident in what follows. Using the expansions (39), the generating functional (38) expands as

[TABLE]

Changing the dummy indices $i$ and $n$ by setting $i=r-n$ with $r=0,\,1,\,2,\,...$ and (a new) $n=0,\,1,\,2,\,...,\,r$ equation (40) becomes

[TABLE]

Our approach is to use (41) and compute the summation term over $n$, which, according to (8), is by definition equal to $f_{r,w}$ (with the $k$ index dropped). With $f_{r,w}$ at hand, using (21) we will eventually obtain the required probability formula.

A comment is in order: both, the numerator and the denominator of the last fraction of the Gamma functions in (41) diverge because $n$, $w$ and $r$ are positive integers with $r\geq n$. This is why we have replaced the factorials with Gamma functions as an analytic continuation of the factorials. These infinities are fictitious and are a consequence of interchanging the summation orderings. Similar infinities will appear in sections 6 and 7. Hence, in order to proceed we need the behavior of the Gamma functions near their poles; the corresponding expansions are documented in appendix A. As we will show, the end result, when taking the limits carefully, will be finite and will reproduce exactly the correct coefficients such that the recursion (20) and the boundary condition $p(r,0)=0$ are both fulfilled, which is what matters at the end.

Using (A2a) in order to cancel the poles of the Gamma functions555We only consider the poles resulting from $r\to\mathbb{Z}_{\geq 0}$ because we ignore poles resulting from the dummy summation index $n\to\mathbb{Z}_{\geq 0}$. Another way to see it is to only consider the pole of the variable $-(r-n)\leq 0$. Furthermore, we ignore the pole in $\Gamma(-2(r-n)-n-w)$ coming from $w\to\mathbb{Z}_{\geq 0}$ (and from $n\to\mathbb{Z}_{\geq 0}$) because it is not compensated by a similar $w$ (or $n$) pole in the numerator and hence, it does not contribute a finite part. the last fraction of the Gamma’s in (41) is exchanged according to

[TABLE]

Simplifying then the $\Gamma(1+2r-2n)$ that appears in the numerator and in the denominator of the resulting step, keeping track of the correct overall sign as $(-1)^{n}$, and changing the dummy index $n$ setting $n\to r-n$, transforms the summand of (41) into $(-1)^{r+n}\Gamma(1+r+w+n)/\left((1+n)\Gamma(1+w)\Gamma(1+r-n)\Gamma^{2}(1+n)\right)$ where as before the (new) $n=0,\,1,\,...,\,r$. Performing then the summation of the resulting expression over $n$ using the summation identity

[TABLE]

proved in Appendix D, equation (41) eventually becomes666In Appendix C, we will briefly provide a second derivation of (44).

[TABLE]

The last step is to use (8) in order to identify the $f_{r,w}$ coefficient from (44) and multiply it by the overall ansatz factor from (21) to obtain

[TABLE]

Equation (45) is the main result of the section and reproduces the result of [1] obtained through a completely different (inductive) approach.

5 Deriving the generic generating functional for any number of remaining white balls $k$, Rule III

Having solved Problem 3, which uses all the ingredients of our proposed approach, we gained confidence in order to attack Problem 1. In the following two sections, we simplify the notation by dropping the index III from $Y_{III}^{(k)}(x_{1},x_{2})$, $f^{(k)}_{III;r,w}$ and $p_{III}^{(k)}(r,w)$; we will restore it at the very end on the final formula for $p_{III}^{(k)}(r,w)$.

Starting from (28), we write the differential equation and the boundary conditions for Problem 1, Rule III. Noting that $f^{(k)}_{0,w}=\delta^{k}_{w}$, because $f^{(k)}_{0,w}$ denotes the probability to remain with $k$ white balls if we start with no red balls and with $w$ white balls (see (2b)), we obtain

[TABLE]

where we have indicated explicitly that $z$ and $u$ are functions of $x_{1}$ and of $x_{2}$ (see (27)).

The boundary condition $z\to\infty$ (game starts without any red balls) is a boundary case where the game ends before it (even) begins. Then, the probability to remain with $k$ white balls should be equal to one if we start with $k$ white balls and zero otherwise 777The argument takes into account the ansatz factor in (21), which also reduces to one, when $r=0$.. The boundary condition in the generating functional representation, equation (46b), is a manifestation of equation (2b), which is the boundary condition in the probability representation.

The solution of (46) that satisfies the boundary condition turns out to be

[TABLE]

where in (47a) we simplified using (37), and where ${}_{2}F_{1}$ in (47b) is a hypergeometric function.

It is a straightforward but also a tedious process to show that (47) satisfies the differential equation (46) and a computer program such as “Mathematica” is recommended. The fact that $Y^{(0)}$ fulfills the boundary condition can be verified by taking first the limit $\lim{x_{2}\to-\sqrt{z}}(1-\frac{1}{u})$ followed by the limit $\lim{z\to\infty}$. These limits are a consequence of (27) and of the boundary condition (46b).

In order to see that the boundary condition for $Y^{(k>0)}$ is also satisfied we take the same limits as the limits we took for $Y^{(0)}$. Then, the overall multiplicative factor of (47b) decays as $-(u-1)^{-3-2k}u^{2}/2\times 1/z^{1+k}$. Expanding the first term of the curly bracket (not the one involving the ${}_{2}F_{1}$) to leading order in $z$ shows that this quantity grows as $-2(u-1)^{2k+2}/u^{k+1}\times z^{k+1}$. Next, we move to the second term in the curly bracket, that is the hypergeometric function whose argument grows as $-(1-u)^{2}/u\times z$. Taking into account that if $k$ is a positive integer, we find that the ${}_{2}F_{1}\left(-1-2k,-k;-2k;-x\right)$ becomes a terminating polynomial of degree $k$, which grows as $(2k+1)/(k+1)x^{k}$ as $x\to\infty$. Thus, the whole term involving the ${}_{2}F_{1}$ grows as $2(1-u)^{2k+2}/u^{k+2}\times z^{k+1}$. Therefore, the curly bracket grows as $-2(u-1)^{2k+3}/u^{k+2}\times z^{k+1}$ and given the overall pre-factor $-(u-1)^{-3-2k}u^{2}/2\times 1/z^{1+k}$ it is concluded that (46b) converges to $1/u^{k}$ as should. We have thus just shown that the boundary condition is fulfilled $\forall\,k>0$.

In the following section we will finalize the solution of Problem 1, Rule III. Using (47), we will compute explicitly the probabilities $p^{(k)}(r,w)$ for the first few values of $k$ and from there we will derive the formula for any $k$.

It is also interesting to note that the alternative boundary conditions

[TABLE]

yield the same exact solutions provided by equations (46). In particular, the boundary conditions defined by (48) are the continuum version of the discrete boundary conditions provided by equation (2d).

These boundary conditions imply that as $u\to\infty$ (game starts without any white balls), which due to (27) is equivalent to $x_{2}\to-\sqrt{x_{1}}-\frac{1}{\sqrt{x_{1}}}$, the probability to remain with $k>0$ white balls is zero. If, on the other hand, $k=0$ then the probability to remain with no white balls starting without any white balls, must be equal to one independently on the initial number of red balls 888The argument takes into account the ansatz factor in (21), which also reduces to one, when $w=0$.. Indeed, expanding the expression $x_{1}/(x_{1}-1)$ in inverse powers of $x_{1}$ yields unit coefficients 999We remind the reader that $x_{1}=z$ and that the variable $z$ corresponds to the red balls.. We reserve the details, which show that the boundary conditions (48) are satisfied by (46), for appendix B.

In case the reader wonders about the uniqueness of the problem, the answer is that the two boundary conditions are equivalent because they define the same exact problem. In particular, one may choose to evolve the initial data starting from the surface $z\to\infty$ or to evolve the initial data starting from the surface $u\to\infty$, as long as the chosen boundary surface is known. In this problem, it happens that we know both boundaries. In what follows, things will become clearer and the simultaneous fulfillment of both boundary conditions will be checked once explicit formulas are obtained.

6 The general probability formula for Rule III

Having developed all the necessary machinery we are now in the position to answer the question of Problem 1, Rule III. The strategy we will apply is as follows. We first consider (47) and take the limits $k\to 0,\,1\,,\,2$. Then, noting from (27) that as $z,\,u\to\infty$ we have $x_{1}\to\infty$ and $x_{2}\to-\infty$, we Laurent expand the solutions at infinity, that is in inverse powers of $x_{1}=z$ and $x_{2}$. Then, we re-expand $x_{2}$ in inverse powers of $z$ and $u$ in an analogous way as we worked in Problem 3, section 4. The final step is to identify the coefficients of $1/z^{r}$ and $1/u^{w}$ as the $f_{r,w}^{(k)}$’s of (8) and from there, using (21), we will obtain the $p^{(k)}(r,w)$. Using then the explicit expressions for $p^{(0)}(r,w)$, $p^{(1)}(r,w)$ and $p^{(2)}(r,w)$ and by studying the pattern, we will make an ansatz for the generic solution $p^{(k)}(r,w)$ for any $k$. The last step is to check that the ansatz fulfills the initial recursion we begun with, namely equation (20), and the condition $p^{(k)}(0,w)=\delta_{w}^{k}$. It is also noted that in principle we could had reached the same answer without making any anstaz at all. That is through an even more systematic approach; by Laurent expanding the general PDE solutions (47) for arbitrary $k$.

6.1 The case $k=0$

Despite this case can be obtained using the strategy outlined in the introduction of the current section, there is a way around it using the results of section 4. In particular, the probability of Problem 1, Rule III for $k=0$ is equal to the complement of the probability of Problem 3, which is the probability that the last ball is red; that it no white balls remain at the end of the game. Thus, using (45), we find

[TABLE]

In the view of (2a) and (2c), equation (49) can be written compactly as

[TABLE]

As a cross check of (50), we substitute (27) into (47a) and expand the result at $z=\infty$ and at $u=\infty$ up to $4^{th}$ order to obtain

[TABLE]

We then compare previous expansion with the quantity

[TABLE]

observing an exact matching order to order in $z$ and in $u$. It is also notable that both boundary conditions (46b) and (48) are satisfied. Indeed, checking the condition (46b) is obvious. On the other hand, as $u\to\infty$, the expansion reduces to $1+1/z+1/z^{2}+1/z^{3}+....=z/(z-1)=x_{1}/(x_{1}-1)$ in accordance with the (alternative) boundary condition in (48). Thus, as discussed below (48), we can see explicitly the equivalence of the two boundary conditions.

Equation (50) will also be obtained as the $k\to 0$ limit of the general solution $p^{(k)}(r,w)$.

6.2 The case $k=1$

Taking the limit in (47b) as $k\to 1$ yields

[TABLE]

Using (27), the first line of (53) expands as

[TABLE]

while the second line, after some algebra, expressed as an expansion in $x_{1}$ and $x_{2}$ yields

[TABLE]

From here and on we work analogously to equations (39)-(41). In particular101010The explanations that follow assume an identical indexing in $i$, $n$ $r$ and $w$ as the indexing of section 4., we expand $\frac{1}{x_{2}^{2i-1}}$ using (39b), exchange the index $i$ with $i=r-n+1$ where $r=1,\,2,\,...$ and $n=0,\,1,\,2,\,...\,r-1$ as before to finally obtain

[TABLE]

We note that for $r\geq 1$, $n<r$ and $w\geq 0$ the numerator and the denominator in the last fraction of the Gammas in (56) diverges. Using then (A2b), which allows us to replace the last fraction of the Gamma functions with $(-1)^{n+w}\Gamma(1+2r-n+w)/\Gamma(1+2(r-n))$ 111111We only consider the poles resulting from $r\to\mathbb{Z}_{\geq 0}$ because we ignore poles resulting from the dummy summation index $n\to\mathbb{Z}_{\geq 0}$. Furthermore, we ignore the pole in $\Gamma(-2r+n-w)$ coming from $w\to\mathbb{Z}_{\geq 0}$ because it is not compensated by a similar $w$ pole in the numerator and hence, it does not contribute a finite part., and equation (8), which defines the $f^{(k)}_{r,w}$ coefficients, we obtain

[TABLE]

where in the second equality we changed the dummy index according to $n\to r-n$ and where $|_{r,w}$ denotes the coefficient of $\frac{1}{z^{r}}\frac{1}{u^{w}}$ in the preceding bracket. Equation (57) can be partitioned in four cases according to

[TABLE]

where (58a) is obtained from (57) for $(r,w)=(1,0)$, and it has exactly the form we need in order to cancel the undesired $1/(3z)$ term coming from (56). Equation (58c) is obtained by performing the summation using the identity

[TABLE]

which can be proved working along the lines of equations (D1)-(D5) with minor modifications. We note that in the limit $(r,w)\to(0,1)$, the right hand side of equation (58c) tends to one121212The limit $r=0,\,w=1$ in (58c) must be taken according to equation (2) (see also comments that follow (2)). In particular, when taking the first limit $w\to 1$, both of the $(r-1)!$ terms that appear in the numerator and the denominator of (58c) cancel out yielding $f^{(1)}_{0,1}=1$., and hence, (in theory) it provides the $1/u$ term (see boundary conditions, (46b)); in reality, the $1/u$ term is actually coming from \left(\frac{1}{u}-\frac{1}{3z}\right)\big{|}_{r=0,w=1}. Thus, equation (58b) can be absorbed in (58c) in the view of the limits ordering of (2a), and by extending the range of validity suitably.

Collecting all the terms of (53) using (58) and the observations we just made above, we finally obtain

[TABLE]

where we extend the summations on $r$ and on $w$ according to $r,\,w\geq 0$ in the view of the fact that $1/(r-1)!$ and $1/(w-1)!$ tend separately to zero as $r\to 0$ and $w\to 0$ respectively. As noted earlier, when $(r,w)=(0,1)$ we get $f^{(k)}_{0,1}=\delta_{1}^{k}$ (see (2) and the discussion below (2) in order to see how these limits should be taken). The reason we expect $f^{(1)}_{0,1}=1$ is because this coefficient encodes the $1/u$ term (see (46b) and (54)) and hence, equation (60) provides the right coefficient with value equal to one, which is precisely the probability to remain with one white ball if we start with one white ball and no red balls. To rephrase, if we start with no red balls $(r=0)$ and $w\neq 1$ white balls the probability to remain with $k=1$ white balls is zero, otherwise if $w=1$ the probability is equal to one, and this is precisely the meaning of $f^{(k)}_{0,1}=\delta_{1}^{k}$.

The last step is to use (8) in order to identify $f^{(k)}_{r,w}$ from (60), multiply by the coefficient of (21) and do the necessary simplifications to obtain

[TABLE]

which also includes the cases $w=0,\,\,r\neq 1$ and $(r,w)=(0,1)$ in the view of (2d) and (2b) respectively. One may verify that equation (61) satisfies the initial recursion we begun with, namely equation (20) and the boundary condition (2b), and in fact, all the equations (2b)-(2e).

Comparing (50) with (61) we already start to see a pattern forming. In the next section, once the $k=2$ case is computed, the pattern will become obvious.

6.3 The case $k=2$

Taking the limit in (47b) as $k\to 2$ yields

[TABLE]

Before computing the general $f^{(k)}_{r,w}$ term, we expand (62) up to $4^{th}$ order in $1/z$ and $1/u$ as a way to cross-check our calculations along the way. The result is

[TABLE]

As expected, according to the boundary condition (46b), the lowest order in $1/u$ is $1/u^{2}$, which basically says that unless we start with at least two white balls, the probability to end up with two white balls is zero. More specifically, we know that if $(r,w)=(0,2)$ we should get a probability equal to one, which is precisely equal to the coefficient of the $1/z^{0}1/u^{2}$ term in the expansion. We also note that as $u\to\infty$, $Y^{(2)}(z,u)\to 0$ in agreement with the alternative boundary condition (48).

Returning now to equation (62) and expanding it yields

[TABLE]

where the sum corresponds to the last two lines of (62). We note that we already have the desired $1/u^{2}$ term. Working then analogously to sections 4 and 6.2, and using the Gamma function identities

[TABLE]

and expanding $x_{2}$ as in (39b) we find that (64) yields

[TABLE]

We note that for any $w$, when $r=0$, and hence $n=0$, the summation term over $n$ becomes trivial and also it yields a zero result. The reason is because the numerator has only a single pole in $r$ coming from $\csc(2\pi(r-n))$ while the denominator has a second order pole in $r$, one coming from $\Gamma(r-n)$ and one coming from $\Gamma(-2r-w+n)$. Hence, the summation over $r$ should begin from $r=1$.

As it turns out the summation over $n$ is, up to overall Gamma function factors, a regularized hypergeometric function with integer coefficients that reduces to a hypergeometric function ${}_{2}F_{1}$. In particular, (66) yields

[TABLE]

where in the second equality we used equation (A2b) on $\Gamma(-2r-w)$ and the fact that as $r\to j\in\mathbb{Z}$ then $\pi\csc(2\pi r)\approx\frac{1}{2(r-j)}$. These two facts allow the cancellation of the poles in the ratio $\pi\csc(2\pi r)/\Gamma(-2r-w)$ and its replacement with $-\Gamma(2r+w+1)(-1)^{w}$ 131313We only consider the pole resulting from $r\to\mathbb{Z}_{\geq 0}$ and we ignore the pole in $\Gamma(-2r-w)$ coming from $w\to\mathbb{Z}_{\geq 0}$ because this pole is not compensated by a similar $w$ pole in the numerator and hence, it does not contribute a finite part.. Given also that ${}_{2}F_{1}\left(-4,0;-2;x\right)=1$, ${}_{2}F_{1}\left(-4,0;-3;x\right)=1$, ${}_{2}F_{1}\left(-5,-1;-4;x\right)=1-5/4x$ and ${}_{2}F_{1}\left(-5,-1;-5;x\right)=1-x$, the curly bracket in (67) yields the following terms

[TABLE]

Equation (68f) is not obvious but it will become soon due to equation (71b), which implies that for $w=\{0,1\}$ and $r>2$, the $\Gamma(w-1)$ term in the denominator diverges without any compensating factor coming from the numerator; thus the whole term is equal to zero. It is noted that the first three sub-equations of (68) have precisely the right form in order to cancel the unnecessary terms from (67) (see (63)), which thus, in the view of (8), yields

[TABLE]

Next, we use the Gauss theorem on the hypergeometric function ${}_{2}F_{1}$, which states that

[TABLE]

Strictly speaking in our case $c-a-b=2-w\leq 0$ for $w\geq 2$ and hence, the Gamma functions, of both, the numerator and the denominator, will diverge. Hence, we take the limits of the Gamma functions carefully by taking care of the poles that appear using (A2)141414While taking the required limits, we cross-check by matching the $f^{(2)}_{r,w}$’s with the coefficients of (63).. The end result is finite and is given by

[TABLE]

Equation (71) reproduces the $\frac{1}{z^{r}}\frac{1}{u^{w}}$ coefficients of (63), and this provides confidence that the equation is correct; it can be written compactly as

[TABLE]

where we have extended the range of $r$ and of $w$ by observing that: (i) If $r=0$ and $w=2$, taking the limits in (72) according to (2a), yields $f^{(2)}_{0,2}=1$ in agreement with (2b). (ii) Else if $r=0$ and $w>2$ we have one $r$ pole in the denominator coming from $\Gamma(r)$ that is not compensated by a similar pole in the numerator and hence, the result is zero. (iii) Else if $r=0$ and $w=0,\,1$, despite the ratio $\Gamma(r+w-2)/\Gamma(r)$ is finite, the $1/\Gamma(w-1)$ term tends to zero and hence, $f^{(2)}_{0,1}=f^{(2)}_{0,2}=0$. (iv) Else, if $w=0,\,1$ and $r\geq 1$ the result is also zero, including the cases $(r,w)=\{(1,0),\,(1,1)\}$, in the view of the orderings of the limits of (2a) and the fact that the $1/\Gamma(w-1)$ in (72) goes to zero as $w\to 0,\,1$. Given that we are dealing with $k=2$ and $w=0,\,1$, we see that case (iv) is consistent with equation (2e) as should.

The last step is to multiply (72) by the coefficient of (21) and do the necessary simplifications to obtain

[TABLE]

It can be checked that equation (73) satisfies the initial recursion (20) and the boundary condition (2b), and in fact, all the equations (2b)-(2e). This completes the derivation of $p_{III}^{(2)}(r,w)$. In the following section, we will provide the general solution $p_{III}^{(k)}(r,w)$ $\forall k$.

6.4 The general solution

Having computed $p^{(k)}(r,w)$ for $k=0,\,1$ and $2$, we are now in position to attempt for a general ansatz solution. Observing equations (50), (61) and (73) one may guess that the solution should have the form

[TABLE]

up to an overall constant that does not depend on $r$ and $w$. The constant is then specified by the requirement that $p^{(k)}(0,k)=1$ in the view of the boundary condition (2b) from where we infer that the solution must be

[TABLE]

The final step is to verify that (75) is the required probability solution. Indeed, it is a matter of straightforward algebra to show that (75) satisfies the initial recursion we begun with, namely equation (20), and all the equations (2b)-(2e).

We do the following two cross-checks: (i) We first note that the solution reproduces the cases $k=0,\,1,\,2$, equations (50), (61) and (73). (ii) Also, as it is shown in Appendix E, the probability formula is normalized and hence, it satisfies

[TABLE]

which basically says that when the game ends, we will surely end up with a number of white balls between zero and the initial number $w$. We note that the normalization equation (76) comes out automatically and this key fact serves as another cross check of the correctness of (75). At this stage, the solution is considered as complete.

Two important observations can be made: (a) It is notable that (75) is, up to the rescaling factor $\frac{k!r!(r+w+1)}{(r+k+1)!}$, the same as the solution of Problem 2, equation (15). This implies that this overall factor, in a sense, encodes the additional complication of adding red balls back into the box. (ii) One can also show that the (true) generating probability functional $\tilde{Y}^{(k)}(z,u)=\tilde{Y}^{(k)}_{III}(z,u)$ defined in (29) satisfies the PDE (32) with $p^{(k)}(r,w)=p_{III}^{(k)}(r,w)$ given by the probability formula (75).

We also compute the probability generating functional for fixed $r$ and $w$. It is is given by

[TABLE]

which involves a $\,{}_{3}F_{2}$ hypergeometric function rather than a $\,{}_{2}F_{1}$, which was the case for Problem 2 (see (19)).

6.4.1 Maximal probability

In this section we investigate the maxima of $p_{III}^{(k)}(r,w)$ as a function of $k$ for fixed values of $r$ and of $w$. We start by constraining $k$ solving the two inequalities

[TABLE]

which provide the set of $k$’s for which $p_{III}^{(k)}(r,w)$ has local maxima (if any). The inequalities (78a) and (78b) imply

[TABLE]

respectively. Evidently the system of inequalities has no solution. This hinds that, for $r>0$, the maximal probability is at the boundary cases $p_{III}^{(0)}(r>0,w)$ or $p_{III}^{(w)}(r>0,w)$. In fact, we guess that the $p_{III}^{(0)}(r>0,w)$ should be the required maximal probability. In order to show the claim, we take the difference $p_{III}^{(k)}(r,w)-p_{III}^{(1+k)}(r,w)$ and simplify to obtain

[TABLE]

where $|C|$ is a positive constant of ratios of factorials. The difference $p_{III}^{(k)}(r,w)-p_{III}^{(1+k)}(r,w)$ is minimum when $k=w-1$, and given that $r\geq 1$, it implies that the difference is always positive and hence, $p_{III}^{(k)}(r,w)>p_{III}^{(1+k)}(r,w)$. Therefore, when $r>0$ we find that $p_{III}^{(k)}(r,w)$ is a monotonically decreasing function of $k$. On the other hand, we know that for $r=0$, $p_{III}^{(k)}(0,w)=\delta_{w}^{k}$. Thus, we have just proved

[TABLE]

To conclude, the maximal probability occurs for: (i) $k=w$ if $r=0$ and is given by $p_{III}^{(k=w)}(0,w)=1$. (ii) Else, the maximal probability occurs for $k=0$ and the corresponding probability $p_{III}^{(k=0)}(r,w)$ is given by equation (50).

6.4.2 Limiting cases

It is also interesting to investigate the behavior of (75) in the limits $r\to\infty$ for $w$ fixed and any $k$, and for $w\to\infty$ for $r$ and $k$ fixed. For this purpose we use the more convenient equation (74) from where it is deduced that

[TABLE]

where $\prod_{i=0}^{k-1}$ for $k=0$ is defined to be zero.

The asymptotic expansion (82a) says that in the limit $r\to\infty$ with $w$ kept fixed, the probability $p_{III}^{(k>0)}(r,w)$ decays as $\frac{1}{r^{2k}}$. Moreover, if $k=0$, the probability $p_{III}^{(k=0)}(r,w)$ tends to one, which means that when the red balls are much more than the white balls, the game will (most likely) end without any white balls. Both of these two asymptotic results behave as expected. On the other hand, the asymptotic expansion (82b) in the limit $w\to\infty$ with $r$ and $k$ kept fixed yields less expected results. In particular, the probability $p_{III}^{(k)}(r,w)$ (to leading order in $w$) becomes independent on the initial number of the white balls $w$.

7 The general probability formula for Rule IV

In this section we solve Problem 1, Rule IV. The ideas are similar as those of Problem 1, Rule III and hence, the derivations are sketchy and much shorter compared to those of sections 5 and 6. In particular, the proofs between several steps are similar or even identical to the aforementioned sections and therefore, they are omitted.

7.1 The PDE, the boundary conditions and the solution, Rule IV

The starting point is the PDE (28) where the coefficients in the sums in the right hand side of the equation are all equal to $f_{0,w}^{(k)}=\delta_{0}^{k}$ (see (3d)). The (physical) reason is because $f_{0,w}^{(k)}$ is (proportional to) the probability to remain with $k$ red balls if we start without any and this should yield probability equal to one if $k=0$ and equal to zero otherwise. Summing then the two series over $w$ we find that the result is zero $\forall k$, including $k=0$. Hence, the PDE reduces to

[TABLE]

where the boundary condition (83b) is explained as follows. The condition $x_{2}\to-\sqrt{x_{1}}-\frac{1}{\sqrt{x_{1}}}$, according to (27), is equivalent to the condition that all the white balls are out of the box ($u\to\infty$). Then, in this case, the only solution for $k$ red balls to remain is if inside the box exist exactly $k$ red balls. In other words, the lowest order term in the Laurent expansion of $Y^{(k)}(x_{1},x_{2})$ should be $\sim 1/z^{k}$. The solution then to (83) is given by151515Analogously to the Rule III case, we could had imposed the alternative boundary conditions $\lim_{x_{1}\to\infty}Y^{(k)}(x_{1},x_{2})=\frac{u}{u-1}\delta_{0}^{k}$ and still obtain the same $Y^{(k)}(x_{1},x_{2})$ given by equation (84).

[TABLE]

7.2 The probabilities for the first few cases, Rule IV

In this section we Laurent expand (84) for $k=0,\,\,1,$ and $2$ and obtain explicit formulas for $p_{IV}^{(k)}(r,w)$. In particular, we first expand in inverse powers of $x_{2}$ and then re-expand in inverse powers of $z$ and of $u$ identifying the coefficients with the $f^{(k)}_{r,w}$’s (see (21)).

7.2.1 The case $k=0$

In fact, the solution to this case is already found in section (4) and is (almost) given by equation (45). The precise answer is

[TABLE]

where we extend the applicability of the formula in order to include $(r,w)=(0,0)$ in the view of (3c) and the discussion below equation (3). One can check that equation (85) satisfies the recursion (20), and also fulfills the boundary condition (3b).

7.2.2 The case $k=1$

It can be shown that (84) for $k=1$ expands as

[TABLE]

The next step is to follow exactly the same steps that led from equation (55) to equation (60) to eventually obtain

[TABLE]

One can check that equation (87) satisfies the recursion (20), and also fulfills the boundary condition (3b).

7.2.3 The case $k=2$

It can be shown that (84) for $k=2$ expands as

[TABLE]

The next step is to work as in the $k=1$ case by following exactly the same steps that led from equation (55) to equation (60) to eventually obtain

[TABLE]

One can check that equation (89) satisfies the recursion (20), and also fulfills the boundary condition (3b).

7.3 The general solution

Equations (85), (87) and (89) motivate the following ansatz for the general solution for any k

[TABLE]

which can be re-written as

[TABLE]

One may check that (91) satisfies both, the recursive equation (20) and the boundary condition (3b), and in fact, all the consistency-check equations (3b)-(3e). As a cross check, we verified that equation (91) reproduces (85), (87) and (89) and most importantly, it satisfies the normalization condition

[TABLE]

We note that the normalization equation (92) comes out automatically. This key fact serves as another cross check of (91). At this stage, the solution is considered as complete.

Two important observations can be made: (a) It is notable that (91) is, up to the rescaling factor $\frac{(2k+1)r!(r+w+k)!}{(r+k+1)!(w+r)!}$, the same as the solution of Problem 2, equation (15) with the (expected) reflection $r\leftrightarrow w$. This implies that this overall factor, in a sense, encodes the additional complication of adding red balls back into the box. (ii) One can also show that the (true) generating probability functional $\tilde{Y}^{(k)}(z,u)=\tilde{Y}^{(k)}_{IV}(z,u)$ defined in (29) satisfies the PDE (33) with $p^{(k)}(r,w)=p_{IV}^{(k)}(r,w)$ given by the probability formula (91).

We conclude the section by providing the probability generating functional

[TABLE]

for fixed $r$ and $w$. The generating functional involves a linear combination of $\,{}_{3}F_{2}$ hypergeo- metric functions rather than a single $\,{}_{3}F_{2}$, which was the case for Problem 1, Rule III (see (77)).

7.3.1 Maximal probability

This section investigates the set of $k$’s for which $p_{IV}^{(k)}(r,w)$ is maximized. For the boundary case $w=0$ we have $p_{IV}^{(k)}(r,0)=\delta_{r}^{k}$ (see (3b)) and hence, the probability is maximized for $k=r$ with maximal value $p_{IV}^{(r)}(r,0)=1$. The other boundary case is when $r=0$ in which case $p_{IV}^{(k)}(0,w)=\delta_{0}^{k}$ (see (3d)) and hence, the probability is maximized for $k=0$ with maximal value $p_{IV}^{(0)}(0,0)=1$. Otherwise, if $w>0$ and $r>0$, we work analogously to section 6.4.1 and we find that the set of $k$’s that maximize (locally) the $p_{IV}^{(k)}(r,w)$ is specified by the two inequalities

[TABLE]

We observe that the inequalities are linear in $k$, which means that (modulus degeneracies) there exists only one unique maximum. Before investigating the maxima, we define $k_{0}$ by

[TABLE]

Given that the investigation here is for $r\geq 1$ and $w\geq 1$, it is deduced that $k_{0}>1$. We now want to constrain $k_{0}$ from above. As it can be shown, there are two cases: Case 1. If $w=1$, then $k_{0}=r+1$ and hence, the maximum occurs for $k=r$. Case 2. Else if $w>1$, then $k_{0}<r+1$. Case 2 is partitioned into two sub-cases: Case 2a. If $k_{0}\in\mathbb{Z}_{+}$ then there are two maxima that correspond to $k=k_{0}-1$ and to $k=k_{0}$, which occur with equal probabilities. Case 2b. Otherwise if $k_{0}\notin\mathbb{Z}_{+}$ then there exists a unique maximum that corresponds to the integer between $k_{0}-1$ and $k_{0}$; that is for $k=\left\lfloor{k_{0}}\right\rfloor$ (the integer part of $k_{0}$). This is always the case for $r=1$ and $w>1$; the maximum occurs for $k=1$ because $k_{0}\notin\mathbb{Z}_{+}$ with $k_{0}\in(1,2)$.

Collecting all cases together, we have just proved

[TABLE]

This ends our investigation of the maxima of $p_{IV}^{(k)}(r,w)$.

7.3.2 Limiting cases

In this section we investigate the behavior of (91) in the limits $w\to\infty$ for $r$ fixed and any $k$, and for $r\to\infty$ for $w$ and $k$ fixed. For this purpose we use the more convenient equation (90) from where it is deduced that

[TABLE]

The asymptotic expansion (97a) says that in the limit $r\to\infty$ with $w$ and $k$ kept fixed, the probability $p_{IV}^{(k)}(r,w)$ decays as $\frac{1}{r^{2}}$ independently on $k$ 161616Modulus the $2k+1$ overall factor. More precisely, in this limit, the probabilities grow linearly in $k$ according to $(2k+1)p_{IV}^{(0)}(r,w)$; they grow as integer multiples of a fundamental quantity, the probability $p_{IV}^{(0)}(r,w)$. contrary to the case of Rule III, equation (82a). Moreover, the asymptotic expansion (97b) in the limit $w\to\infty$ with $r$ kept fixed is also interesting. In particular, the probability $p_{IV}^{(k)}(r,w)$ (to leading order in $w$) becomes independent on the initial number of the white balls $w$ analogously to the behavior of the Rule III case, equation (82b). In particular, if $r=1$, then $p_{IV}^{(0)}(1,w)\approx p_{IV}^{(1)}(1,w)\approx\frac{1}{2}$ for $w\gg 1$.

8 Discussion

Wow!! That has been a long but joyful journey for such a rather simple problem. One could think that our approach amounts in shooting a bug with a bazooka! While such a view might be right, in this work we have seen a concrete example of how a recursion equation arising from a probability problem, without a systematic method to approach, was transformed into a differential equation that was solved using standard PDE methods. Hence, we have seen a concrete example of how difference equations, such as those arising from discrete random processes, and for which the tools in the literature are less developed, could be transformed into differential equations, where the literature is rich and well studied. In particular, the approach we followed has been algorithmic and it may be summarized by the following recipe:

Step 1. Usage of the total probability law in order to derive a difference equation (rec-

ursion) in one or more variables. Impose suitable discrete boundary conditions.

Step 2. Apply z-transformations in order to transform the recursion into a differential

equation for the probability generating functional.

Step 3. Translate the discrete level boundary conditions of Step 1 into continuum

boundary conditions at the (P)DE level.

Step 4. Solve the (P)DE enforcing the boundary conditions of step 3.

Step 5. Laurent expand the solutions taking care of any fictitious infinities that may

appear by taking suitable limits. Hence, obtain the required probabilities.

This recipe is adaptable to large classes of probability problems, and generally for problems involving difference equations. For instance, adapting the ideas of Problem 1, Rule III, to those of Rule IV, whose results have been obtained fast and straightforwardly, has been effortless.

Retracing the steps (for Problem 1, Rule III), we started from a probability question and derived a recursion, equation (20), with two variables using the law of total probability. Then, through a 2D $z-$transformation (8), we eventually derived a first order 2D PDE (25) with suitable source terms and boundary conditions (46) whose solution, up to a factor (see (21)), provides the generating probability functional of the problem at hand for any number of remaining white balls $k$. Given that a first order 2D PDE is equivalent to a set of two first order ODEs, the PDE, through suitable coordinate transformations (see (27)), was reduced into two decoupled first order ODEs, equations (26) and (25). This step, in a sense, decouples the 2D problem into two 1D problems of suitable scaling variables, which are functions of the initial ones. While the two scaling variables are decoupled at the DE level, they are yet coupled in a non-trivial way through the boundary conditions (i.e. see (47)). Solving then the differential equation, and hence, obtaining the general generating functional (47), we performed a closed form Laurent expansion. Using the formula from the expansion, which provides the coefficients $f^{(k)}(r,w)$ (see (8) and (21)), we were eventually able to find the general probability formula, equation (75). Finally, we checked that the probability formula satisfies the initial recursion equation we begun with, and its boundary conditions (2b), thus completing the solution of the problem. An analogous approach has been followed in Problem 1, Rule IV.

It is notable how a relatively simple probability problem to state and to understand, has involved such a heavy computational machinery. The most challenging step has been to Laurent expand the generating functional and to deal with the fictitious infinities that appeared on the way due to the interchanging of the summation orderings. These infinities showed up as poles of Gamma functions; this part has been considerably more involved than solving the PDE itself.

We have also seen that the probability generating functionals involve hypergeometric functions (see (77)) or linear combinations of hypergeometric functions (see (93)) of type ${}_{3}F_{2}$. These results generalize the simpler version of the problem, which does not require placing balls back into the box, and, which involves a ${}_{2}F_{1}$ distribution instead (see (19)). Hence, in a sense, the additional complication of adding balls back into the box is captured by extending the ${}_{2}F_{1}$ type of the probability generating functional into a ${}_{3}F_{2}$ type.

As a bonus, we found the power series solution of the rather complicated family of PDEs given by equation (32), which satisfy the same boundary conditions $\tilde{Y}_{III}^{(k)}(z\to\infty,u)=1/u^{k}$ as the boundary conditions satisfied by $Y_{III}^{(k)}(z,u)$ (see (46b)). The solution is given by $\tilde{Y}_{III}^{(k)}(z,u)=\sum_{r,w=0}^{\infty}p_{III}^{(k)}(r,w)\frac{1}{z^{r}}\frac{1}{u^{w}}$. Likewise, for Rule IV, $\tilde{Y}_{IV}^{(k)}(z,u)=\sum_{r,w=0}^{\infty}p_{IV}^{(k)}(r,w)\frac{1}{z^{r}}\frac{1}{u^{w}}$ is the solution of the PDE (33), which satisfy the same boundary conditions $\tilde{Y}_{IV}^{(k)}(z,u\to\infty,u)=1/z^{k}$ as the boundary conditions satisfied by $Y_{IV}^{(k)}(z,u)$ (see (83b) and (27)).

It would be interesting to find simple probability problems, such as the ones presented in this work, whose differential equation representation is the same as that of known problems from other areas of mathematics or physics or even from finance. That would provide a sort of duality between a physically meaningful and involved problem, and a rather simple probability problem such as a box of balls of various colors.

Acknowledgments

I would like to thank A. Kryftis for keep challenging me and in particular, A. Anastasiou (LSE and Un. of Cyprus) for reading the manuscript and for providing insightful comments, and D. Christofides (UCLAN, Cyprus) for providing the solution of appendix F, and also S. Agapiou (Un. of Cyprus), S. Hormann (Graz Un. of Technology, Austria), T. Bruss (Universite Libre de Bruxelles), and E. Mossel and S. Sheffield (MIT) for referring me to useful sources. I also thank my friends at Vincent House (London), and especially G. Vogiatzi for encouraging me to publish this work, and my colleagues at JPMorgan Chase, and in particular, S. El Hamoui, A. Eriksson, M. Green, J. Lorenzen and S. Mcgarvie for engaging in stimulating discussions171717Disclaimer: The current paper consists of a personal work curried out by A. Taliotis, and it expresses his own personal views on the particular topic. In particular, this work is not, in any way, endorsed by or related with JPMorgan Chase (the firm) or its employees or stakeholders or any interests or entities represented by the firm.. A warm thank you to Staxto for “being here”, during the long nights after the JPMC office hours, while preparing this write-up. I would like to express my deep gratitude to my teacher and good friend Y. Kovchegov from The Ohio State University for patiently showing me, among other, how to carry on long and tedious calculations and how to, in a quantum field theoretical language, “renormalize” (deal with) infinities, and eventually obtain finite and meaningful results. Lastly, I would like to thank Nikolas and Leo-Anastasis for invading into my routine teaching me the beauty in life, and for all the time I have taken away from them.

Appendix A Gamma function expansions near poles

For Gamma function arguments $\in 0\cup\mathbb{Z}_{-}$, the following key equation provides the residues

[TABLE]

Using (A1), the following expansions, which are needed in the intermediate steps, can be derived.

[TABLE]

Appendix B Alternative boundary conditions, Problem 1, Rule III

The fact that $Y^{(0)}$ fulfills the boundary condition (48) for $k=0$ can be verified through a simple substitution. In order to see that the boundary condition for $Y^{(k>0)}$ is also satisfied one needs to note that181818We define $t\equiv-\sqrt{x_{1}}-\frac{1}{\sqrt{x_{1}}}$ and study the limit of $Y^{(k>0)}$ as $x_{2}\to t$. the overall multiplicative factor of (47b) decays as $(x_{2}-t)^{1+2k}$. Expanding the first term of the curly bracket (product of three factors) shows that this quantity grows as $(x_{2}-t)^{-k-1}$. Next, we move to the second term in the curly bracket, which is the term involving the hypergeometric function, and whose argument grows as $(x_{2}-t)^{-1}$. Taking into account that $k$ is a positive integer, we deduce that the ${}_{2}F_{1}\left(-1-2k,-k;-2k,x\right)$ is a terminating polynomial of degree $k$ and hence, we conclude that this term grows as $(x_{2}-t)^{-k}$. Thus, the curly bracket grows as $(x_{2}-t)^{-k-1}$ while the overall coefficient decays as $(x_{2}-t)^{1+2k}$ and thus the whole equation decays as $(x_{2}-t)^{k}$. This implies that the boundary condition is fulfilled $\forall\,k>0$.

We thus conclude that the solutions (47) could had been derived from the same differential equation (46a) but with the alternative boundary conditions (48).

Appendix C A second derivation of (44)

The current appendix describes briefly a second derivation of equation (44). Starting from (41) and performing the summation using standard Gamma function identities we obtain

[TABLE]

Using (A2c) and (A2d) for respectively the last two Gamma function ratios, (C1) becomes

[TABLE]

where in arriving in the second equality we used the identity $\Gamma(2z)=\frac{1}{\sqrt{2\pi}}2^{2z-\frac{1}{2}}\Gamma(z)\Gamma(z+\frac{1}{2})$ and dropped the $O\left(\frac{y-w}{x-r},\frac{x-r}{y-w}\right)$ terms assuming some regulation argument. Equation (C2) then can be written in terms of factorials and when it is simplified, it becomes exactly identical to (44). It is notable that unless we took sub-leading corrections in the expansions as in (A2c) and (A2d), we would be off by a factor of 2 with respect to (44).

Appendix D Finite summation identities involving a product of binomials and a rational function

We show that for positive integers $r$ and $w$

[TABLE]

The left hand side of (D1) for integer parameters is

[TABLE]

Next, we use the summation identity

[TABLE]

which implies

[TABLE]

Inserting (D4) in (D2) with $s=m=r+w$ and noting that the second binomial in the most right hand side of (D4) vanishes, we find

[TABLE]

and this completes the proof.

Appendix E Proving the normalization of the probability formula, Rule III

Starting from (75) we prove that the normalization condition (76) applies. Equation (75) yields

[TABLE]

where the summation is, up to overall Gamma function factors, a hypergeometric function ${}_{3}F_{2}$. The next step is to use the identity [3]

[TABLE]

for $a=b=1$, $n=w$, $d=2+r$, which yield $a+b-d-n+1=1-r-w$, and where $x_{n}\equiv\Gamma(x+n)/\Gamma(x)$ is the Pochhammer’s Symbol. This means that the identity is applicable for equation (E1), which after simplifying yields

[TABLE]

completing the proof.

Appendix F Solution to Problem 1, Rule III: a combinatorial-inductive derivation

In this appendix, a combinatorial-inductive derivation of (75), which is shorter and more intuitive, is given. On the other hand, this derivation is less systematic and it exclusively refers to the particular problem (Problem 1). In particular, it may not be generalized and adapted to other similar problems. An analogous derivation, as the one presented here using a similar logic, exists for Problem 1, Rule IV and its presentation is therefore omitted.

We consider sequences containing $r$ red ball removals out of $r$ and $w-k$ white ball removals out of $w$, ending with a red ball removal. Each removal occurs with probability $r^{2}/(r+w)^{2}$ for a red and with probability $w(2r+w)/(r+w)^{2}$ for a white ball respectively, where $r$ and $w$ are the number of the red and of the white balls at the time of the given removal. Defining with

[TABLE]

the probability of a given sequence of removals has the generic form

[TABLE]

It could happen that sometimes the sub-product of terms between successive white ball removals, i.e. the terms between $\frac{(w-i+1)(2r_{i}+(w-i+1))}{(r_{i}+(w-i+1))^{2}}$ and $\frac{(w-i)(2r_{i+1}+(w-i))}{(r_{i+1}+(w-i))^{2}}$, $i=1,\,2,\,...,\,w-k-1$, including the first sub-product of terms until the $\frac{w(2r_{1}+w)}{(r_{1}+w)^{2}}$ factor, could collapse to unity. Equation (F2) can be written as a product of two factors according to

[TABLE]

where the first factor is the term in the square bracket, which simplifies to $\frac{(r!)^{2}k!w!}{((r+w)!)^{2}}$. Interestingly, all the $r_{i}$ dependance drops out from this first term and hence, this factor becomes a multiplicative overall factor to the sequence (F2).

Using (F1) and (F3), the required probability is then given by

[TABLE]

where we have set $w=k+s$ and where we have defined

[TABLE]

Starting from the definition of $g^{(k)}(r,s)$, equation (F5), it is straightforward to show

[TABLE]

So the problem reduces in computing $g^{(k)}(r,s=w-k)$ from (F5) subject to the boundary cases (F6). We proceed by taking the difference $g^{(k)}(r,s)-g^{(k)}(r-1,s)$ from where we find

[TABLE]

where in the second equality we have assumed a renaming of the dummy $r_{i}$ index according to $r_{i}\rightarrow r_{i-1}$, $i=2,\,3\,...,\,s$. Thus we have arrived at the recursion191919If we had used the alternative ansatz $p_{III}^{(k)}(r,w)=\frac{(r!)^{2}k!w!}{((r+w)!)^{2}}f^{(k)}(r,w)$ instated of the ansatz of (21), equation (20) would had yielded the much simpler difference equation $f^{(k)}(r,w)=f^{(k)}(r-1,w)+(2r+w)f^{(k)}(r,w-1)$ rather than the recursion provided by equation (22).

[TABLE]

subject to equations (F6). In the remaining part we will only sketch the derivation rather than provide the details. While we could proceed by induction, we follow a more straightforward approach.

Using (F6) and (F8) for $s=2$ it is deduced that

[TABLE]

from where it is guessed that

[TABLE]

Equation (F10) can then be substituted in (F8), and with the help of (F6), one can verify that (F8) is indeed fulfilled. The next step is to work for $s=3$ using (F6) and (F10) to obtain

[TABLE]

from where it is guessed that

[TABLE]

Equation (F12) can then be verified through a substitution in (F8) with the help of (F10).

Observing then equations (F10) and (F12) we make the general ansatz solution

[TABLE]

which we verify by direct substitution in the recursion (F8).

Finally, using (F4) and (F13) we obtain

[TABLE]

reproducing equation (75), and this completes the alternative derivation of $p_{III}^{(k)}(r,w)$.

Remark: We could had started from (22) and subsequently worked along the same lines as between equations (F8)-(F14). Conversely, starting from (F8), we could had worked along the lines between equation (22) and equation (75); that is using the PDE approach. Either approach has its advantages: using the method of this appendix we reached the result faster. On the other hand, the approach here has involved educated guessing, which (in principle) has not been the case during the initial derivation presented by the main body of this paper. That approach has been purely systematic and it may, in principle, be generalized and be applied to other cases, as that approach is to a large extend problem independent.

Appendix G Problem 4: A challenging problem from the literature

The problem presented in [2] is as follows.

Problem 4

A bag contains $m>0$ black balls and $n>0$ white balls. A sequence of balls from the bag is discarded in the following manner:

(i) A ball is chosen at random and discarded. (ii) Another ball is chosen at random from the remainder. If its color is different from the last it is replaced in the bag and the process repeated from the beginning (i.e. (i)). If the second ball is the same color as the first it is discarded and we proceed from (ii). Thus the balls are sampled and discarded until a change in color occurs, at which point the last ball is replaced and the process starts afresh.

*The question is: what is the probability that the final ball should be black? *

For the record, the authors show that the requested probability is as simple as $1/2$.

Appendix H Python simulation code and some results

The python codes that simulate the removals of Problem 1 are provided below.

The first block of the code that follows (Problem 1, Rule III) imports several necessary Python packages. The “newListAfterRemoval” function provides the updated content of the box once a ball is removed. The function “redWhite_Problem1_Rule_IIII($r,w$,epochs)” refers to Problem 1, Rule III, and it takes care of everything else. The arguments $r$ and $w$ refer to the initial red ($r$) and the initial white ($w$) balls respectively. The variable “epochs” is a positive integer and defines the number of times that the experiment will be repeated. In particular, this function provides the simulated probabilities (Rule III) contrasting them with the theoretical ones. This part of the code is supplemented by explanations (see comments in green fonds).

Likewise, the function “redWhite_Problem1_Rule_IV($r,w$,epochs)” refers to Problem 1, Rule IV, and it takes care of everything else. The arguments in the function are analogous to those of Rule III. In particular, this function provides the simulated probabilities (Rule IV) contrasting them with the theoretical ones. This part of the code is not supplemented by explanations as these are analogous to those of “redWhite_Problem1_Rule_IIII($r,w$,epochs)”.

See pages 1 of Rule_III_code2.pdf

Below we present the simulated and the theoretical results for a few cases for Problem 1, Rule III.

Next, we present the code and the results for Problem 1, Rule IV.

See pages 1 of Rule_IV_code2.pdf

Below we present the simulated and the theoretical results for a few cases for Problem 1, Rule IV.