连续时间马尔可夫链(CTMC)

连续时间马尔可夫链(continuous-time Markov chain, CTMC)是一类应用广泛的随机过程。

一个CTMC 的定义包含两部分：

• 状态空间(state space)： ，即 的所有可能取值的集合。
• 转移率矩阵： ，刻画 在不同状态之间转移的“速率”。

假设当前 处于的状态为 ，那么过了很短一段时间 后，

即以概率 转移到其他状态 。

那么其呆在原地的概率则为：

我们设：

即为 从状态 转出的“速率”。

也可以把 理解为 的无穷小生成算符，即对任意的函数 :

的转移概率矩阵定义为：

即在时间 从 出发，在时间 落在 的概率。

CTMC最重要的是马尔可夫性(Markov property)，即未来CTMC的走势只与当前状态有关，而与历史无关。

据此我们可以推导出 满足的微分方程组。

首先注意到在 不依赖于时间的情况下， 只与 ， ， 有关，即：

这种情况下，我们只需要计算 。

根据马尔可夫性，我们有：

其中：

令 以及 为单位矩阵，我们就有：

那么：

这是一个常微分方程组。

同样的，通过：

我们也能得到：

据此我们可以得到转移概率矩阵：

如果初始状态的分布为 （行向量），其中，

则 时刻的状态分布为：

当 ，如果CTMC的分布达到一个稳态。应有：

令 ，即有：

矩阵 的第 行是从 状态出发最终形成的稳态分布，其可以通过求解以下线性方程组得到。

假设当前状态为 ，当前时间为 ，下次转移的等待时间设为 ，其分布可以通过状态转移的无记忆性推出。

即：

因此:

注意到 ，那么：

即 符合指数分布。

下一次转移之后的状态分布为：

生灭过程是一种特殊的CTMC，每次转移只能到达相邻状态。

在种群增长模型中，“生”代表种群数量增加一个，“灭”代表种群数量减少一个，“生”、“灭”的转移率由出生率和死亡率决定。

根据以上推导，我们可以模拟生灭过程的状态转移过程。

# 模拟一个生灭过程 (Birth-death process)

import numpy as np

b = 0.1 # 出生率 (per capita)
c = 0.05 # 死亡率 (per capita)

def birth_rate(N):
    return b*N # 总和出生率

def death_rate(N):
    return c*N # 总和死亡率

def simulation(N0, T):
    # 给定N0: 初始种群大小
    # 模拟时间T之前的种群增长
    
    t, N = 0, N0 # 初始化时间和种群大小
    
    tp = np.array([t]) # 记录时间
    Np = np.array([N]) # 记录不同时间的种群数量
    
    while t < T and N>1:
        # 产生下次出生或死亡的等待时间
        tau = np.random.exponential(1/(birth_rate(N)+death_rate(N)))
                
        t = t + tau # 更新当前时间
        
        # 判断此次事件具体是出生还是死亡
        U = np.random.uniform(size=1)
        if U < birth_rate(N)/(birth_rate(N)+death_rate(N)): # 出生
            N = N + 1
            
        else: # 死亡
            N = N - 1
        
        tp = np.append(tp, t) # 记录时间
        Np = np.append(Np, N) # 记录种群数量
    
    return tp, Np

# 实验
N0, T = 10, 10 # 初始数量和模拟时间跨度

tp, Np = simulation(N0, T) # 模拟

# 展示结果
import matplotlib.pyplot as plt

plt.step(tp, Np, where='post')
plt.plot(tp, Np, 'C0o', alpha=0.5)

plt.xlabel('time')
plt.ylabel('Population size')

plt.savefig('Birth-death.jpg')
plt.show()

一条模拟路径如下：