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