g

gqzhang

V1

2022/12/06阅读:104主题:默认主题

连续时间马尔可夫链(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 = 1010 # 初始数量和模拟时间跨度

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()

一条模拟路径如下:

分类:

数学

标签:

数学

作者介绍

g
gqzhang
V1