Matrix updating through time

One matrix multiplication gives the distribution after one step. Repeated multiplication follows the entire probability distribution through time without selecting random events.

The unique purpose of this lesson

Lesson 10 constructed a row-stochastic transition matrix \(P\) and calculated:

\[ \boldsymbol\pi_1=\boldsymbol\pi_0P. \]

We now repeat that update:

\[ \boldsymbol\pi_0 \longrightarrow \boldsymbol\pi_1 \longrightarrow \boldsymbol\pi_2 \longrightarrow\cdots. \]

The new tasks are to store every distribution, attach it to the correct time, verify probability conservation and interpret how probability flows towards extinction.

Continue with the same small SIS model

Use states \(\{0,1,2,3\}\) and the row-stochastic matrix:

\[ P= \begin{pmatrix} 1&0&0&0\\ 0.1&0.7&0.2&0\\ 0&0.2&0.6&0.2\\ 0&0&0.3&0.7 \end{pmatrix}. \]

The initial condition remains certain at \(I_0=1\):

\[ \boldsymbol\pi_0=(0,1,0,0). \]

Keeping the model unchanged lets us focus entirely on repeated distribution updating.

Two equivalent calculation methods

Iterative updating

\[ \boldsymbol\pi_{n+1} =\boldsymbol\pi_nP. \]

Apply one update inside a loop and store each result. This is convenient when every intermediate distribution is required.

Matrix powers

\[ \boldsymbol\pi_n =\boldsymbol\pi_0P^n. \]

Calculate the distribution at a particular step directly from the \(n\)-step transition matrix \(P^n\).

These methods must agree. For example:

\[ \boldsymbol\pi_2 =(\boldsymbol\pi_0P)P =\boldsymbol\pi_0P^2. \]

What does \(P^n\) mean?

The entry \((P^n)_{ij}\) is the probability of being in state \(j\) after \(n\) steps when the chain begins in state \(i\).

\(P^n\) does not mean raising every entry of \(P\) separately to the power \(n\). It means multiplying the matrix by itself \(n\) times:

\[ P^3=P\cdot P\cdot P. \]

In Python, np.linalg.matrix_power(P, n) performs matrix powers.

Store all distributions correctly

distributions = [pi_current.copy()]

for step in range(number_of_steps):
    pi_current = pi_current @ P
    distributions.append(pi_current.copy())

Initial storageStore \(\boldsymbol\pi_0\) before updating.

MultiplyCalculate one new distribution.

ReplaceThe result becomes the current distribution.

AppendStore a separate copy for this time step.

Why use `.copy()`?

A NumPy array is a mutable object: its contents can be changed. pi_current.copy() creates a separate array containing the current values.

In this particular loop, pi_current @ P creates a new array, so appending without .copy() would also work. We use .copy() to state the intention clearly and to keep storage safe if the code is later changed to update arrays in place.

Attach steps to physical time

The matrix advances the process by one fixed interval \(\Delta t\). At step \(n\):

\[ t_n=n\Delta t. \]

With \(\Delta t=0.5\) day, step 10 corresponds to day 5. The matrix exponent counts steps, whereas the time axis records biological time.

Probability conservation through time

Every stored distribution must satisfy:

\[ \sum_{i=0}^{N}\pi_{n,i}=1. \]

For a row-stochastic matrix:

\[ \boldsymbol\pi_{n+1}\mathbf 1 = \boldsymbol\pi_nP\mathbf 1 = \boldsymbol\pi_n\mathbf 1 =1, \]

because \(P\mathbf 1=\mathbf 1\) expresses that every row of \(P\) sums to 1.

Computers store decimals approximately, so calculated sums may differ from 1 by a tiny rounding amount. Use np.allclose rather than requiring every floating-point sum to equal exactly 1.

How probability flows into extinction

State 0 is absorbing because row 0 is \((1,0,0,0)\). Once probability reaches state 0, later updates keep it there.

\[ \Pr(I_n=0)=\pi_{n,0}. \]

Therefore the probability of being extinct by step \(n\) cannot decrease with \(n\).

The other state probabilities need not change monotonically. Probability can enter and leave transient states 1, 2 and 3 during different updates.

An exact mean from the distribution

Once the full distribution is known, the mean infectious count at step \(n\) is:

\[ \mathbb E[I_n] = \sum_{i=0}^{N}i\pi_{n,i}. \]

In Python, the dot product distributions @ states calculates this weighted average at every stored step.

The mean need not be a possible individual state. For example, \(\mathbb E[I_n]=0.8\) does not mean that a trajectory contains 0.8 of a person. It averages over the probabilities of integer-valued states.

Exact distribution versus one trajectory

Matrix result

At every step, it gives probabilities for all four states and the exact mean implied by the finite-state model.

One sample path

At every step, it occupies one integer state. It may differ considerably from the distribution mean.

The code generates one illustrative path using the same rows of \(P\). It is not used to estimate the matrix probabilities. Estimation from many independent trajectories is the purpose of Lesson 12.

Interactive Python laboratory

Click Run code to update the exact distribution for 40 steps, compare iteration with \(P^n\), and display state probabilities, extinction probability, probability conservation, the exact mean and one sample path.

Interactive PythonExact DTMC distribution through time

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Same row-stochastic SIS matrix as Lesson 10
P = np.array([
    [1.0, 0.0, 0.0, 0.0],
    [0.1, 0.7, 0.2, 0.0],
    [0.0, 0.2, 0.6, 0.2],
    [0.0, 0.0, 0.3, 0.7]
])

states = np.arange(4)
dt = 0.5
number_of_steps = 40

# Certain initial state I_0 = 1
pi_initial = np.array([0.0, 1.0, 0.0, 0.0])
pi_current = pi_initial.copy()

# Validate P and the initial distribution
if not np.allclose(P.sum(axis=1), 1.0):
    raise ValueError("Rows of P must sum to 1.")

if (P < -1e-12).any() or (P > 1 + 1e-12).any():
    raise ValueError("P contains an invalid probability.")

if not np.isclose(pi_initial.sum(), 1.0):
    raise ValueError("Initial probabilities must sum to 1.")

# Store pi_0, pi_1, ..., pi_number_of_steps
distributions = [pi_current.copy()]

for step in range(number_of_steps):
    pi_current = pi_current @ P
    distributions.append(pi_current.copy())

distributions = np.vstack(distributions)
step_values = np.arange(number_of_steps + 1)
time_values = step_values * dt

# Validate every stored distribution
probability_sums = distributions.sum(axis=1)

if not np.allclose(probability_sums, 1.0):
    raise ValueError("Probability was not conserved.")

if (distributions < -1e-12).any():
    raise ValueError("A negative state probability appeared.")

# Compare iteration with a matrix power at step 10
chosen_step = 10
pi_from_iteration = distributions[chosen_step]
pi_from_power = (
    pi_initial
    @ np.linalg.matrix_power(P, chosen_step)
)

if not np.allclose(pi_from_iteration, pi_from_power):
    raise ValueError("Iteration and matrix power disagree.")

# Derived exact quantities
extinction_probability = distributions[:, 0]
expected_I = distributions @ states
maximum_sum_error = np.max(np.abs(probability_sums - 1))

# Selected distributions for a readable table
selected_steps = [0, 1, 5, 10, 20, 40]
selected_table = pd.DataFrame(
    distributions[selected_steps],
    columns=["Pr(I=0)", "Pr(I=1)", "Pr(I=2)", "Pr(I=3)"]
)
selected_table.insert(0, "time", np.array(selected_steps) * dt)
selected_table.insert(0, "step", selected_steps)

print("Selected exact distributions")
print(selected_table.round(6).to_string(index=False))
print()
print(f"Maximum probability-sum error: "
      f"{maximum_sum_error:.3e}")
print()
print(f"Step {chosen_step} from iteration:")
print(np.round(pi_from_iteration, 8))
print(f"Step {chosen_step} from P^{chosen_step}:")
print(np.round(pi_from_power, 8))

# One illustrative trajectory from the same matrix
seed = 7
rng = np.random.default_rng(seed)
current_state = 1
sample_path = [current_state]

for step in range(number_of_steps):
    current_state = rng.choice(
        states,
        p=P[current_state]
    )
    sample_path.append(current_state)

# Figure 1: all exact state probabilities
plt.figure(figsize=(9, 5.2))
for state in states:
    plt.plot(
        time_values,
        distributions[:, state],
        marker="o",
        markersize=3,
        label=f"Pr(I={state})"
    )
plt.xlabel("Time (days)")
plt.ylabel("Probability")
plt.title("Exact SIS state probabilities through time")
plt.ylim(0, 1.03)
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

# Figure 2: absorption probability and conservation
fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))

axes[0].plot(
    time_values,
    extinction_probability,
    color="#dc2626",
    linewidth=2
)
axes[0].set_xlabel("Time (days)")
axes[0].set_ylabel("Pr(I = 0)")
axes[0].set_title("Probability of extinction by time")
axes[0].set_ylim(0, 1.03)
axes[0].grid(alpha=0.25)

axes[1].plot(
    time_values,
    probability_sums,
    color="#16a34a",
    linewidth=2
)
axes[1].axhline(1, color="black", linestyle="--")
axes[1].set_xlabel("Time (days)")
axes[1].set_ylabel("Sum of state probabilities")
axes[1].set_title("Probability conservation")
axes[1].set_ylim(0.999999, 1.000001)
axes[1].grid(alpha=0.25)

fig.tight_layout()

# Figure 3: exact mean versus one random trajectory
plt.figure(figsize=(9, 4.8))
plt.plot(
    time_values,
    expected_I,
    color="#7c3aed",
    linewidth=2.5,
    label="Exact mean E[I]"
)
plt.step(
    time_values,
    sample_path,
    where="post",
    color="#f59e0b",
    alpha=0.8,
    label=f"One trajectory, seed {seed}"
)
plt.xlabel("Time (days)")
plt.ylabel("Infectious count")
plt.title("A distribution mean is not an individual path")
plt.yticks(states)
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

Output

Run the code to see the result.

Understand the new code

Code	Meaning
`distributions = [pi_current.copy()]`	Begins storage with the initial distribution at step 0.
`np.vstack(distributions)`	Stacks the stored one-dimensional vectors as rows of one two-dimensional array.
`distributions[:, 0]`	Selects every time row and column 0, giving extinction probability through time.
`np.linalg.matrix_power(P, 10)`	Calculates the 10-step transition matrix \(P^{10}\).
`distributions @ states`	Calculates the probability-weighted mean infectious count at every step.
`np.max(np.abs(...))`	Finds the largest absolute probability-conservation error.
`rng.choice(states, p=P[current_state])`	Selects one next state using the probability row for the current state. It is a compact alternative to writing cumulative intervals explicitly.

Interpret the long-time pattern carefully

For this finite SIS chain, state 0 is absorbing and every positive state can eventually reach it. Consequently, extinction probability approaches 1 as the number of steps becomes large.

\[ \boldsymbol\pi_n \longrightarrow (1,0,0,0) \qquad\text{as }n\to\infty. \]

At a finite time such as day 20, the extinction probability may still be below 1. “Eventually certain” does not mean “already certain by the current observation time.”

Try these experiments

Increase number_of_steps to 100. Observe the approach to the absorbing distribution.
Change pi_initial to [0, 0, 1, 0]. Compare early extinction probabilities.
Use the mixed initial distribution [0.2, 0.5, 0.3, 0].
Change the trajectory seed. Confirm that the exact matrix probabilities and mean do not change.
Change chosen_step. Verify that iterative updating and the matching matrix power still agree.
Temporarily change a row of \(P\) so it does not sum to 1. Read the validation error.

What this lesson has—and has not—done

You can now:

repeat \(\boldsymbol\pi_{n+1}=\boldsymbol\pi_nP\) through time;
store and validate every exact state distribution;
use \(\boldsymbol\pi_n=\boldsymbol\pi_0P^n\);
follow extinction probability and probability conservation;
calculate the exact mean from a finite distribution; and
distinguish an exact distribution mean from one random path.

Not yet: one illustrative trajectory does not estimate probabilities. Lesson 12 simulates many independent trajectories and constructs empirical distributions for models whose full transition matrix may be too large to update conveniently.

Matrix updating through time

The unique purpose of this lesson

Continue with the same small SIS model

Two equivalent calculation methods

Iterative updating

Matrix powers

What does \(P^n\) mean?

Store all distributions correctly

Why use .copy()?

Attach steps to physical time

Probability conservation through time

How probability flows into extinction

An exact mean from the distribution

Exact distribution versus one trajectory

Matrix result

One sample path

Interactive Python laboratory

Understand the new code

Interpret the long-time pattern carefully

Try these experiments

What this lesson has—and has not—done

Why use `.copy()`?