Birth-Death Processes, Lattice Path Combinatorics,
Continued Fractions, and Orthogonal Polynomials

Fabrice Guillemin

France Telecom, CNET, Lannuon, Breizh

Algorithms Seminar

February 2, 1998

[summary by Philippe Flajolet and Fabrice Guillemin]

A properly typeset version of this document is available in postscript and in pdf.

If some fonts do not look right on your screen, this might be fixed by configuring your browser (see the documentation here).

Abstract
Classic works of Karlin-McGregor and Jones-Magnus have established a fully general correspondence between birth-death processes and continued fractions of the Stieltjes-Jacobi type together with their associated orthogonal polynomials. This fundamental correspondence can be revisited in the light of the otherwise known combinatorial correspondence between weighted lattice paths and continued fractions. For birth-death processes, this approach separates clearly the formal apparatus from the analytic-probabilistic machinery and neatly delineates those parameters that are amenable to a treatment by means of continued fractions and orthogonal polynomials.



1   Birth-Death Processes

Consider a particle initially in state 0 that, at any given time, may change to another state 1 (where it stays), with rate l. This means that the probability of a state change in an interval of time of length dt is l dt. Then, the probability p0(t) that the particle is still in state 0 at time t satisfies
p0(t+dt)-p0(t)=-l p0(t) dt
or p0'(t)=-l p0(t), whose solution is an exponential distribution,
p0(t)=e
-l t
 
.

Similarly, a particle initially in state 0 that may change either to state 1 with rate l or to state -1 with rate µ will satisfy (pj(t) is the probability of being in state j at time t)
p0(t)=e
-(l+µ)t
 
,     p1(t)=
l
l+µ
(1-e
-(l+µ)t
 
),      p-1(t)=
µ
l+µ
(1-e
-(l+µ)t
 
).
The interpretation is obvious: the particle stays in state 0 for a random amount of time with an exponential distribution of rate l+µ and then changes to states -1,+1 with probabilities equal to to l/(l+µ) and µ(l+µ).

In a general birth-death process a particle can be in any state in {0,1,2,...} and when in state j, it can only change to state j+1 at rate lj or to state j-1 at rate µj. By analogy with the model of an evolving population (whose size is represented by the state), the lj are called birth rates and the µj death rates. The general problem is to understand the evolution of a process given values (or properties) of its birth and death rates; see [12, Ch. 4] for an excellent introduction.

Let pn(t) be the probability of being in state n at time t. An essential rôle is played by the coefficients
pn =
l0l1··· ln-1
µ1µ2··· µn
.
Indeed, a classical result asserts that the process is ergodic (the expected time to return from each state to itself is finite) if and only if
 
å
n³1
pn<¥,    
 
å
n³0
1
lnpn
=+¥.
(The first condition ensures the existence of an invariant measure for the embedded discrete-time Markov chain; the second one guarantees that, in the continuous-time process, the particle is not absorbed at infinity in finite time.) In that case, one has
pn:=
 
lim
t®¥
pn(t)=
pn
 
å
n³1
pn
,
where these quantities represent the long run probability of being in state n.

More puzzling is the nonstationary behaviour of the process that is described by the infinite-dimensional differential system
p'j(t)=l j-1pj-1(t)-(lj+µj)pj(t)+µj+1pj+1(t),      pj(0)=dj,0.     (1)
Although finite-dimensional versions are ``easy'' and reduce to combinations of exponentials, it is precisely the infinite-dimensional character of the system that renders its analysis interesting.

In a series of important papers, Karlin and McGregor [10, 11] have developed a general connection between the fundamental system (1) and an associated family of orthogonal polynomials. Later, Jones and Magnus constructed a direct continued fraction representation; see [8, 9].

This summary is an account of Guillemin's lecture (see [5, 6]), as well as of later developments. The point of view that is adopted here consists in relating the combinatorial theory of lattice paths to birth-death processes in the following way: (i) trajectories of birth-death processes are precisely lattice paths; (ii) lattice paths have generating functions expressed as continued fractions; (iii) the Laplace transform expresses the main parameters of birth-death processes as weighted lattice paths to which the combinatorial theory applies.

2   Lattice Paths and Continued Fractions

It is known that the formal theory of continued fraction expansions for power series is identical to the combinatorial theory of weighted lattice paths; see [1, 2, 4]. Define a path u=(U0,U1,...,Un) to be a sequence of points in the lattice N×N such that if Uj=(xj,yj), then xj=j and |yj+1-yj|=1. If successive points are connected by edges, then an edge can only be an ascent (a: yj+1-yj=+1), a descent (b: yj+1-yj=-1), or a level step (c: yj+1-yj=0). Thus a path is always nonnegative and by a horizontal translation, one may always assume that x0=0. A path can be encoded by a word with a,b,c representing the three types of steps. What we call the standard encoding is such a word in which each step a,b,c is subscripted by the value of the y-coordinate of its associated point. For instance,
w=a0a1a2b3c2c2a2b3b2b1a0c1
encodes a path that connects the source U0=(0,0) to the destination U12=(12,1). We freely identify a path u defined as a sequence of points, its word encoding w, and the corresponding monomial.

We consider various geometric conditions that may be imposed on paths: Hk,l is the collection of all paths that connect a source at altitude k to a destination at altitude l, H[£ h] denotes paths of height (maximal altitude) at most h, etc.
Theorem 1   The collection H0,0 of all paths has generating function
H0,0=
1
1-c0-
a0b1
1-c1-
a1b2
1-c2-
a2b3
·  
 · 
  ·

Proof. It suffices to observe that (1-f)-1=1+f+f2+··· generates symbolically all the sequences with components f. For instance, in H0,0, the expressions
1
1-c0
,     
1
1-c0-a0b1
,    
1
1-c0-
a0b1
1-c1
    (2)
generate successively paths composed from c0 level steps only, paths of height at most 1 without c1 steps, all paths of height at most 1. The complete continued fraction representation is easily built by stages in a similar fashion.


In particular, the collection of all paths from level 0 to level 0 with height at most h is
H0,0[< h]=
Ph
Qh
,     (3)
a rational fraction, whose numerators and denominators, Ph,Qh, each satisfy the recurrence
yh+1=(1-ch)yh-ah-1bhyh-1,
with Q-1=P0=0, Q0=P1=1. (Linear fractional transformations are 2×2 matrices in disguise!)

Well-known path decompositions, like those based on first or last time at which levels are reached, can then be used provided they are combinatorially ``unambiguous''. This and simple manipulations on linear fractional transformations give access to many geometric constraints in addition to (2) and (3). We cite here some representative identities from [1, 2],
H0,h-1[<h]=
a0a1··· ah-1
Qh
,      H0,k=
1
b1b2··· bk
( QkH0,0-Pk ) ,     (4)
Hk,l=
Qk
a0··· ak-1b1··· bl
( QlH0,0-Pl ) ,     (5)
where the latter holds provided k£ l.

The forms (2), (3) (4), (5) can be converted into bona fide counting generating functions of paths weighted multiplicatively by means of the combinatorial morphism,
c(ak)=akz,     c(bk)=bk z,     c(ck)=gkz.
In that case, the continued fraction (2) becomes the general fraction of the J-type (for Jacobi); see  [7, 9, 13].

3   The Connection

We illustrate here in its simplest form the many-faceted connection between birth-death processes and continued fractions. It was apparently first stated explicitly by Jones and Magnus but it is implicit in earlier works of Karlin and McGregor. The connection goes through the probabilities pi,j(t) of being in state j at time t starting from state i and the Laplace transforms,
Pi,j(s) = ó
õ
¥


0
pi,j(t) e-st  dt.

Theorem 2   The Laplace transform of the probability of return to the origin satisfies
P0,0(s)=
1
l0+s-
l0µ1
l1+µ1+s-
l1µ2
·  
 · 
  ·
.
We offer here two proofs. A third proof that is based on ``uniformization of time'' can also be given but is omitted in this note.

Proof.[Proof 1]Take the Laplace transform of the fundamental system (1) (so that pj(t)=p0,j(t)) and use the induced relations on the ratios P0,r/P0,r+1. This proof is the most direct but the least illuminating from a structural standpoint. In particular, this proof does not provide an immediate grasp on the question of deciding which parameters are amenable to continued fraction representations.


Proof.[Proof 2] Examine the times at which the (continuous time) birth-death process {Lt} changes states. This defines an embedded (discrete time) Markov chain {Yn}. Then the set of trajectories of the chain {Yn} is exactly the family of lattice paths of Section 2. The method consists in splitting the probabilities by conditioning according to all legal trajectories. The different sojourn times are independent by the nature of the process (the strong Markov property satisfied by {Lt}). Also, sums of independent random variables correspond to products of Laplace transforms. Thus, the Laplace transform of the probability in the continuous model of following a path w has a product form; for instance, to w=a0a1b2a1, there corresponds the transform
æ
ç
ç
è
l0
l0+µ0
l1
l1+µ1
µ2
l2+µ2
l1
l1+µ1
ö
÷
÷
ø		 · æ
ç
ç
è
l0+µ0
s+l0+µ0
l1+µ1
s+l1+µ1
l2+µ2
s+l2+µ2
l1+µ1
s+l1+µ1
ö
÷
÷
ø		 .
Thus, the Laplace transform P0,0(s) is, apart from a fudge factor of 1/(s+l0), a sum over all paths lattice from zero to zero weighted multiplicatively by the probabilistic morphism,
c'(aj)=
lj
s+lj+µj
,      c'(bj)=
µj
s+lj+µj
,     (6)
with c'(cj)=0. In other words, one has P0,0(s)=c'(H0,0)1/s+l0, and the statement follows.


The same method applies to the computation of transition probabilities, the analysis of maximum height, and so on. For instance, the probability of reaching state k has
P0,k(s)=
1
µ1µ2··· µk
( Ak(s)P0,0(s)-Bk(s) ) ,
where Ak/Bk is the kth convergent of the continued fraction that represents P0,0, so that Ak,Bk are simple variants of c'(Pk),c'(Qk).

Orthogonality.
In the case of paths, the reciprocals of the Qh polynomials, Qh(z)=zhc(Q)(z-1) are formally orthogonal with respect to a measure defined its moments,
L [zn] º ó
õ zn dµ (z) = H0,0,n.     (7)
Formal aspects of paths and orthogonality are detailed in Godsil's book [3].

A similar orthogonality property then holds for the probabilistic counterparts Ah,Bh of the Pk,Qk polynomials. This provides alternative expressions of various probabilistic quantities in terms of scalar products involving the measure µ of (7). One can rederive in this way, via the combinatorial theory, a number of formulæ originally discovered by Karlin and McGregor. For instance, one has
pm,n(t)=pn ó
õ
¥


0
e-tx qm(x)qn(x)  dµ(x),
where the qk polynomials (closely related to the Bk and Qk) satisfy the recurrence lnqn+1 +(x-ln-µn)qn +µn qn-1=0.

4   So What?

The original motivation for the talk comes from the need to elucidate the behaviour of certain queueing systems in the context of telecommunication applications. For instance, the single server queue (M/M/1) is modelled by lj=r, µj=1, while the infinite server queue (M/M/¥) corresponds to lj=r, µj=j. (Models of population growth lead to considering different types of weights, like lj=(j+1)r, µj=j.) More specifically, the problem is to quantify parameters of some simple statistical multiplexing scheme that describe the quality of service on an ATM link. The relevant model is that of the M/M/¥ queue and parameters are to be analysed, like the duration q of an excursion above some level c, the volume V of lost information, or the number of bursts C in a busy period.

Each parameter leads to a specific continued fraction representation. By Theorem 2, the basic continued fraction of the M/M/¥ process is
1
s+r-
1r
s+1+r-
2r
·  
 · 
  ·
.
This is recognizable as an instance of Gauß's continued fraction associated to a quotient of contiguous hypergeometric functions. The numerator and denominator polynomials are the Poisson-Charlier polynomials that are orthogonal with respect to the Poisson measure.

The quantity V (area) leads to challenging asymptotics questions both for the M/M/¥ queue and for the M/M/1 queue. A simple modification of the basic techniques of this note shows that the bivariate Laplace transform with (s,u) ``marking'' (t,V) is obtained by the modified morphism,
c''(aj)=
lj
s+ju+lj+µj
,      c''(bj)=
µj
s+ju+lj+µj
.
In the case of area under the M/M/1 queue, quotients of continuous Bessel functions make an appearance. Stripped of its probabilistic context, the corresponding problem of tail estimation then admits a purely analytic formulation:
Problem. Let A(x) be a function whose Laplace transform is
~
A
 
(s)=
1
(s)1/2
J
 
n(s)+1
æ
ç
ç
è
2(r)1/2
s
ö
÷
÷
ø
J
 
n(s)
æ
ç
ç
è
2(r)1/2
s
ö
÷
÷
ø
,      n(s)=(1+r)/s,
with Jn a Bessel function, and r>0 a parameter. Show that, for some constants c1,c2, one has
ó
õ
¥


x
A(y)  dy ~ c1 x-1/4e
-c2(x)1/2
 
,      (x®+¥).
Under plausible analytic or probabilistic conjectures, precise (and useful!) quantitative conclusions can be drawn. See the papers by Guillemin and Pinchon [5, 6] for full developments.

References

[1]
Flajolet (P.). -- Combinatorial aspects of continued fractions. Discrete Mathematics, vol. 32, 1980, pp. 125--161.

[2]
Flajolet (P.), Françon (J.), and Vuillemin (J.). -- Sequence of operations analysis for dynamic data structures. Journal of Algorithms, vol. 1, 1980, pp. 111--141.

[3]
Godsil (C. D.). -- Algebraic Combinatorics. -- Chapman and Hall, 1993.

[4]
Goulden (Ian P.) and Jackson (David M.). -- Combinatorial Enumeration. -- John Wiley, New York, 1983.

[5]
Guillemin (Fabrice) and Pinchon (Didier). -- Continued fraction analysis of the duration of an excursion in an M/M/¥ system. Journal of Applied Probability, vol. 35, n°1, 1998, pp. 165--183.

[6]
Guillemin (Fabrice) and Pinchon (Didier). -- Excursions of birth and death processes, orthogonal polynomials, and continued fractions. -- Preprint, 1998. To appear in Journal of Applied Probability. 21 pages.

[7]
Henrici (Peter). -- Applied and Computational Complex Analysis. -- John Wiley, New York, 1977. 3 volumes.

[8]
Jones (William B.) and Magnus (Arne). -- Application of Stieltjes fractions to birth-death processes. In Saff (E. B.) and Varga (Richard S.) (editors), Padé and rational approximation. pp. 173--179. -- New York, 1977. Proceedings of an International Symposium held at the University of South Florida, Tampa, Fla., December 15-17, 1976.

[9]
Jones (William B.) and Thron (W. J.). -- Continued Fractions: Analytic Theory and Applications. -- Addison--Wesley, 1990, Encyclopedia of Mathematics and its Applications, vol. 11.

[10]
Karlin (S.) and McGregor (J. L.). -- The differential equations of birth-and-death processes, and the Stieltjes moment problem. Transactions of the American Mathematical Society, vol. 85, 1957, pp. 489--546.

[11]
Karlin (Samuel) and McGregor (James). -- The classification of birth and death processes. Transactions of the American Mathematical Society, vol. 86, 1957, pp. 366--400.

[12]
Karlin (Samuel) and Taylor (Howard). -- A First Course in Stochastic Processes. -- Academic Press, 1975, second edition.

[13]
Perron (Oskar). -- Die Lehre von der Kettenbrüchen. -- Teubner, 1954, vol. 2.
This document was translated from LATEX by HEVEA.