关于连续时间马尔可夫链拟合与模拟流程及代码实现的正确性验证问询

嘿，先给你点个赞！能从数据模拟到CTMC的拟合、模拟全流程自己动手实现，这个学习劲头太值得肯定了。咱们一步步拆解你的思路和代码，看看整体框架的正确性，以及需要调整的细节：

一、数据模拟部分 #

你用R模拟的CTMC观测数据是符合逻辑的：

• 确保了无自环转移（每次新状态都和当前状态不同），这完全符合CTMC的特性（在连续时间里，过程不会在同一时刻停留在原状态）；
• 用指数分布模拟状态停留时间，这也是CTMC的核心假设（停留时间服从指数分布，无记忆性）。

对应的模拟代码整理如下：

# Define the possible states
states <- c("S1", "S2", "S3")

# Initialize the sequence of state transitions with the first state
transitions <- "S1"

# Loop to generate 99 more state transitions
for (i in 2:100) {
    # Get the current state
    current_state <- transitions[i-1]
    # Generate a new state that is different from the current state
    new_state <- current_state
    while (new_state == current_state) {
        new_state <- sample(states, 1)
    }
    # Add the new state to the sequence of state transitions
    transitions <- c(transitions, new_state)
}

# Set the rate parameter of the exponential distribution
lambda <- 0.5
# Simulate 100 observations from the exponential distribution
times <- rexp(100, rate = lambda)

# final data: no consecutive identical states (compliant with CTMC)
my_data = data.frame(id = 1:100, transitions, times)

二、CTMC拟合步骤的正确性分析 #

你的整体框架完全符合CTMC的拟合逻辑，但有两个关键细节需要修正：

1. 嵌入离散马尔可夫链概率计算（Step 1） #

你的思路是对的：$p_{ij} = \frac{n_{ij}}{n_i}$，其中$n_{ij}$是从状态i到j的转移次数，$n_i$是从状态i的总转出次数。

但代码里的exit_counts计算有误：

# 你的错误代码：
exit_counts <- colSums(count_matrix)
# 正确应该是行求和：因为count_matrix的行是出发状态，列是到达状态
exit_counts <- rowSums(count_matrix)

修正后，transition_probs的计算才是正确的$p_{ij}$。

2. 状态离开率$\lambda_i$的计算（Step 2） #

你对MLE的公式理解反了！指数分布率参数的MLE是转出次数除以总停留时间，即$\lambda_i = \frac{n_i}{t_i}$，其中$n_i$是从状态i的转出次数，$t_i$是在状态i的总停留时间。

你的代码里写反了，修正如下：

# 你的错误代码：
lambda_1 <- sum(my_data$times[my_data$transitions == "S1"]) / nrow(my_data[my_data$transitions == "S1", ])
# 正确代码：
n_1 <- exit_counts["S1"]  # 从S1的转出次数
t_1 <- sum(my_data$times[my_data$transitions == "S1"])
lambda_1 <- n_1 / t_1

# 同理修正lambda_2和lambda_3：
n_2 <- exit_counts["S2"]
t_2 <- sum(my_data$times[my_data$transitions == "S2"])
lambda_2 <- n_2 / t_2

n_3 <- exit_counts["S3"]
t_3 <- sum(my_data$times[my_data$transitions == "S3"])
lambda_3 <- n_3 / t_3

3. 速率矩阵Q的构建（Step 3） #

你的思路完全正确：

• 当$i \neq j$时，$Q_{ij} = \lambda_i p_{ij}$
• 对角线元素$Q_{ii} = -\sum_{k \neq i} Q_{ik}$（确保每行和为0）

只要修正了Step 1和Step 2的错误，这部分代码就能正确运行。

4. 时间依赖转移概率计算（Step 4） #

你用expm包计算矩阵指数$P(t) = e^{Qt}$，再通过初始分布$p_0$计算t时刻的状态分布，这个方法是CTMC中最常用的数值解法，完全正确。

对应的模拟和可视化代码也没问题，整理如下：

library(Matrix)
library(expm)
library(ggplot2)

# 初始分布：从S1出发，概率为1
p0 <- c(1, 0, 0)
t <- 5.4
P_t <- expm(rate_matrix * t)
p_t <- p0 %*% P_t

# 绘制状态概率随时间变化的曲线
t_range <- seq(0, 100, length.out=1000)
p_t <- t(sapply(t_range, function(t) p0 %*% expm(rate_matrix * t)))
df <- data.frame(t = t_range, p1 = p_t[,1], p2 = p_t[,2], p3 = p_t[,3])

ggplot(df, aes(x = t)) +
    geom_line(aes(y = p1, color = "State 1")) +
    geom_line(aes(y = p2, color = "State 2")) +
    geom_line(aes(y = p3, color = "State 3")) +
    labs(x = "Time", y = "State Probability", color = "State")

三、平稳分布计算的正确性 #

你通过寻找速率矩阵Q的特征值为0的特征向量来计算平稳分布，这个思路是对的：CTMC的平稳分布$\pi$满足$\pi Q = 0$，对应Q的转置矩阵特征值为0的特征向量。

不过需要注意：特征向量可能包含虚部，需要取实部后再归一化，修正后的代码如下：

# Calculate the eigenvectors and eigenvalues of the rate matrix's transpose
eig <- eigen(t(rate_matrix))
# Extract the eigenvector corresponding to the eigenvalue closest to zero
pi_complex <- eig$vectors[, which.min(abs(eig$values))]
# Take real part and normalize
pi <- Re(pi_complex)
pi <- pi / sum(pi)

总结 #

你的整体思路框架完全正确，只要修正上述两个关键细节（exit_counts的计算、lambda_i的MLE公式），整个流程就能正确拟合CTMC并模拟状态轨迹啦！

备注：内容来源于stack exchange，提问作者stats_noob