MODELLING QUASISTATIONARITY IN EVANESCENT PROCESSES
by
Phil Pollett
Department of Mathematics
The University of Queensland
(Download the following image (as postscript)):
Figure 1. Quasistationarity of an explosive population process
[Simulation |
java code]
QUASISTATIONARITY
Processes appear to be stationary over any reasonable time scale, yet they are evanescent.
Examples:
AN AUTO-CATALYTIC REACTION
Consider the following reaction scheme:
where X is a catalyst. Suppose that there are two stages, namely
Let number of X molecules at time t.
Suppose that the concentration of A is held constant; let a be the number of molecules of A. The state space is and the transition rates, , are given by
(Download the following image (as postscript)):
Figure 2. Auto-catalytic reaction simulation
[Simulation |
java code]
AN EPIDEMIC MODEL
(Ridler-Rowe)
The state at time t:
- number of susceptibles
- number of infectives
State space:
Transition rates :
if , then
(Download the following image (as postscript)):
Figure 3. State transition diagram for
Ridler-Rowe epidemic model.
Clearly,
is an irreducible transient class.
Ridler-Rowe (1967): Q is regular (nonexplosive) and absorption occurs with probability 1.
(Download the following image (as postscript)):
Figure 4. Simulation of the Ridler-Rowe epidemic model
[Simulation |
java code]
(Download the following image (as postscript)):
Figure 5. Simple contour plot of the usual diffusion
approximation for the Ridler-Rowe epidemic model
(to be imagined as superimposed on Figure 4)
(Download the following image (as postscript)):
Figure 6. Surface plot of the usual diffusion
approximation for the Ridler-Rowe epidemic model
QUASISTATIONARY DISTRIBUTIONS
Notation. Let be the state of the system at time t and suppose that , where S, the state space, is a discrete set. Let be the matrix of transition rates (q-matrix), so that , for , represents the transition rate from state i to state j and , where represents the transition rate out of state i.
Assumptions. Take 0 to be the sole absorbing state, that is, , and suppose that , where is an irreducible class. In order that there be a positive probability of ever reaching 0, given that the process starts in C, suppose that for at least one .
Let , .
Definition. Let be a probability distribution over C and let
Then, is a quasistationary distribution if
That is, if the process has as its initial distribution, then is a quasistationary distribution if the state probabilities at time t, conditional on the process being in C at t, are the same for all t.
If is a quasistationary distribution then, under certain conditions,
no matter what the value of .
A CHARACTERIZATION
If is a quasistationary distribution, then clearly
where . It is easy to show that g satisfies:
and . Thus, , for some . Since the converse is also true, we have the following result:
Proposition. A probability distribution, , over C is a quasistationary distribution if and only if, for some , is a -invariant measure, that is
HOW DO WE DETERMINE from Q?
Rewrite (1) as
and remember that is the right-hand derivative of near 0. So, on dividing by t and letting , we get (formally)
or, equivalently,
Accordingly, we shall say that is a -invariant measure for Q whenever (2) holds.
Proposition. If is a quasistationary distribution then, for some , is a -invariant measure for Q.
IS THE CONVERSE TRUE?
Suppose that, for some , is a -invariant measure for Q, that is
Is a quasistationary distribution?
Sum this equation over : we get (formally)
Theorem. Let be a probability distribution over C and suppose that is a -invariant measure for Q. Then, with equality if and only if is a quasistationary distribution.
AN EXAMPLE
The birth-death-catastrophe process. Let and suppose that
where , a>0, for at least one value of and . Thus, jumps occur at a constant ``per capita'' rate and, at a jump time, a birth occurs with probability a, or otherwise a catastrophe occurs, the size of which is determined by the probabilities , .
Clearly, 0 is an absorbing state and is an irreducible class.
So, does the process admit a quasistationary distribution?
Substitute the transition rates in (2):
and, for ,
If we try a solution of the form , the first equation tells us that , where
and, on substituting both of these quantities in the second equation, we find that . This latter equation has a unique solution, , on . Thus, by setting we can obtain a positive solution, , to (2), such that whenever . Further, it is easy to show that .
The condition is satisfied only in the subcritical case, that is, when the drift, D, given by
is strictly negative; incidentally, this also guarantees that absorption occurs with probability 1.
Proposition. The subcritical birth-death-catastrophe process has a geometric quasistationary distribution , , given by
where is the unique solution to on the interval .
COMPUTATIONAL METHODS
Finite S : Mandl (1960) showed that the restriction of Q to C has eigenvalues with negative real parts and the one with maximal real part (called above) is real and has multiplicity 1, and, the corresponding left eigenvector, , has positive entries; this is, of course, a -invariant measure for Q (unique up to constant multiples). Since S is finite, the quasistationary distribution, , exists and is given by
Infinite S : Truncate the restriction of Q to an matrix, , and construct a sequence, , of eigenvectors and hope that this converges to a -invariant measure, m, for Q, et cetera.
HOW SHOULD WE EVALUATE m?
Consider the Ridler-Rowe epidemic model.
First truncate C to
and restrict Q to . Use the transformation to convert the restricted q-matrix into an matrix, , where .
Evaluation of the eigenvectors of Q is not completely trivial. For example, if (as well shall assume) N=100, that is , Q needs 400 Mbytes of storage.
However, for N large, this matrix is large and sparse. Does this help?
THE ARNOLDI METHOD
We need to solve , where .
Using an initial estimate of x, the basic Arnoldi method produces an (upper-Hessenberg) matrix and an matrix with
and such that if is an eigenvector of , then, for m large, is close to an eigenvector of A.
Solve for using standard (dense-matrix) methods.
AN ITERATIVE ARNOLDI METHOD
Take m small (we found that m=20 worked best). Then, using an initial estimate, , of the eigenvector x, apply the basic Arnoldi method (to obtain and ) and set to be the dominant eigenvalue of if this is real, or set equal to zero otherwise.
Now solve
with z chosen at random and repeat the procedure with a new initial estimate, given by
Continue until the residual, , is sufficiently small.
(Download the following image (as postscript)):
Figure 7. Surface plot of the quasistationary distribution
for the Ridler-Rowe epidemic model
(Download the following image (as postscript)):
Figure 8. Contour plot of the quasistationary distribution
for the Ridler-Rowe epidemic model
(Download the following image (as postscript)):
Figure 8. Comparison between the quasistationary distribution
and the diffusion approximation (contour plot) for the
Ridler-Rowe epidemic model
(Download the following image (as postscript)):
Figure 10. Graph of the residual norm versus the
number of outer iterations for the
Iterative Arnoldi Method
(Download the following image (as postscript)):
Figure 11. Graph of the residual norm versus the
number of outer iterations for various
values of m, the number of inner iterations