probability distribution (in this case rate = 2.0) which we use to model the number of
they cannot have any progeny. A bacterial population of one would have the following for
the first ten transition probabilities, T
1,0
to T
1,10
: 0.1353, 0.2707, 0.2707, 0.18046,
0.09028, 0.0361, 0.0120, 0.003, 0.0008, 0.0001, assuming a Poisson distribution with rate
2.0. The calculation of transition probabilities for a population of two (T
2,j
) would take
more work, because you have separate outcomes for each individual in the population and
these would have to be multiplied for the total number of combinations. However,
simulating population growth with this kind of model you fortunately do not need to know
the exact value of T
j,k
for all states j and k but instead consider individuals by themselves
and calculate the number of progeny for each.
We can simulate the whole process in Python, given a particular number of individuals
in one generation, and with a Poisson distribution randomly determining how many
progeny one bacterium will have. We can calculate the number of bacteria in the next
generation by considering how many offspring each of the individuals in the current
generation has. We define a function to generate the size of the population using a
scipy.stats random variable object. Note that we use the .rvs() method to draw the required
number of random samples from this distribution, which in this case would be the
numbers of progeny for each bacterium. The population of the next generation is then
simply the summation of the numbers of progeny.
def getNextGenPop (currentPop, randVar):
progeny = randVar.rvs(size=currentPop)
nextPop = progeny.sum()
return nextPop
A random variable object with a Poisson distribution is generated with the required rate
(average number of progeny in the generation):
from scipy.stats import poisson
rate = 1.02 # A deliberately low rate
poissRandVar = poisson(rate)
The Markov process is then simulated by repeatedly recalculating the populations for
the subsequent generations. At the low reproductive rate which we use here, you will note
that the population sometimes dies out.
pop = 25
for i in range(100):
pop = getNextGenPop(pop, poissRandVar)
print("Generation:%3d Population:%7d" % (i, pop))
Returning to the general situation, the transition matrix determines the properties of the
Markov chain. For example, if we define the initial distribution I of the probability of the
state j to be:
I
j
= Pr(State
0
= j)
then the probability for each state at the next position in the chain is the summation over
all the possible starting states multiplied by the conditional probability to get to the next
state. In other words we can get to the next state via multiple routes from different initial
states, each with a potentially different probability:
Pr(State
1
= j) = Pr(State
1
= j | State
0
= k) Pr(State
0
= k) =
∑
k
I
k
T
k,j
Thus the distribution at chain position 1 can be determined from that at chain position 0
just by multiplication by the transition matrix (T). It can also be shown, for example, that:
Pr(State
n
= k | State
0
=j) = (T
n
)
j,k
where T
n
is the nth power of the transition matrix. Thus in theory we can exactly calculate
the distribution of states at any subsequent position (State
n
), given the starting state (at
position zero) and multiplying by the matrix T
n
. Accordingly we have:
Pr(State
n
= j) = Pr(State
n
= j | State
0
= k) Pr(State
0
= k) =
∑
k
I
k
(T
n
)
k,j
This also allows us to calculate the probability of a whole sequence of specified states
(from i
0
to i
n
) from the starting state and the product of the individual transition matrix
elements which are selected from the knowledge of the states:
Pr(State
0
= i
0
, State
1
= i
1
, …, State
n
= i
n
) = I
i
0
T
i
0
i
1
T
i
n-1
i
n
Lastly, analyses of Markov chains sometimes refer to the equilibrium distribution. An
equilibrium distribution, π
i
, for a Markov chain is a distribution of states which does not
vary as the chain evolves for subsequent probabilistic trials, i.e. applying the transition
matrix regenerates the previous distribution (πT = π). If an equilibrium distribution exists,
then it can be shown, under weak assumptions, that all distributions approach the
equilibrium distribution as the number of trials gets large, i.e. I × T
n
approaches π as n
gets large, for any starting state. This can be important for biological and predictive
systems where the equilibrium distribution represents a long-term description of a system
when it has had time to ‘settle down’.
Do'stlaringiz bilan baham: 
