基于Python的AGCT碱基转移矩阵概率计算问题求助

问题分析与修正说明 #

问题1：20步后处于各状态的概率 #

原代码这部分逻辑正确。通过初始分布向量与转移矩阵的20次幂相乘，得到的就是20步后的状态概率分布。

问题2：20步全程处于状态G的概率 #

原代码错误地使用了(np.linalg.matrix_power(P,20)[1,1])**20，这是完全错误的。正确逻辑是：从初始状态G出发，每一步都必须保持在G，即每次转移都是G→G，共20次转移。每次转移的概率为pii=0.97，因此总概率为pii ** 20。

问题3：20步后处于G且过程中至少访问过其他状态的概率 #

原代码的思路正确（用20步后处于G的概率减去全程处于G的概率），但由于问题2的计算错误，导致结果不正确。修正问题2的计算后即可得到正确结果。

修正后的代码 #

import numpy as np

# 转移矩阵定义
pij = 0.01
pii = 0.97
P = np.array([[pii, pij, pij, pij], 
              [pij, pii, pij, pij], 
              [pij, pij, pii, pij], 
              [pij, pij, pij, pii]])

# 初始概率分布：从状态G出发（索引1）
init_dist = np.array([0, 1, 0, 0])

# 问题1：计算20步后各状态的概率
prob_20_steps = np.linalg.matrix_power(P, 20) @ init_dist
print("20步后各状态概率：")
print(f"P(A) = {prob_20_steps[0]:.6f}")
print(f"P(G) = {prob_20_steps[1]:.6f}")
print(f"P(C) = {prob_20_steps[2]:.6f}")
print(f"P(T) = {prob_20_steps[3]:.6f}")

# 问题2：20步全程处于状态G的概率
prob_stay_g_all_steps = pii ** 20
print(f"\n20步全程处于G的概率：{prob_stay_g_all_steps:.6f}")

# 问题3：20步后处于G且过程中至少访问过其他状态的概率
prob_g_with_other_visited = prob_20_steps[1] - prob_stay_g_all_steps
print(f"20步后处于G且过程中至少访问过其他状态的概率：{prob_g_with_other_visited:.6f}")

运行结果说明 #

执行修正后的代码，会得到如下结果（保留6位小数）：

20步后各状态概率：
P(A) = 0.139497
P(G) = 0.581510
P(C) = 0.139497
P(T) = 0.139497

20步全程处于G的概率：0.543794
20步后处于G且过程中至少访问过其他状态的概率：0.037716

• 问题1的结果符合马尔可夫链的收敛趋势：随着步数增加，各状态概率逐渐向均匀分布（0.25）靠近，但20步后仍保留初始状态G的较高概率。
• 问题2的结果是连续20次保持G状态的概率，由于每次保持的概率较高（0.97），该值仍较大。
• 问题3的结果是离开G后又返回G的概率，数值较小，符合预期。

内容的提问来源于stack exchange，提问作者Allisa Newmann