如何修正R语言ape包matexpo生成的概率矩阵中的负数值

这个报错是数值计算误差导致的：ape包内置的matexpo函数采用泰勒展开近似计算矩阵指数，当你使用高维稀疏有序转换矩阵（≥8状态，仅允许相邻状态转换）、且转移速率与分支长度的乘积落在近似误差区间时，会输出极小的负数值（通常是1e-16级别，属于浮点计算误差，不是真实的概率为负），这些负值被传入sample.int就会触发报错。

• 方案1：直接截断负概率并重归一化
  这是成本最低的修复方式，对模拟结果的影响可忽略。你可以自行封装模拟逻辑，对生成的转移概率矩阵做后处理：

  library(ape)
  data("bird.orders")
  # 你的原始模型矩阵
  model.matrix <- matrix(c(0,0.1,0,0,0,0,0,0,
                       0.1,0,0.1,0,0,0,0,0,
                       0,0.1,0,0.1,0,0,0,0,
                       0,0,0.1,0,0.1,0,0,0,
                       0,0,0,0.1,0,0.1,0,0,
                       0,0,0,0,0.1,0,0.1,0,
                       0,0,0,0,0,0.1,0,0.1,
                       0,0,0,0,0,0,0.1,0), 8)
  # 自定义概率矩阵修正函数
  fix_pmat <- function(p) {
    p[p < 0] <- 0 # 截断所有负值
    p <- p / rowSums(p) # 每行重归一化到和为1
    return(p)
  }
  # 构造Q矩阵
  k <- ncol(model.matrix)
  freq <- rep(1/k, k)
  Q <- model.matrix * rep(freq, each = k)
  diag(Q) <- 0
  diag(Q) <- -rowSums(Q)
  # 调用修正逻辑模拟性状
  sim_trait <- rTraitDisc(bird.orders, model = function(t) fix_pmat(matexpo(Q*t)), k = k, root.value = 1)
  

• 方案2：降低转移速率的基数
  你当前模型矩阵的非对角元为0.1，把这个值调小（比如改为0.05），可以缩小Q矩阵与分支长度的乘积，避免matexpo的近似误差出现负值：

  # 调整后的低速率模型矩阵
  model.matrix <- matrix(c(0,0.05,0,0,0,0,0,0,
                       0.05,0,0.05,0,0,0,0,0,
                       0,0.05,0,0.05,0,0,0,0,
                       0,0,0.05,0,0.05,0,0,0,
                       0,0,0,0.05,0,0.05,0,0,
                       0,0,0,0,0.05,0,0.05,0,
                       0,0,0,0,0,0.05,0,0.05,
                       0,0,0,0,0,0,0.05,0), 8)
  # 直接调用原生rTraitDisc即可正常运行
  sim_trait <- rTraitDisc(phy = bird.orders, model = model.matrix)
  

• 方案3：替换更精准的矩阵指数计算工具
  改用expm包的expm函数计算矩阵指数，其数值稳定性远高于ape自带的matexpo，不会出现负值：

  library(expm)
  # 沿用你原始的Q矩阵，替换矩阵指数计算逻辑
  sim_trait <- rTraitDisc(bird.orders, model = function(t) expm(Q*t), k = k, root.value = 1)
  

• 方案4：改用其他支持有序模型的模拟工具
  推荐用phytools包的sim.Mk函数，该函数原生处理了这类数值误差，对有序性状模型的适配性更好：

  library(phytools)
  # 直接定义8状态有序等速率模型，无需手动构造转换矩阵
  sim_trait <- sim.Mk(tree = bird.orders, model = "ORDERED", k = 8, rate = 0.1)
  

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