- Email: [email protected]
Distributions in a class of Poissonized urns with an application to Apollonian networks
Statistics and Probability Letters 115 (2016) 1–7
Contents lists available at ScienceDirect
Statistics and Probability Letters journal homepage: www.elsevier.com/locate/stapro
Distributions in a class of Poissonized urns with an application to Apollonian networks Panpan Zhang ∗ , Hosam M. Mahmoud Department of Statistics, The George Washington University, Washington, DC 20052, USA
article
info
Article history: Received 17 October 2015 Received in revised form 23 March 2016 Accepted 25 March 2016 Available online 8 April 2016 MSC: 60G20 65Q30 35A99 Keywords: Pólya urn Pólya process Partial differential equation Stochastic recurrence
abstract We study a class of Pólya processes that underlie terminal nodes in a random Apollonian network. We calculate the exact first and second moments of the number of terminal nodes by solving ordinary differential equations. These equations are derived from the partial differential equation governing the process. In fact, the partial differential equation yields a stochastic hierarchy of moment equations, which can be bootstrapped to get higher moments from the equations that have been solved for lower moments. We also show that the number of terminal nodes, when appropriately scaled, converges in distribution to a gamma random variable via the method of moments. The asymptotic results can be obtained using classic methods of branching processes. The manuscript explores the potential of an alternative method capable of producing exact moments and rates of convergence. © 2016 Elsevier B.V. All rights reserved.
1. Introduction Urn schemes growing in discrete time have applications in many fields. Early applications include modeling epidemiological contagion (the Pólya-Eggenberger scheme) and modeling gas diffusion and heat transfer (the Ehrenfest model). Throughout the last century, many more applications were proposed, literally hundreds, and it is not possible to concisely list them all, but we direct the interested reader to the numerous applications discussed in Johnson and Kotz (1977) and Mahmoud (2008). Associated with a Pólya urn scheme is a process obtained by embedding in continuous time. Processes growing in continuous time are more realistic and applicable than those growing in discrete time. In practice, there is no guarantee that the epochs of process will be equispaced. They also address a different question—that of the status after a certain amount of (real) time, rather than the status after a certain number of events. For better readability, we call the evolution in continuous time the Pólya process (as in Balaji and Mahmoud, 2006; Sparks and Mahmoud, 2013), and reserve the term ‘‘scheme’’ for urns evolving in discrete time. The Pólya process was introduced in Athreya and Karlin (1968) (though it was not named that) as a mathematical transform (Poissonization) to shed light on urns evolving in discrete time. Athreya and Karlin (1968) reports that dePoissonization, to obtain the inverse transform, is highly nontrivial, and until today a rigorous inverse transform remains elusive.
∗
Corresponding author. Tel.: +1 3472858737. E-mail addresses: [email protected] (P. Zhang), [email protected] (H.M. Mahmoud).
http://dx.doi.org/10.1016/j.spl.2016.03.023 0167-7152/© 2016 Elsevier B.V. All rights reserved.
2
P. Zhang, H.M. Mahmoud / Statistics and Probability Letters 115 (2016) 1–7
2. The Pólya process A two-color Pólya process is an urn containing balls of up to two colors, say white and blue, and evolving over time according to certain rules. Balls are drawn at epochs in time with independent Exp(1) interarrival times.1 When a white ball is drawn, it is placed back in the urn together with a white balls and b blue balls; when a blue ball is drawn, it is placed back in the urn together with c white balls and d blue balls. In general these numbers can be positive, zero or negative (interpret adding a negative number of balls to mean taking that many balls out of the urn). An urn process is called tenable, if it is always possible to execute the rules; that is, there never comes a time when the urn becomes empty or the rules call for removing a number of balls of a color more than their number in the urn. It is customary to represent these dynamics by the replacement matrix
a c
b . d
The rows of this matrix are indexed with white and blue, respectively from top to bottom; the columns are indexed with white and blue, respectively from left to right. To be precise about the dynamics of evolution in continuous time, we add a few words of explanation. Each ball in the process is endowed with a clock that rings in exponential time Exp(1). All the clocks are independent of each other and of any other random variables related to the past. At a renewal point (an epoch), we execute the rules associated with its color instantaneously. All new balls added come endowed with their own independent clocks. By the memorylessness of the exponential interarrival time, the Pólya process is Markovian. However, in general, it is an inhomogeneous jump process, as the rates keep changing owing to ball additions. We are interested in the class of two-color Pólya processes with the replacement matrix
0 k
k−1 , −1
(1)
for integer k ≥ 3. The study is motivated by an application from the growth in real time of Apollonian networks. We call this class of Pólya processes the Apollonian process. 3. Organization The rest of this paper is organized as follows. In Section 4, we introduce a discrete-time Apollonian network. A coloring scheme of its cliques connects it to Pólya urns. The Poissonized counterpart is discussed in Section 5. In Section 5.1, we take up the mean and variance of the number of terminal nodes at time t by finding the joint moment of white and blue cliques defined in Section 4. Higher asymptotic moments are derived in Section 5.2, leading to a conclusion that the scaled number of terminal nodes at time t asymptotically has a Gamma distribution. Section 6 concludes the paper with remarks comparing the behavior of terminal nodes in both discrete- and continuous-time Apollonian networks. We also mention connections to branching processes that can alternatively be utilized to get the asymptotics. 4. Apollonian networks Apollonian networks are graphs based on triangulation. More generally, Apollonian networks are constructed on subdivision of simplexes of higher order. As random graphs, Apollonian networks have the following features: scale free, small world, Euclidean, space filling and with matching graphs. Each of these features has a significant impact in current random network research. We refer the reader to Andrade et al. (2005) for motivation and important applications of Apollonian networks. For purposes of initial introduction, we shall use triangulation, but our intent is to study general simplexes. Consider an active cycle of length 3 (with three edges, and three vertices labeled with three zeros to indicate they are the initial nodes). The term ‘‘active’’ indicates that it is a ‘‘recruiter’’ of other self-similar graphs. A random Apollonian network (of index 3) is a random graph that grows as follows. At the first step, we insert a new vertex (labeled with 1) at an arbitrary position inside the original cycle, connect it with three edges to the three initial nodes (labeled with 0), and we then discard the recruiter. We thus have three new active cycles after the first step of insertion, and the initial cycle becomes inactive and it does not recruit any more. In the second and subsequent steps, we choose an active cycle at random (all active cycles being equally likely), and insert a new vertex inside the cycle. We then connect the newly inserted vertex with three edges to the three nodes of the chosen cycle, creating three new active cycles. Lastly, we deactivate the recruiting cycle. At step n ≥ 1, the newly added vertex is labeled with n. For a general index k ≥ 3, the formal definition is the following. A random Apollonian network of index k is a random graph that grows out of a k-clique with vertices labeled with zeros. At the nth step, a k-clique is chosen at random. We insert
1 Exp(1) refers to an exponential random variable with parameter 1.
P. Zhang, H.M. Mahmoud / Statistics and Probability Letters 115 (2016) 1–7
3
a new vertex in the graph, connect it with k edges to the k vertices of the chosen clique, and we then discard the recruiting clique. For more motivational information about Apollonian networks, we refer the reader to Zhang and Mahmoud (2016). In various evolutionary random graphs, the number of terminal nodes—nodes that appeared but did not participate yet in further recruiting—is a significant property. For instance, leaves are the terminal nodes in random trees, and they were extensively studied in many tree families; see for example Janson (2005). In an Apollonian network, we consider a node terminal, if all the cliques that are incident with this node have not yet recruited any other cliques. (k) Let Xn be the number of terminal nodes in an Apollonian network with index k after n node insertions. We devise a (k) two-color Pólya urn strategy to determine the distribution of Xn . Let us color cliques incident with terminal nodes with white (W ), and all other cliques with blue (B). Consider the network at insertion step 0 to be the initial condition of an urn, having no white cliques, and one blue. Note that no two terminal nodes ever appear in the same clique. So, a count of white cliques (Wn ) at insertion step n translates directly into a count of terminal nodes, namely Xn(k) =
Wn
. (2) k When a white clique recruits, the terminal node in it is no longer terminal, and the recruiter turns inactive. The other k − 1 white cliques incident with this chosen terminal node turn into blue. However, all k newly born cliques are white; on net there is no change in the number of white cliques. On the other hand, if a blue clique recruits, we lose one blue clique as the recruiter becomes inactive; meanwhile, k new active white cliques appear. Therefore, we have a corresponding Pólya urn scheme with the replacement matrix as given in (1). This is an instance of the well-studied Bagchi–Pal urn (see Bagchi and Pal, 1985), and we can conduct a probabilistic analysis based on a known urn theory (see Athreya and Karlin, 1968; Bagchi and Pal, 1985). (k)
Proposition 4.1 (Zhang and Mahmoud, 2016). Let Xn age n. We then have (k)
Xn
n
a.s.
−→
k−1 2k − 1
be the number terminal nodes in an Apollonian network of index k at
,
and (k) k−1 Xn − 2k n −1
√
n
D
−→ N
0,
k(k − 1)2
.
(2k − 1)2 (3k − 1)
5. Poissonized Apollonian networks We next consider the analogue of Apollonian networks obtained by embedding processes in continuous time, which may be a more realistic reflection of applications. In other words, we are Poissonizing the urn with replacement matrix (1). In an Apollonian process generated from an Apollonian network of index k, let X (k) (t ) be the number of terminal nodes at time t, and W (t ) and B(t ) be the number of respectively white and blue cliques based on the color code we introduced in the previous section. The Apollonian process is the two-component vector
φ(t , u, v) := E e
W (t ) B(t )
. Let
uW (t )+v B(t )
be the joint moment generating function of
W (t ) B(t )
. For a general tenable two-color Pólya process, a certain partial
differential equation (PDE) is known; see Balaji and Mahmoud (2006). It is derived for a Pólya process with a deterministic replacement matrix. The equation is generalized in Mahmoud (2008) to Pólya processes with random replacement matrices. For deterministic replacement matrices, the PDE is
∂φ(t , u, v) ∂φ(t , u, v) ∂φ(t , u, v) + 1 − eau+bv + 1 − ecu+dv = 0. ∂t ∂u ∂v This functional equation has only been solved in a very limited number of cases, as for instance the case of forward and backward diagonal processes (see Balaji and Mahmoud, 2006), and the Ehrenfest process (see Balaji et al., 2006). The quest continues to find more classes of Pólya processes amenable to analysis via this method. As we shall see, the Apollonian process is one such case. Specialized to the number of terminal nodes of an Apollonian network, the PDE is
∂φ(t , u, v) ∂φ(t , u, v) ∂φ(t , u, v) + 1 − e(k−1)v + 1 − eku−v = 0. ∂t ∂u ∂v
(3)
Toward ordinary differential equations for moments, for any nonnegative integers i and j, we take the partial derivative i times for u and j times for v , respectively, for each term in (3). Setting u = v = 0, we obtain a recurrence for the mixed joint
4
P. Zhang, H.M. Mahmoud / Statistics and Probability Letters 115 (2016) 1–7
moments E W i (t )Bj (t ) as follows:
d
j−1 j E W (t )B (t ) = (k − 1)j−s E W i+1 (t )Bs (t )
dt
i
j
s
s =0
j−1 i i−1 r i j i i−r r j−s i−r s +1 + (−1) k E W (t )B (t ) + k E W (t )Bj+1 (t ) .
r
r =0 s=0
s
r =0
(4)
r
5.1. The mean and variance For small values of i and j, We are able to determine the exact moments of X (k) (t ) by solving the ordinary differential recurrence equations for W (t ). We illustrate this paradigm with the mean and variance where we have 1 ≤ i + j ≤ 2. Proposition 5.1. Let X (k) (t ) be the number of terminal nodes in a Poissonized Apollonian network of index k at time t. We have
E X (k) (t ) =
1
Var X (k) (t ) =
2k − 1 k−1
(2k − 1) −
1
e(k−1)t −
2
2k − 1
e−kt , k−1
e2(k−1)t −
1
(2k − 1)(3k − 1) k2 + 2k − 1 − e−2kt . (2k − 1)2 (3k − 1)
e−kt
(2k − 1)
e(k−1)t +
2k
(2k − 1)2
e−t
Proof. We determine the expectation and variance of X (k) (t ) via the expectation and variance of W (t ). We compute the first moment of W (t ) by a theorem derived in Mahmoud (2008) as follows:
0 E W (t ) =e k E B(t )
⊤
k−1 −1
W (0) , B(0)
t
where ⊤ is the transpose of a matrix, and W (0) and B(0) are initial numbers of white and blue cliques, respectively. We thus have
0 E W (t ) k−1 =e E B(t )
k −1
t
0 . 1
Reading off the components of this vector one by one, we have k k e(k−1)t − 2k − 1 2k − 1 k − 1 (k−1)t k E B(t ) = e + 2k − 1 2k − 1
E W (t ) =
e−kt , e−kt .
We resort to the established recurrence (4) for the joint moments of W (t ) and B(t ) to compute the second moment of W (t ). We simplify the recurrence for the three cases {i = 2, j = 0}, {i = 1, j = 1}, and {i = 0, j = 2}, and obtain the following matrix form
2 0 E W (t ) E W (t )B(t ) = k − 1 dt 0 E B2 (t )
d
2k −1 2(k − 1)
2 k2 E B(t ) 0 E W (t ) k E W (t )B(t ) + −k E B(t ) . −2 E B2 (t ) (k − 1)2 E W (t ) + E B(t )
This is a nonhomogeneous differential linear equation for
E[W 2 (t )] E[W (t )B(t )] E[B2 (t )]
. Notice that the coefficient matrix
is invertible. Hence, for all k ≥ 3, we have a general unique solution to this type of linear equation:2 k t
2 E W (t ) E W (t )B(t ) = e 0 E B2 (t )
0
−1 0
2 See Edwards and Penney (2008) for details.
2k −1 2k − 1
0 k (t −s) −2
k2 E B(s) −k E B(s) ds 2 (k − 1) E W (s) + E B(s)
0 2k 0 k − 1 −1 k 0 2k − 1 −2
P. Zhang, H.M. Mahmoud / Statistics and Probability Letters 115 (2016) 1–7
0
k
−1
+e
2k
0 k t E W 2 (0) −2 E W (0)B(0) , E B2 (0)
−1 2k − 1
0
where the initial vector condition is
5
E[W 2 (0)] E[W (0)B(0)] E[B2 (0)]
0 =
0 1
. We plug in the results of E W (t ) and E B(t ) , compute the
matrix exponential and the integration, and get the following simplified solutions for E W 2 (t ) , E W (t )B(t ) and E B2 (t ) . We obtain
k3
E W 2 (t ) =
k2 (k − 1)
e2(k−1)t −
(2k − 1)
e(k−1)t +
2k2 (k − 1) −t e (2k − 1)2
(2k − 1)(3k − 1) k3 (k − 1) − e−kt − e−2kt , 2k − 1 (2k − 1)2 (3k − 1) k2 (k − 1) 2(k−1)t k(k − 1) (k−1)t k(k − 1) −t k3 (k − 1) − − + e e e e−2kt , E W (t )B(t ) = (2k − 1)2 3k − 1 (2k − 1)2 (2k − 1)2 (3k − 1) k(k − 1)2 2(k−1)t (k − 1)(k2 − 3k + 1) (k−1)t 2k(k − 1)2 −t E B2 (t ) = e − e − e (2k − 1)2 (2k − 1)(3k − 1) (2k − 1)2 k2 k3 (k − 1) + e−kt − e−2kt . 2k − 1 (2k − 1)2 (3k − 1) 2
k2
We translate W (t ) to X (k) (t ) by their linear relationship (2), and we get the stated results.
5.2. Higher mixed moments We call E W i (t )Bj (t ) := mi,j (t ) a mixed moment of order i + j. In Section 5.1, we computed the mixed moments of orders 1 and 2. The exact calculation of mixed moments of higher order is indeed complicated, thus we solve the recurrence (4) only asymptotically. We write the differential equation (4) for the mixed moments in the following form
j−1 j
m′i,j (t ) =
s
s =0
(k − 1)j−s mi+1,s (t ) +
j−1 i i j r =0 s =0
r
s
(−1)j−s ki−r mr ,s+1 (t ) +
r =0
= j(k − 1)mi+1,j−1 (t ) − jmi,j (t ) + ikmi−1,j+1 (t ) +
j −2 j s=0
+
0≤r ≤i, 0≤s≤j−1 (r ,s)̸=(i,j−1)
i
j
r
s
(−1)j−s ki−r mr ,s+1 (t ) +
i −1 i
s
i −2 i r =0
r
r
ki−r mr ,j+1 (t ).
(k − 1)j−s mi+1,s (t )
ki−r mr ,j+1 (t ).
It is now clear that the differential equation for a particular mixed moment can be written as a linear combination of mixed moments of the same order plus a function of all lower order moments. So, we can solve simultaneously for all the mixed moments of the same order. The development will involve multifactorials, which are extensions (q-analogs) of the usual factorials and double factorials. For example, the 3-analog of the factorial of n is n!!! = n(n − 3) . . . (a + 3)a, where a = 3, if n is a multiple of 3; otherwise, a = n(mod 3). It is evident that if q is large, we will need long strings of exclamation marks. Therefore, we borrow a notation from formal languages; a string of q exclamation marks will be denoted by !(q) . In other words, formally speaking, the multifactorial of n with order q is defined recursively as n!
(q)
=
n, n((n − q)!(q) ),
if 0 < n ≤ q; , if n > q,
for any positive integers n and q. Subsequently, we have
E W i (t )Bj (t ) =
ki (k − 1)j ((k − 1)(i + j − 1) + 1)!(k−1) (i+j)(k−1)t e + O e(i+j−1)(k−1)t . (2k − 1)i+j
(5)
6
P. Zhang, H.M. Mahmoud / Statistics and Probability Letters 115 (2016) 1–7
Theorem 5.1. Let X (k) (t ) be the number of terminal nodes in an Apollonian network of index k at time t. As t → ∞, we have
(2k − 1)X (k) (t ) e(k−1)t
D
−→ Gamma
1 k−1
,k − 1 .
Proof. Setting j = 0 in (5), we have
E W i (t ) ∼
ki ((k − 1)(i − 1) + 1)!(k−1) i(k−1)t e , (2k − 1)i
which is equivalent to
E
(2k − 1)W (t )
i
ke(k−1)t
→ ((k − 1)(i − 1) + 1)!(k−1) .
Notice that
((k − 1)(i − 1) + 1)!(k−1) = ((k − 1)(i − 1) + 1)((k − 1)(i − 2) + 1) × · · · × 1 1 1 1 i−2+ ··· × (k − 1)i = i−1+ k−1 k−1 k−1 Γ i + k−1 1 (k − 1)i , = Γ k−1 1 which is the ith moment of a Gamma random variable with shape parameter k−1 1 and scale parameter k − 1. Translating W (t ) to X (k) (t ) via (2), we complete the proof. 6. Concluding remarks In the discrete-time Apollonian processes that arise in underlying Apollonian networks, the number of terminal nodes √ (k) Xn asymptotically has a normal distribution upon shifting to the center and scaling by n. By contrast, in the continuoustime Apollonian processes that arise in the Poissonized Apollonian networks, we only need to appropriately scale the corresponding random variable X (k) (t ) to obtain a limit, and that limit is a Gamma distributed random variable. We see that embedding in real time produces markedly different results. By methods identical to those we used to find the asymptotic Gamma distribution of white cliques in Poissonized Apollonian processes, we can find a Gamma distribution for the number of blue cliques,
(2k − 1)B(k) (t ) D 1 −→ Gamma ,k − 1 . (k − 1)e(k−1)t k−1 From Eq. (5), we see that
E
(2k − 1)W (t )
ke(k−1)t
(2k − 1)B(k) (t ) (k − 1)e(k−1)t
→ k.
Since the variance of each of the two limiting random variables is k − 1, we deduce that the correlation coefficient between them is 1. Therefore, the two limiting random variables are linearly related. Now we see that, even though each of the two limiting random variables has a Gamma distribution, the joint distribution does not fit any of the known types of bivariate Gamma distributions—as Balakrishnan and Lai (2009) lists 19 bivariate Gamma distributions, each of which has correlation strictly less than 1. We point out that the asymptotic results in this paper can be obtained by combining results from two sources: results from Janson (2004) and classic results of branching processes in Athreya and Ney (1972). Theorem 3.1 in Janson (2004) states existentially that, for a certain class of Pólya white–blue urn schemes, e−λ1 t
W (t ) a.s. −→ X B(t )
v1 , 1 − v1
for some random variable X , where λ1 is the principal (larger) eigenvalue, and
v1
1 − v1
is the corresponding eigenvector of
the transposed replacement matrix. The class for which this result holds is governed by conditions that guarantee tenability and irreducibility. We call the urn associated with Apollonian process the ‘‘Apollonian urn’’. This Apollonian urn falls in this class; after a calculation of the eigenvalues and eigenvectors, it can be ascertained that e−(k−1)t
X W (t ) a.s. −→ B(t ) 2k − 1
for some random variable X .
k
k−1
,
P. Zhang, H.M. Mahmoud / Statistics and Probability Letters 115 (2016) 1–7
7
On the other hand, the total number of balls in the Apollonian urn,
τ (t ) = W (t ) + B(t ), is shown in Janson (2004) to follow the relation
τ (t )
a.s.
e(k−1) t
−→ Gamma
1 k−1
,k − 1 .
(An existential proof is discussed in Athreya and Ney, 1972.) The variable X is determined to be this Gamma random variable. It is curious to note that the ratio of white to blue balls converges: W (t ) B(t )
a.s.
−→
k k−1
,
just like the discrete case, where W (n) B(n)
a.s.
−→
k k−1
.
Limit ratios are the same in both the discrete Apollonian urn scheme and the Apollonian process, even though the distributions are quite different—they are normal in the discrete case, and Gamma in the process. The present manuscript explores the potential on an alternative differential method. The proposed method is capable of producing exact moments (hence rates of convergence to lower order terms). References Andrade, J.S.J., Herrmann, H.J., Andrade, R.F., da Silva, L.R., 2005. Apollonian networks: Simultaneously scale-free, small world, euclidean, space filling, and with matching graphs. Phys. Rev. Lett. 94, 018702. http://dx.doi.org/10.1103/PhysRevLett.94.018702. URL: http://journals.aps.org/prl/abstract/10.1103/PhysRevLett.94.018702. Athreya, K., Karlin, S., 1968. Embedding of urn schemes into continuous time markov branching processes and related limit theorems. Ann. Math. Stat. 39, 1801–1817. http://dx.doi.org/10.1214/aoms/1177698013. Athreya, K., Ney, P., 1972. Branching Process. Springer-Verlag. Bagchi, A., Pal, A.K., 1985. Asymptotic normality in the generalized pólya–eggenberger urn model, with an application to computer data structures. SIAM J. Algebr. Discrete Methods 6, 394–405. http://dx.doi.org/10.1137/0606041. URL: http://epubs.siam.org/doi/abs/10.1137/0606041. Balaji, S., Mahmoud, H., 2006. Exact and limiting distributions in diagonal Pólya processes. Ann. Inst. Statist. Math. 58, 171–185. http://dx.doi.org/10.1007/s10463-006-0093-1. URL: http://www.springerlink.com/content/4347u64681j64568/abstract/. Balaji, S., Mahmoud, H., Watanabe, O., 2006. Distributions in the ehrenfest process. Statist. Probab. Lett. 76, 666–674. http://dx.doi.org/10.1016/j.spl.2005.09.013. URL: http://www.sciencedirect.com/science/article/pii/S0167715205003676. Balakrishnan, N., Lai, C., 2009. Continuous Bivariate Distributions, second ed. Springer. Edwards, C., Penney, D., Differential Equations and Linear Algebra, third ed., 2008. Janson, S., 2004. Functional limit theorems for multitype branching processes and generalized pólya urns. Stochastic Process. Appl. 110, 177–245. http://dx.doi.org/10.1016/j.spa.2003.12.002. URL: http://www.sciencedirect.com/science/article/pii/S0304414903001790. Janson, S., 2005. Asymptotic degree distribution in random recursive trees. Random Struct. Algorithms 26, 69–83. http://dx.doi.org/10.1002/rsa.20046. URL: http://onlinelibrary.wiley.com/doi/10.1002/rsa.20046/abstract. Johnson, N., Kotz, S., 1977. Urn Models and Their Application: An Approach to Modern Discrete Probability Theory. Wiley. Mahmoud, H., 2008. Pólya Urn Models. CRC Press. Sparks, J., Mahmoud, H., 2013. Pahses in the two-color tenable zero-balanced Pólya processes. Statist. Probab. Lett. 83, 265–271. http://dx.doi.org/10.1016/j.spl.2012.08.020. URL: http://www.sciencedirect.com/science/article/pii/S0167715212003252. Zhang, P., Mahmoud, H., 2016. The degree profile and weight in Apollonian networks and k-trees. Adv. Appl. Probab. 48, 163–175.
Contents lists available at ScienceDirect
Statistics and Probability Letters journal homepage: www.elsevier.com/locate/stapro
Distributions in a class of Poissonized urns with an application to Apollonian networks Panpan Zhang ∗ , Hosam M. Mahmoud Department of Statistics, The George Washington University, Washington, DC 20052, USA
article
info
Article history: Received 17 October 2015 Received in revised form 23 March 2016 Accepted 25 March 2016 Available online 8 April 2016 MSC: 60G20 65Q30 35A99 Keywords: Pólya urn Pólya process Partial differential equation Stochastic recurrence
abstract We study a class of Pólya processes that underlie terminal nodes in a random Apollonian network. We calculate the exact first and second moments of the number of terminal nodes by solving ordinary differential equations. These equations are derived from the partial differential equation governing the process. In fact, the partial differential equation yields a stochastic hierarchy of moment equations, which can be bootstrapped to get higher moments from the equations that have been solved for lower moments. We also show that the number of terminal nodes, when appropriately scaled, converges in distribution to a gamma random variable via the method of moments. The asymptotic results can be obtained using classic methods of branching processes. The manuscript explores the potential of an alternative method capable of producing exact moments and rates of convergence. © 2016 Elsevier B.V. All rights reserved.
1. Introduction Urn schemes growing in discrete time have applications in many fields. Early applications include modeling epidemiological contagion (the Pólya-Eggenberger scheme) and modeling gas diffusion and heat transfer (the Ehrenfest model). Throughout the last century, many more applications were proposed, literally hundreds, and it is not possible to concisely list them all, but we direct the interested reader to the numerous applications discussed in Johnson and Kotz (1977) and Mahmoud (2008). Associated with a Pólya urn scheme is a process obtained by embedding in continuous time. Processes growing in continuous time are more realistic and applicable than those growing in discrete time. In practice, there is no guarantee that the epochs of process will be equispaced. They also address a different question—that of the status after a certain amount of (real) time, rather than the status after a certain number of events. For better readability, we call the evolution in continuous time the Pólya process (as in Balaji and Mahmoud, 2006; Sparks and Mahmoud, 2013), and reserve the term ‘‘scheme’’ for urns evolving in discrete time. The Pólya process was introduced in Athreya and Karlin (1968) (though it was not named that) as a mathematical transform (Poissonization) to shed light on urns evolving in discrete time. Athreya and Karlin (1968) reports that dePoissonization, to obtain the inverse transform, is highly nontrivial, and until today a rigorous inverse transform remains elusive.
∗
Corresponding author. Tel.: +1 3472858737. E-mail addresses: [email protected] (P. Zhang), [email protected] (H.M. Mahmoud).
http://dx.doi.org/10.1016/j.spl.2016.03.023 0167-7152/© 2016 Elsevier B.V. All rights reserved.
2
P. Zhang, H.M. Mahmoud / Statistics and Probability Letters 115 (2016) 1–7
2. The Pólya process A two-color Pólya process is an urn containing balls of up to two colors, say white and blue, and evolving over time according to certain rules. Balls are drawn at epochs in time with independent Exp(1) interarrival times.1 When a white ball is drawn, it is placed back in the urn together with a white balls and b blue balls; when a blue ball is drawn, it is placed back in the urn together with c white balls and d blue balls. In general these numbers can be positive, zero or negative (interpret adding a negative number of balls to mean taking that many balls out of the urn). An urn process is called tenable, if it is always possible to execute the rules; that is, there never comes a time when the urn becomes empty or the rules call for removing a number of balls of a color more than their number in the urn. It is customary to represent these dynamics by the replacement matrix
a c
b . d
The rows of this matrix are indexed with white and blue, respectively from top to bottom; the columns are indexed with white and blue, respectively from left to right. To be precise about the dynamics of evolution in continuous time, we add a few words of explanation. Each ball in the process is endowed with a clock that rings in exponential time Exp(1). All the clocks are independent of each other and of any other random variables related to the past. At a renewal point (an epoch), we execute the rules associated with its color instantaneously. All new balls added come endowed with their own independent clocks. By the memorylessness of the exponential interarrival time, the Pólya process is Markovian. However, in general, it is an inhomogeneous jump process, as the rates keep changing owing to ball additions. We are interested in the class of two-color Pólya processes with the replacement matrix
0 k
k−1 , −1
(1)
for integer k ≥ 3. The study is motivated by an application from the growth in real time of Apollonian networks. We call this class of Pólya processes the Apollonian process. 3. Organization The rest of this paper is organized as follows. In Section 4, we introduce a discrete-time Apollonian network. A coloring scheme of its cliques connects it to Pólya urns. The Poissonized counterpart is discussed in Section 5. In Section 5.1, we take up the mean and variance of the number of terminal nodes at time t by finding the joint moment of white and blue cliques defined in Section 4. Higher asymptotic moments are derived in Section 5.2, leading to a conclusion that the scaled number of terminal nodes at time t asymptotically has a Gamma distribution. Section 6 concludes the paper with remarks comparing the behavior of terminal nodes in both discrete- and continuous-time Apollonian networks. We also mention connections to branching processes that can alternatively be utilized to get the asymptotics. 4. Apollonian networks Apollonian networks are graphs based on triangulation. More generally, Apollonian networks are constructed on subdivision of simplexes of higher order. As random graphs, Apollonian networks have the following features: scale free, small world, Euclidean, space filling and with matching graphs. Each of these features has a significant impact in current random network research. We refer the reader to Andrade et al. (2005) for motivation and important applications of Apollonian networks. For purposes of initial introduction, we shall use triangulation, but our intent is to study general simplexes. Consider an active cycle of length 3 (with three edges, and three vertices labeled with three zeros to indicate they are the initial nodes). The term ‘‘active’’ indicates that it is a ‘‘recruiter’’ of other self-similar graphs. A random Apollonian network (of index 3) is a random graph that grows as follows. At the first step, we insert a new vertex (labeled with 1) at an arbitrary position inside the original cycle, connect it with three edges to the three initial nodes (labeled with 0), and we then discard the recruiter. We thus have three new active cycles after the first step of insertion, and the initial cycle becomes inactive and it does not recruit any more. In the second and subsequent steps, we choose an active cycle at random (all active cycles being equally likely), and insert a new vertex inside the cycle. We then connect the newly inserted vertex with three edges to the three nodes of the chosen cycle, creating three new active cycles. Lastly, we deactivate the recruiting cycle. At step n ≥ 1, the newly added vertex is labeled with n. For a general index k ≥ 3, the formal definition is the following. A random Apollonian network of index k is a random graph that grows out of a k-clique with vertices labeled with zeros. At the nth step, a k-clique is chosen at random. We insert
1 Exp(1) refers to an exponential random variable with parameter 1.
P. Zhang, H.M. Mahmoud / Statistics and Probability Letters 115 (2016) 1–7
3
a new vertex in the graph, connect it with k edges to the k vertices of the chosen clique, and we then discard the recruiting clique. For more motivational information about Apollonian networks, we refer the reader to Zhang and Mahmoud (2016). In various evolutionary random graphs, the number of terminal nodes—nodes that appeared but did not participate yet in further recruiting—is a significant property. For instance, leaves are the terminal nodes in random trees, and they were extensively studied in many tree families; see for example Janson (2005). In an Apollonian network, we consider a node terminal, if all the cliques that are incident with this node have not yet recruited any other cliques. (k) Let Xn be the number of terminal nodes in an Apollonian network with index k after n node insertions. We devise a (k) two-color Pólya urn strategy to determine the distribution of Xn . Let us color cliques incident with terminal nodes with white (W ), and all other cliques with blue (B). Consider the network at insertion step 0 to be the initial condition of an urn, having no white cliques, and one blue. Note that no two terminal nodes ever appear in the same clique. So, a count of white cliques (Wn ) at insertion step n translates directly into a count of terminal nodes, namely Xn(k) =
Wn
. (2) k When a white clique recruits, the terminal node in it is no longer terminal, and the recruiter turns inactive. The other k − 1 white cliques incident with this chosen terminal node turn into blue. However, all k newly born cliques are white; on net there is no change in the number of white cliques. On the other hand, if a blue clique recruits, we lose one blue clique as the recruiter becomes inactive; meanwhile, k new active white cliques appear. Therefore, we have a corresponding Pólya urn scheme with the replacement matrix as given in (1). This is an instance of the well-studied Bagchi–Pal urn (see Bagchi and Pal, 1985), and we can conduct a probabilistic analysis based on a known urn theory (see Athreya and Karlin, 1968; Bagchi and Pal, 1985). (k)
Proposition 4.1 (Zhang and Mahmoud, 2016). Let Xn age n. We then have (k)
Xn
n
a.s.
−→
k−1 2k − 1
be the number terminal nodes in an Apollonian network of index k at
,
and (k) k−1 Xn − 2k n −1
√
n
D
−→ N
0,
k(k − 1)2
.
(2k − 1)2 (3k − 1)
5. Poissonized Apollonian networks We next consider the analogue of Apollonian networks obtained by embedding processes in continuous time, which may be a more realistic reflection of applications. In other words, we are Poissonizing the urn with replacement matrix (1). In an Apollonian process generated from an Apollonian network of index k, let X (k) (t ) be the number of terminal nodes at time t, and W (t ) and B(t ) be the number of respectively white and blue cliques based on the color code we introduced in the previous section. The Apollonian process is the two-component vector
φ(t , u, v) := E e
W (t ) B(t )
. Let
uW (t )+v B(t )
be the joint moment generating function of
W (t ) B(t )
. For a general tenable two-color Pólya process, a certain partial
differential equation (PDE) is known; see Balaji and Mahmoud (2006). It is derived for a Pólya process with a deterministic replacement matrix. The equation is generalized in Mahmoud (2008) to Pólya processes with random replacement matrices. For deterministic replacement matrices, the PDE is
∂φ(t , u, v) ∂φ(t , u, v) ∂φ(t , u, v) + 1 − eau+bv + 1 − ecu+dv = 0. ∂t ∂u ∂v This functional equation has only been solved in a very limited number of cases, as for instance the case of forward and backward diagonal processes (see Balaji and Mahmoud, 2006), and the Ehrenfest process (see Balaji et al., 2006). The quest continues to find more classes of Pólya processes amenable to analysis via this method. As we shall see, the Apollonian process is one such case. Specialized to the number of terminal nodes of an Apollonian network, the PDE is
∂φ(t , u, v) ∂φ(t , u, v) ∂φ(t , u, v) + 1 − e(k−1)v + 1 − eku−v = 0. ∂t ∂u ∂v
(3)
Toward ordinary differential equations for moments, for any nonnegative integers i and j, we take the partial derivative i times for u and j times for v , respectively, for each term in (3). Setting u = v = 0, we obtain a recurrence for the mixed joint
4
P. Zhang, H.M. Mahmoud / Statistics and Probability Letters 115 (2016) 1–7
moments E W i (t )Bj (t ) as follows:
d
j−1 j E W (t )B (t ) = (k − 1)j−s E W i+1 (t )Bs (t )
dt
i
j
s
s =0
j−1 i i−1 r i j i i−r r j−s i−r s +1 + (−1) k E W (t )B (t ) + k E W (t )Bj+1 (t ) .
r
r =0 s=0
s
r =0
(4)
r
5.1. The mean and variance For small values of i and j, We are able to determine the exact moments of X (k) (t ) by solving the ordinary differential recurrence equations for W (t ). We illustrate this paradigm with the mean and variance where we have 1 ≤ i + j ≤ 2. Proposition 5.1. Let X (k) (t ) be the number of terminal nodes in a Poissonized Apollonian network of index k at time t. We have
E X (k) (t ) =
1
Var X (k) (t ) =
2k − 1 k−1
(2k − 1) −
1
e(k−1)t −
2
2k − 1
e−kt , k−1
e2(k−1)t −
1
(2k − 1)(3k − 1) k2 + 2k − 1 − e−2kt . (2k − 1)2 (3k − 1)
e−kt
(2k − 1)
e(k−1)t +
2k
(2k − 1)2
e−t
Proof. We determine the expectation and variance of X (k) (t ) via the expectation and variance of W (t ). We compute the first moment of W (t ) by a theorem derived in Mahmoud (2008) as follows:
0 E W (t ) =e k E B(t )
⊤
k−1 −1
W (0) , B(0)
t
where ⊤ is the transpose of a matrix, and W (0) and B(0) are initial numbers of white and blue cliques, respectively. We thus have
0 E W (t ) k−1 =e E B(t )
k −1
t
0 . 1
Reading off the components of this vector one by one, we have k k e(k−1)t − 2k − 1 2k − 1 k − 1 (k−1)t k E B(t ) = e + 2k − 1 2k − 1
E W (t ) =
e−kt , e−kt .
We resort to the established recurrence (4) for the joint moments of W (t ) and B(t ) to compute the second moment of W (t ). We simplify the recurrence for the three cases {i = 2, j = 0}, {i = 1, j = 1}, and {i = 0, j = 2}, and obtain the following matrix form
2 0 E W (t ) E W (t )B(t ) = k − 1 dt 0 E B2 (t )
d
2k −1 2(k − 1)
2 k2 E B(t ) 0 E W (t ) k E W (t )B(t ) + −k E B(t ) . −2 E B2 (t ) (k − 1)2 E W (t ) + E B(t )
This is a nonhomogeneous differential linear equation for
E[W 2 (t )] E[W (t )B(t )] E[B2 (t )]
. Notice that the coefficient matrix
is invertible. Hence, for all k ≥ 3, we have a general unique solution to this type of linear equation:2 k t
2 E W (t ) E W (t )B(t ) = e 0 E B2 (t )
0
−1 0
2 See Edwards and Penney (2008) for details.
2k −1 2k − 1
0 k (t −s) −2
k2 E B(s) −k E B(s) ds 2 (k − 1) E W (s) + E B(s)
0 2k 0 k − 1 −1 k 0 2k − 1 −2
P. Zhang, H.M. Mahmoud / Statistics and Probability Letters 115 (2016) 1–7
0
k
−1
+e
2k
0 k t E W 2 (0) −2 E W (0)B(0) , E B2 (0)
−1 2k − 1
0
where the initial vector condition is
5
E[W 2 (0)] E[W (0)B(0)] E[B2 (0)]
0 =
0 1
. We plug in the results of E W (t ) and E B(t ) , compute the
matrix exponential and the integration, and get the following simplified solutions for E W 2 (t ) , E W (t )B(t ) and E B2 (t ) . We obtain
k3
E W 2 (t ) =
k2 (k − 1)
e2(k−1)t −
(2k − 1)
e(k−1)t +
2k2 (k − 1) −t e (2k − 1)2
(2k − 1)(3k − 1) k3 (k − 1) − e−kt − e−2kt , 2k − 1 (2k − 1)2 (3k − 1) k2 (k − 1) 2(k−1)t k(k − 1) (k−1)t k(k − 1) −t k3 (k − 1) − − + e e e e−2kt , E W (t )B(t ) = (2k − 1)2 3k − 1 (2k − 1)2 (2k − 1)2 (3k − 1) k(k − 1)2 2(k−1)t (k − 1)(k2 − 3k + 1) (k−1)t 2k(k − 1)2 −t E B2 (t ) = e − e − e (2k − 1)2 (2k − 1)(3k − 1) (2k − 1)2 k2 k3 (k − 1) + e−kt − e−2kt . 2k − 1 (2k − 1)2 (3k − 1) 2
k2
We translate W (t ) to X (k) (t ) by their linear relationship (2), and we get the stated results.
5.2. Higher mixed moments We call E W i (t )Bj (t ) := mi,j (t ) a mixed moment of order i + j. In Section 5.1, we computed the mixed moments of orders 1 and 2. The exact calculation of mixed moments of higher order is indeed complicated, thus we solve the recurrence (4) only asymptotically. We write the differential equation (4) for the mixed moments in the following form
j−1 j
m′i,j (t ) =
s
s =0
(k − 1)j−s mi+1,s (t ) +
j−1 i i j r =0 s =0
r
s
(−1)j−s ki−r mr ,s+1 (t ) +
r =0
= j(k − 1)mi+1,j−1 (t ) − jmi,j (t ) + ikmi−1,j+1 (t ) +
j −2 j s=0
+
0≤r ≤i, 0≤s≤j−1 (r ,s)̸=(i,j−1)
i
j
r
s
(−1)j−s ki−r mr ,s+1 (t ) +
i −1 i
s
i −2 i r =0
r
r
ki−r mr ,j+1 (t ).
(k − 1)j−s mi+1,s (t )
ki−r mr ,j+1 (t ).
It is now clear that the differential equation for a particular mixed moment can be written as a linear combination of mixed moments of the same order plus a function of all lower order moments. So, we can solve simultaneously for all the mixed moments of the same order. The development will involve multifactorials, which are extensions (q-analogs) of the usual factorials and double factorials. For example, the 3-analog of the factorial of n is n!!! = n(n − 3) . . . (a + 3)a, where a = 3, if n is a multiple of 3; otherwise, a = n(mod 3). It is evident that if q is large, we will need long strings of exclamation marks. Therefore, we borrow a notation from formal languages; a string of q exclamation marks will be denoted by !(q) . In other words, formally speaking, the multifactorial of n with order q is defined recursively as n!
(q)
=
n, n((n − q)!(q) ),
if 0 < n ≤ q; , if n > q,
for any positive integers n and q. Subsequently, we have
E W i (t )Bj (t ) =
ki (k − 1)j ((k − 1)(i + j − 1) + 1)!(k−1) (i+j)(k−1)t e + O e(i+j−1)(k−1)t . (2k − 1)i+j
(5)
6
P. Zhang, H.M. Mahmoud / Statistics and Probability Letters 115 (2016) 1–7
Theorem 5.1. Let X (k) (t ) be the number of terminal nodes in an Apollonian network of index k at time t. As t → ∞, we have
(2k − 1)X (k) (t ) e(k−1)t
D
−→ Gamma
1 k−1
,k − 1 .
Proof. Setting j = 0 in (5), we have
E W i (t ) ∼
ki ((k − 1)(i − 1) + 1)!(k−1) i(k−1)t e , (2k − 1)i
which is equivalent to
E
(2k − 1)W (t )
i
ke(k−1)t
→ ((k − 1)(i − 1) + 1)!(k−1) .
Notice that
((k − 1)(i − 1) + 1)!(k−1) = ((k − 1)(i − 1) + 1)((k − 1)(i − 2) + 1) × · · · × 1 1 1 1 i−2+ ··· × (k − 1)i = i−1+ k−1 k−1 k−1 Γ i + k−1 1 (k − 1)i , = Γ k−1 1 which is the ith moment of a Gamma random variable with shape parameter k−1 1 and scale parameter k − 1. Translating W (t ) to X (k) (t ) via (2), we complete the proof. 6. Concluding remarks In the discrete-time Apollonian processes that arise in underlying Apollonian networks, the number of terminal nodes √ (k) Xn asymptotically has a normal distribution upon shifting to the center and scaling by n. By contrast, in the continuoustime Apollonian processes that arise in the Poissonized Apollonian networks, we only need to appropriately scale the corresponding random variable X (k) (t ) to obtain a limit, and that limit is a Gamma distributed random variable. We see that embedding in real time produces markedly different results. By methods identical to those we used to find the asymptotic Gamma distribution of white cliques in Poissonized Apollonian processes, we can find a Gamma distribution for the number of blue cliques,
(2k − 1)B(k) (t ) D 1 −→ Gamma ,k − 1 . (k − 1)e(k−1)t k−1 From Eq. (5), we see that
E
(2k − 1)W (t )
ke(k−1)t
(2k − 1)B(k) (t ) (k − 1)e(k−1)t
→ k.
Since the variance of each of the two limiting random variables is k − 1, we deduce that the correlation coefficient between them is 1. Therefore, the two limiting random variables are linearly related. Now we see that, even though each of the two limiting random variables has a Gamma distribution, the joint distribution does not fit any of the known types of bivariate Gamma distributions—as Balakrishnan and Lai (2009) lists 19 bivariate Gamma distributions, each of which has correlation strictly less than 1. We point out that the asymptotic results in this paper can be obtained by combining results from two sources: results from Janson (2004) and classic results of branching processes in Athreya and Ney (1972). Theorem 3.1 in Janson (2004) states existentially that, for a certain class of Pólya white–blue urn schemes, e−λ1 t
W (t ) a.s. −→ X B(t )
v1 , 1 − v1
for some random variable X , where λ1 is the principal (larger) eigenvalue, and
v1
1 − v1
is the corresponding eigenvector of
the transposed replacement matrix. The class for which this result holds is governed by conditions that guarantee tenability and irreducibility. We call the urn associated with Apollonian process the ‘‘Apollonian urn’’. This Apollonian urn falls in this class; after a calculation of the eigenvalues and eigenvectors, it can be ascertained that e−(k−1)t
X W (t ) a.s. −→ B(t ) 2k − 1
for some random variable X .
k
k−1
,
P. Zhang, H.M. Mahmoud / Statistics and Probability Letters 115 (2016) 1–7
7
On the other hand, the total number of balls in the Apollonian urn,
τ (t ) = W (t ) + B(t ), is shown in Janson (2004) to follow the relation
τ (t )
a.s.
e(k−1) t
−→ Gamma
1 k−1
,k − 1 .
(An existential proof is discussed in Athreya and Ney, 1972.) The variable X is determined to be this Gamma random variable. It is curious to note that the ratio of white to blue balls converges: W (t ) B(t )
a.s.
−→
k k−1
,
just like the discrete case, where W (n) B(n)
a.s.
−→
k k−1
.
Limit ratios are the same in both the discrete Apollonian urn scheme and the Apollonian process, even though the distributions are quite different—they are normal in the discrete case, and Gamma in the process. The present manuscript explores the potential on an alternative differential method. The proposed method is capable of producing exact moments (hence rates of convergence to lower order terms). References Andrade, J.S.J., Herrmann, H.J., Andrade, R.F., da Silva, L.R., 2005. Apollonian networks: Simultaneously scale-free, small world, euclidean, space filling, and with matching graphs. Phys. Rev. Lett. 94, 018702. http://dx.doi.org/10.1103/PhysRevLett.94.018702. URL: http://journals.aps.org/prl/abstract/10.1103/PhysRevLett.94.018702. Athreya, K., Karlin, S., 1968. Embedding of urn schemes into continuous time markov branching processes and related limit theorems. Ann. Math. Stat. 39, 1801–1817. http://dx.doi.org/10.1214/aoms/1177698013. Athreya, K., Ney, P., 1972. Branching Process. Springer-Verlag. Bagchi, A., Pal, A.K., 1985. Asymptotic normality in the generalized pólya–eggenberger urn model, with an application to computer data structures. SIAM J. Algebr. Discrete Methods 6, 394–405. http://dx.doi.org/10.1137/0606041. URL: http://epubs.siam.org/doi/abs/10.1137/0606041. Balaji, S., Mahmoud, H., 2006. Exact and limiting distributions in diagonal Pólya processes. Ann. Inst. Statist. Math. 58, 171–185. http://dx.doi.org/10.1007/s10463-006-0093-1. URL: http://www.springerlink.com/content/4347u64681j64568/abstract/. Balaji, S., Mahmoud, H., Watanabe, O., 2006. Distributions in the ehrenfest process. Statist. Probab. Lett. 76, 666–674. http://dx.doi.org/10.1016/j.spl.2005.09.013. URL: http://www.sciencedirect.com/science/article/pii/S0167715205003676. Balakrishnan, N., Lai, C., 2009. Continuous Bivariate Distributions, second ed. Springer. Edwards, C., Penney, D., Differential Equations and Linear Algebra, third ed., 2008. Janson, S., 2004. Functional limit theorems for multitype branching processes and generalized pólya urns. Stochastic Process. Appl. 110, 177–245. http://dx.doi.org/10.1016/j.spa.2003.12.002. URL: http://www.sciencedirect.com/science/article/pii/S0304414903001790. Janson, S., 2005. Asymptotic degree distribution in random recursive trees. Random Struct. Algorithms 26, 69–83. http://dx.doi.org/10.1002/rsa.20046. URL: http://onlinelibrary.wiley.com/doi/10.1002/rsa.20046/abstract. Johnson, N., Kotz, S., 1977. Urn Models and Their Application: An Approach to Modern Discrete Probability Theory. Wiley. Mahmoud, H., 2008. Pólya Urn Models. CRC Press. Sparks, J., Mahmoud, H., 2013. Pahses in the two-color tenable zero-balanced Pólya processes. Statist. Probab. Lett. 83, 265–271. http://dx.doi.org/10.1016/j.spl.2012.08.020. URL: http://www.sciencedirect.com/science/article/pii/S0167715212003252. Zhang, P., Mahmoud, H., 2016. The degree profile and weight in Apollonian networks and k-trees. Adv. Appl. Probab. 48, 163–175.