ANALYSIS OF A CLASS OF CANNIBAL URNS

In this note we study a class of 2 × 2 Pólya-Eggenberger urn models, with ball transition matrix given by M = (cid:0) 0 − b c − d (cid:1) , for arbitrary parameters b , c ∈ N , assuming the balancing condition d = c + b . This urn model serves as a stochastic model in biology describing cannibalistic behavior of populations. The special case b = c = 1 and d = 2 was studied by Pittel [ 16 ] using asymptotic approximation techniques, and more recently by Hwang et al. [ 10 ] using generating functions. We obtain limit laws for the stated class of cannibal urns by using Pittel’s method, and also different techniques.


Urn models with a diminishing content
In this work we study Pólya-Eggenberger urn models with ball transition matrix M given by for arbitrary parameters b, c ∈ N = {1, 2, . . . }, assuming that d = b + c. The dynamics of this urn model can be described as follows. At the beginning the urn contains n white and m black balls. At every discrete time step, we draw a ball at random from the urn, examine its color and put it back into the urn and then add/remove balls according to the observed color by the following rules. If the ball is white, then we remove b black balls from the urn, while if the ball is black, then c white balls are added, and d black balls are removed from the urn. The condition d = b + c implies that the urn is balanced, at every step the total number of balls in the urn decreases by exactly b balls regardless of the color of the drawn ball. We are interested in the distribution of the random variable X n,m , counting the number of white balls, when the process stops, i.e. when the number of black balls is below d, starting with n white and m black balls, respectively, with n ≥ and m ≥ d. Apparently, the process stops after a finite number of steps, and the number of black balls is diminishing throughout the evolution of the urn. This urn model is part of a general setting for Pólya-Eggenberger urns with a diminishing nature, described in detail in [10], see also [9]. For ordinary tenable urn models, where the process of drawing and adding/removing balls can be repeated ad infinitum we refer the reader to [13,14,15], and also for some recent developments to [2,7,5,6,11,12,20,21]. A special instance of the class of urns with ball transition matrix given by (1) is the so-called cannibal urn model, b = c = 1 and d = 2. This urn model has been introduced by R. F. Greene, see [16], as a stochastic model in biology modelling cannibalistic behavior in a population. It can be described as follows: a population consists of non-cannibals and cannibals. At each time instant a non-cannibal is selected as a victim and consumed/removed and after that a member of the remaining population (non-cannibals and cannibals) is selected at random. If the selected individual is a cannibal it remains as a cannibal in the population, but if the selected individual is a non-cannibal it changes its behavior and will become a cannibal. The question is, when starting with n cannibals and m non-cannibals, what is the number of cannibals in the population when all non-cannibals are removed? In the general case we assume that b individuals are consumed at each discrete time step, and that c individuals may become cannibals. The general model improves upon the original model insofar as it takes into account the effects of consuming more than one individual at each time step and how it effects the overall number of cannibals. The original urn model, b = c = 1 and d = 2, has been analyzed by Pittel [16].
His approach is based on asymptotic approximation; more precisely it consists of an asymptotic approximation of the recurrence of the moment generating function of the random variable of interest; we refer to the works of Pittel [17,18,19] concerning the numerous applications of this versatile technique. Pittel obtained a normal limit law for X n,m in a certain growth region of n and m, namely the fraction ρ n,m = n/(n + m) has to be bounded away from one, as N n,m = n + m tends to infinity. Recently, this model was treated by Hwang et al. [10] using generating functions, obtaining limit laws for X n,m in every growth region of n and m. Unfortunately, the generating functions approach of [10] is restricted to the special case b = c = 1 and d = 2. We give a brief discussion of the difficulties of a generating functions approach in the last section of this work.
We will use the approach of Pittel [16] to obtain asymptotic normality of X n,m in the general model, assuming that the fraction ρ n,m = n/(n + m) is bounded away from one, as N n,m = n + m tends to infinity. Furthermore, we can also handle some other cases, where ρ n,m → 1 using different techniques. More precisely, assuming that n → ∞ we prove that the distribution of X n,m is asymptotically normal distributed for n m n 2/3 , is Poisson distributed for m ∼ n, and degenerates for m n. Unfortunately, we could not determine the limiting distribution of X n,m in the region n 2/3 ≤ m n.

Urn models and weighted lattice paths
It is well known that the evolution of urn models can be described by weighted lattice paths, see e.g. Flajolet et al. [5]. Concerning our particular model (1) we have the following description. If the urn contains n white balls and m black balls and we select a white ball (with probability n m+n ), then this corresponds to a step from (m, n) to (m − b, n), to which the weight n m+n is associated; and if we select a black ball (with probability m m+n ), this corresponds to a step from (m, n) to (m − d, n + c) (with weight m m+n ). The weight of a path after t successive draws consists of the product of the weight of every step. By this correspondence, the probability of starting at (m, n) and ending at ( j, k) is equal to the sum of the weights of all possible paths starting at state (m, n) and ending at ( j, k). The expressions so obtained for the probability are, although exact, often less useful for large m or n. Note that, strictly speaking, the description of cannibalistic behavior described above does not fall directly in the scheme of urn models specified by (1); it is a slight modification of the urn model, where the weight of a step from (m, n) to (m − b, n) (and thus the probability) is given by n n+m−b and the weight of a step from (m, n) to (m − d, n + c) is given by m−b m−b+n (instead of weights n m+n and m m+n ). However, both types of weights lead to the same asymptotic behavior, and the proof is readily adapted. Therefore we opted to use the standard Pólya-Eggenberger urn model weights.

Goal
Our main goal is to establish the asymptotic normality of X n,m as n → ∞ and n/(n+ m) is bounded away from one (which excludes the case m = o(n)). Second, we show that assuming that n → ∞ and ρ n,m → 1 the distribution of X n,m is still asymptotically normal for intermediate m such that n m n 2/3 , is Poisson distributed for small m such that m ∼ n, and degenerates for very small m such that m = o( n). The proof technique of the asymptotic normality applied here in the case of N n,m → ∞, with ρ n,m = n/(n + m) bounded away from one, follows closely the asymptotic approximation proof method of Pittel of the original problem [16].

Notations and Overview
We use the notations N = {1, 2, . . . } and N 0 = {0, 1, . . . }, and denote with − → the convergence in distribution. We use the abbreviations N n,m = n + m and ρ n,m = n/N n,m = n/(n + m), and the notations m n, standing for m = o(n), n ∼ m standing for m/n → 1, as n tends to infinity. For sequences (n, m) = (n i , m i ) i∈N , by the formulation "ρ m,n is bounded away from one" we mean that for any > 0 the inequality n i /(n i + m i ) > 1 − ε holds for only finitely many i. Throughout this work we mostly drop the subscript i for the sake of simplicity. The notation n m n 2/3 for n → ∞ means that we consider double sequences (n, m) = (n i , m i ) i∈N with n i → ∞, such that both conditions m = o(n 2/3 ) and n = o(m) are satisfied; this implies that m tends to infinity as n → ∞. In the next section we state the limiting distributions of X n,m , specified according to the growth of n and m. In Section 3 we establish the part concerning the normal limit of X n,m for ρ n,m bounded away from one. Section 4 contains the derivations of the limiting distribution results for intermediate and small m. The last section is devoted to a brief discussion of an alternative approach based on generating functions.

Results
In the following we will state the limiting distribution of the random variable X n,m , counting the number of white balls when the number of black balls is below d, with respect to the general Cannibal urn model with ball transition matrix (1), starting with n white and m black balls. here φ(t) and ψ(t) are given by and

For small m such that m ∼ λn
).

For very small m such that m n the distribution of X
The result stated in the first case of Theorem 1 can be extended as follows.

Theorem 2.
For N n,m = n + m → ∞ with ρ n,m = n/(n + m) being bounded away from one with φ(t) and ψ(t) as given in Theorem 1.
This includes in particular the cases n → ∞, with m at least of order n, and the cases of fixed n ∈ N and m tending to infinity.

Remark 1.
We did not manage to obtain results for the missing region n 2/3 ≤ m n. We expect that the centered and normalized random variable X n,m remains normal distributed in the limit n → ∞, by analogy with the results of [10], where the limiting distributions of the original problem (b = c = 1, d = 2) were obtained using generating functions.

Remark 2.
As observed by Pittel in the special case b = c = 1 and d = 2, we anticipate that E(X n,m ) is for large n compared to m close to n and the variance is close to zero.

Asymptotic normality by asymptotic approximation
In the following we will prove Theorem 2, which includes Theorem 1 case (1), using Pittel's approach [16] to the original problem b = c = 1, d = 2.

Recurrence relations
By definition of Pólya-Eggenberger urn models the random variable X n,m satisfies the following recursive distributional equation, which is obtained by considering the urn after the first draw when starting with n white and m black balls.
where I n,m denotes the indicator variable of the event that a white ball is drawn, with I n,m being independent of (X n,m ) n,m≥0 and (X n,m ) n,m≥0 , which are independent copies of each other. By (2), the Laplace transforms g n,m (z) = E(e zX n,m ) of X n,m satisfy a similar recurrence relation, given as follows.
When b = c one may start with b · n white and b · m black balls, reducing the problem essentially to the original cannibal urn model. In order to prove asymptotic normality of X n,m (suitably centered and normalized) using asymptotic approximation techniques, see [16,17,18,19], we introduce an auxiliary sequence of random variables (Y n,m ) n,m≥0 , whose Laplace transforms h n,m (z) = E(e zY n,m ) are given by with functions φ(t) and ψ(t) as given in Theorem 1. Clearly, the random variables Y n,m are normal distributed with mean N n,m φ(ρ n,m ) and variance N n,m ψ(ρ n,m ). We will prove that h n,m (z) = E(e zY n,m ) almost satisfies a recurrence (3), see Lemma 1. The reason for this, as we shall see, is that φ(t), ψ(t) is the solution of a system of differential equations (8) and (9). We will see that for N n,m = n + m → ∞ such that ρ n,m is bounded away from 1, ψ(ρ n,m ) is bounded away from 0, see Lemma 2, and that the Laplace transforms of X n,m and Y n,m are close enough to imply, see Lemma 5, that X n,m is asymptotically normal distributed with mean N n,m φ(ρ n,m ) and variance N n,m ψ(ρ n,m ). Note that the assumption N n,m = n + m → ∞ such that ρ n,m is bounded away from one includes cases n → ∞ and m at least of order n, but also other cases, i.e. arbitrary but fixed n ∈ N and m → ∞.

An estimate for the Laplace transforms
Lemma 1. For fixed u > 0 and z ∈ R such that N 1/2 n,m |z| ≤ u for N n,m → ∞ the Laplace transforms h n,m (z) = E(e zY n,m ), as defined in (4), according to the φ(t) and ψ(t) given in Theorem 1, satisfy , uniformly over d + 1 ≤ l ≤ N n,m and 0 ≤ k ≤ N n,m − l.

Remark 4.
The functions φ(t) and ψ(t), appearing in the definition (4) of h n,m (z) = E(e zY n,m ), are actually a priori unknown; they are determined via the differential equations (8) and (9) such that the statement of Lemma 1 is true. However, in order to present the proof of the lemma above in a clear manner it is beneficial to reverse the order of argumentation, making the proof more transparent. The condition N 1/2 n,m |z| ≤ u is natural with respect to the fact that we are interested, see Lemma 5, in the normalized quantities h n,m (z/N 1/2 n,m ) and g n,m (z/N 1/2 n,m ).
Proof. In order to obtain the stated asymptotic expansions of k (k+l) h k,l−b (z) and l (k+l) h k+c,l−d (z) we note first that We apply Taylor where the error term is uniform in both k and l due to our previous considerations. Subsequently, we use the auxiliary quantities in order to write (6) as We obtain by expansion at z = 0 The error term, obtained by third order Taylor approximation of exp(z 2 ) and second order approximation of exp(z), is again uniform in both k and l. Note again that by definition 0 ≤ ρ k,l ≤ 1, for all k and l; moreover, we get uniform upper bounds for the quantities A k,l , B k,l , C k,l and D k,l using upper bounds for the absolute values of the functions φ(t), ψ(t) and their derivatives on the unit interval. Furthermore, we use the fact that |z| 3 = ( |z| N k,l ). The coefficient of z is zero, since it is readily checked that the function φ(t), as given in Theorem 1, satisfies the differential equation stated below on the interval [0, 1] with initial condition φ(1) = 1.
Moreover, the coefficient of z 2 also equals zero, since the function ψ(t), as given in Theorem 1, satisfies the differential equation stated below on the interval [0, 1] with initial condition ψ(1) = 0 with inhomogeneous part F (t) given by Consequently, we get the result Next we study the roots of the function ψ(t) in order to ensure that the random variable X * n,m , defined by X * n,m = , is well defined for all n, m ≥ 1.
Note that many computer algebra systems deliver a wrong solution for ψ(t). We remark the special solutions for the cases b = c, and b = 2c, can be obtained from the general solution ψ(t) for b = c, 2c by taking the appropriate limits. An immediate application of Lemma 2 is the following result concerning the value ψ(ρ n,m ).

Lemma 3.
For N n,m = n + m → ∞, such that ρ n,m = n/(n + m) is bounded away from 1, the value ψ(ρ n,m ) is bounded away from 0, with function ψ(t) as given in Theorem 1.

Relating the Laplace transforms
We have proven in Lemma 1 that the Laplace transforms h n,m (z), as defined in (4), of the normal distributed random variables Y n,m almost satisfy the recurrence relation 3 with respect to an error term as given in Lemma (1). The next result relates the Laplace transforms g n,m (z) = E(e zX n,m ) and h n,m = E(e zY n,m ) of the random variables X n,m and Y n,m .

Lemma 4.
For fixed u > 0 and z ∈ R such that N 1/2 n,m |z| ≤ u, we have Proof by induction. We proceed exactly along the lines of Pittel [16] by using induction with respect to the total number of balls N n,m contained in the urn. By definition of φ(t) and ψ(t), see Theorem 1, we have g n,0 (z) = h n,0 (z) for all n ∈ N. There exists a constant λ > 0 such that for all N n,m < d g n,m (z) = exp(nz) ≤ h n,m (z) exp λ|z| .
The starting point for the induction is N 0 , with 1 ≤ N 0 < d. We assume now that the constant λ > 0 is large enough such that , for all k ∈ N and l ∈ N with k + l = N 0 + i b, with 0 ≤ i < j and j = j(N n,m ) with respect to N n,m = j b + N 0 . Then we have for k + l = N n,m = j b + N 0 The other direction is proven as follows. There exists a constant λ > 0 such that for all k + l < d h k,l (z) ≤ g k,l (z) exp λ |z| .
Choosing the same N 0 as before, with 1 ≤ N 0 < d, we assume now that the constant λ > 0 is large enough such that and also according to Lemma 1 for all k ∈ N and l ∈ N with k + l = N 0 + i b, with 0 ≤ i < j. Then we have for k + l = N n,m We can write where Ψ(t) = d d t log Γ(t) denotes the Digamma function, see [1]. Using the asymptotic expansion of the Digamma function Ψ(n + 1) = log n + 1/n + (1/n 2 ) and setting N n,m = N 0 + j b we get ).

Corollary 1.
The random variable X * n,m is asymptotically normal distributed.
Proof. By Curtiss' theorem [4] the pointwise convergence of the sequence of moment generating functions implies the weak convergence of the corresponding sequence of random variables. Hence, by Lemma 5 we obtain the asymptotic normality of the random variable X * n,m .

Proofs for the other regions
In order to prove the results in the cases two to four of Theorem 1, with n tending to infinity and ρ n,m → 1, we proceed differently by studying directly the probabilities P{X n,m = k}. The distributional equation (2) implies the following recurrence relation for the probabilities P{X n,m = k}.
where δ n,k denotes the Kronecker delta function.

The degenerate case
We consider the fourth case of Theorem 1. In order to obtain a closed formula for P{X n,m = n} we use the recursive description of the probabilities P{X n,m = k} as stated in (11). The number of white balls is non-decreasing, in other words P{X n+ ·c,m = n} = 0, for all ≥ 1 and m ≥ 0. Hence, we get for k = n the simple recurrence relation This recurrence relation is readily solved and we get the following explicit expression for P{X n,m = n}, valid for m ≥ 0 and n ≥ 0. .
In order to obtain an asymptotic expansion of P{X n,m = n} = P{X n,m − n = 0} for n → ∞ and m n we will use Stirling's formula for the Gamma function and the following expansion of (1 + x n ) n with respect to n tending to infinity and x = x(n), x n 2/3 , Application of (12) and (13) gives for m = o( n) the asymptotic expansions By combining the stated asymptotic expansion we obtain the result ∼ 1, or equivalently P{X n,m − n = 0} ∼ 1.

The Poisson case
We study the third case of Theorem 1 and derive a closed formula for P{X n,m = n + ck}, for arbitrary k ≥ 0, from the recurrence relation 11 using combinatorial arguments with respect to the weighted lattice path description of the evolution of the urn model, see The k events that a black ball is drawn, which leads to an addition of c white balls and the removal of d black balls, occur during the . The upper bound is due to the fact that we still have to draw k − 1 black balls until the evolution of the urn stops, i.e. m < d. The second draw of a black ball can happen at stage 2 + 1 + 2 of the evolution of the urn, with 0 ≤ 2 ≤ m−kd b − 1 . Similar arguments provide the possible ranges for all of the k draws of black balls. Each of such event of drawing a black ball happens with probability proportional to the number of blacks balls in the urn. Hence, the numerator r(n, n, k) of P{X n,m = n + ck} is given m n 2/3 . It seems difficult to extend this result beyond m ∼ n 2/3 , due to the expansion of the denominator s(n, m, k) with respect to 13. Moreover, the results of [10] suggest that for m = o(n) the expectation satisfies the expansion We will use the following well known result concerning the asymptotic normality of Poisson distributed random variables [3].

Lemma 6.
Assume that the sequence of random variables (X t ) t∈N is Poisson distributed with parameter Λ t , where Λ t → ∞ as t tends to infinity. Then X t , suitably centralized and normalized, is asymptotically normal distributed, x −∞ e − u 2 2 du denotes the standard normal distribution function.
By Lemma 6 and we have to show that (X n,m − n)/c is asymptotically Poisson distributed with parameter Λ = Λ(n, m) = m 2 /(2bn), for k = Λ+ x · Λ, with arbitrary but fixed x ∈ R. Then, an application of Lemma 6 proves the normal limit. Since we already know that the probabilities P{X n,m − n ≤ ck} = r(n, m, k)/s(n, m, k), with r(n, m, k) and s(n, m, k) given by (15) and (14), respectively, the remaining obstacle is to refine our previous results and to obtain a uniform bound on the error term in the expansion of the numerator r(n, m, k) (15) and the denominator s(n, m, k) (14), with respect to 1 ≤ k ≤ Λ + x · Λ and Λ = Λ(n, m) = m 2 /(2bn). First we observe that for all choices of 1 , . . . , k , uniformly for k ≤ Λ + x · Λ, and arbitrary but fixed x ∈ R. We obtain the uniform asymptotic expansion By application of Euler's formula we get uniformly for k ≤ Λ + x · Λ. Consequently, we get the required uniform expansion of the numerator r(n, m, k) (15). Using similar uniform expansion of the denominator s(n, m, k), as given in (14), and combining the two estimates, we finally obtain

Exact solvability -An outlook
As mentioned in the introduction Hwang et al. [10] obtained simple expressions for the exact probability distribution of X n,m in the case of the original cannibal urn model b = c = 1 and d = 2, which allows to give a complete description of the limiting distributions of X n,m . We obtain from (2) the following recurrence for the probability generating function f n,m (v) = E(v X n,m ) = k≥0 P{X n,m = k}v k of X n,m .
Concerning explicit solutions, such recurrences can be treated with ordinary (bivariate) generating functions leading to first order partial differential equations, see [9,10]. However, a big problem is that inhomogeneous part is not determined by the initial conditions, due to the terms f n+c,m−d (v).
The easiest way to overcome this deficit is to restrict the number of white balls to multiples of c and to study modified probability generating functionsf n,m (v) of X cn,m , defined byf n,m (v) = 1 n! f cn,m (v). We have already seen in the normal limit part of the main result, Theorem 1, that the case b = c = 1 and d = 2 is somewhat special. Indeed, for b = c = 1 and d = 2 one may use additional weights n+m n /(m + 1) in order to reduce the problem to a first order partial differential equation, whereas for b = c this is not possible.