枚举R中顺序概率树中的路径子集

Question

为了说明问题，让我们定义以下矩阵（其中NA表示该选项在t期间不可用）

set.seed(1)
x <- matrix(NA, 4, 4, dimnames = list(paste0("t=", seq_len(4)), LETTERS[seq_len(4)]))
x[lower.tri(x, diag = TRUE)] <- rnorm(10)

这给出了一个如下所示的矩阵：

              A           B          C         D
t=1  0.91897737          NA         NA        NA
t=2  0.78213630  0.61982575         NA        NA
t=3  0.07456498 -0.05612874 -1.4707524        NA
t=4 -1.98935170 -0.15579551 -0.4781501 0.4179416

我们的目标是计算每个值在每个时间段 $t$ 中最高的概率，但是，这些值取决于前几个时间段的值。例如，在从周期 t=2 移动到 t=3 并假设 A 最高时，A 仅与 C 相比，而不是 B，因为在 t=2 中假定为更高。我们可以将问题构造成这样的树：

所以对于t=1，概率为 1，对于t=2，我们从 1 个分组中计算 2 个概率，在 t=3 中，我们从 2 个分组中计算 4 个概率（注意一个选项是如何从比较中排除的，因为顺序在t-1) 和t=4 中它不是最高的依赖和固有假设，我们从 4 个分组中计算了 8 个概率。然后，最终概率是构成 8 条路径的每个 t 中概率的乘积。在真正的问题中，t 变得很大，手动识别这些分组变得不可行。

我一直在尝试想出一种聪明的方法来识别这些路径并计算概率。一个想法是为每个可能的模式使用一组“屏蔽矩阵”。这样我就可以简单地将掩码矩阵相乘并执行行操作。但是，随着级别数量的增加，我无法找到一种可靠的方法来填充不同的掩码矩阵。

例如，假设在直到最后一个周期的所有周期中选择A 的模式可以用以下掩码矩阵来描述：

mask <- matrix(c(
1, NA, NA, NA,
1, 1,  NA, NA,
1, NA, 1,  NA,
1, NA, NA, 1
), ncol = 4, byrow = TRUE, dimnames = list(paste0("t=", seq_len(4)), LETTERS[seq_len(4)]))

看起来像这样（在这种情况下，4 个可能的比较之一）：

    A  B  C  D
t=1 1 NA NA NA
t=2 1  1 NA NA
t=3 1 NA  1 NA
t=4 1 NA NA  1

并且我们可以像这样计算每个时期的概率（所有行的总和应为 1）：

exp_x <- exp(x * mask)
sum_exp_x <- rowSums(exp_x, na.rm = TRUE)
pr_x <- exp_x / sum_exp_x

             A         B         C         D
t=1 1.00000000        NA        NA        NA
t=2 0.54048879 0.4595112        NA        NA
t=3 0.82423638        NA 0.1757636        NA
t=4 0.08261824        NA        NA 0.9173818

对于所有可能的路径，如tgrows，是否有一种聪明的方法？还是填充一组掩码矩阵以循环遍历的好方法？我试图避免这个问题失控。完整的路径枚举和消除是否可能是更好的选择，即更快、更健壮？任何帮助、想法和指示都是有帮助的。

Answer 1

这是你想要的吗？

find_path <- function(nperiods, opts = LETTERS[seq_len(period)]) {
  stopifnot(length(opts) == nperiods)
  out <- matrix(nrow = 2 ^ (nperiods - 1L), ncol = nperiods)
  r <- 1L
  recur_ <- function(period, branch, outcome) {
    if (period > length(branch)) {
      out[r, ] <<- opts[branch]
      r <<- r + 1L
      return(NULL)
    }
    for (i in c(outcome, period)) {
      branch[[period]] <- i
      recur_(period + 1L, branch, i)
    }
  }
  recur_(1L, integer(nperiods), NULL)
  out
}

calc_prob <- function(mat) {
  ps <- dimnames(mat)[[1L]]; if (is.null(ps)) ps <- seq_len(nrow(mat))
  ops <- dimnames(mat)[[2L]]; if (is.null(ops)) ops <- seq_len(ncol(mat))
  paths <- find_path(nrow(mat), ops)
  out <- vapply(seq_len(ncol(paths))[-1L], function(i) {
    comp <- ops[[i]]
    comp <- ifelse(paths[, i] == comp, paths[, i - 1L], comp)
    x <- exp(mat[i, paths[, i]])
    y <- exp(mat[i, comp])
    x / (x + y)
  }, numeric(nrow(paths)))
  dimnames(out) <- NULL; out <- cbind(1, out)
  dimnames(out)[[2L]] <- dimnames(paths)[[2L]] <- ps
  list(paths = paths, probs = out)
}

输出

> calc_prob(x) # x is the same lower-triangular matrix as shown in your example.

$paths
     t=1 t=2 t=3 t=4
[1,] "A" "A" "A" "A"
[2,] "A" "A" "A" "D"
[3,] "A" "A" "C" "C"
[4,] "A" "A" "C" "D"
[5,] "A" "B" "B" "B"
[6,] "A" "B" "B" "D"
[7,] "A" "B" "C" "C"
[8,] "A" "B" "C" "D"

$probs
     t=1       t=2       t=3        t=4
[1,]   1 0.5404888 0.8242364 0.08261823
[2,]   1 0.5404888 0.8242364 0.91738177
[3,]   1 0.5404888 0.1757636 0.28985432
[4,]   1 0.5404888 0.1757636 0.71014568
[5,]   1 0.4595112 0.8044942 0.36037495
[6,]   1 0.4595112 0.8044942 0.63962505
[7,]   1 0.4595112 0.1955058 0.28985432
[8,]   1 0.4595112 0.1955058 0.71014568

变量paths 为您提供每个时期t 的所有可能结果； probs 告诉您相应结果的概率。 但是，请注意，这样的概率树会随着周期数的增加而呈指数增长。等式是

其中N 是周期t 内所有可能路径的数量。仅 20 个周期，您将拥有 524288 条不同的路径。如果周期数达到 30，您将有 536870912 条不同的路径，而 R 无法处理该数量的计算。我建议您重新考虑您的预期输出。您是否正在运行具有一些其他约束的模拟，而不仅仅是时间依赖性，以便我们可以进一步修剪一些不必要的路径？或者您可能只需要一些汇总统计信息，例如预期值，这样我们就不必生成所有可能的路径？一定有比仅仅使用这种蛮力方法更好的方法。