时间戳

首页 > 解决方案 > 在R中枚举顺序概率树中的路径子集

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

问题描述

为了说明问题，让我们定义以下矩阵（其中 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的假设中，只是比较而不是因为它被假设为更高。我们可以将问题构造成这样的树：AACBt=2

在此处输入图像描述

因此，t=1对于概率为 1，因为t=2我们从 1 个分组中计算 2 个概率，因此t=3我们从 2 个分组中计算 4 个概率（请注意，由于顺序依赖性和固有假设它不是最​​高的，因此如何从比较中消除一个选项t-1）在 中t=4，我们从 4 个分组中计算出 8 个概率。然后，最终概率是t组成 8 条路径的每个概率的乘积。在真正的问题中，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 个可能的比较中的 1 个）：

    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

随着增长，是否有一种聪明的方法可以针对所有可能的路径执行此操作t？或者填充一组掩码矩阵以循环的好方法？我试图避免这个问题失控。完整的路径枚举和消除是否可能是更好的选择，即更快、更健壮？任何帮助、想法和指示都是有帮助的。

标签： rprobability

解决方案

这是你想要的吗？

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 无法处理该数量的计算。我建议您重新考虑您的预期输出。您是否正在运行具有一些其他约束的模拟，而不仅仅是时间依赖性，以便我们可以进一步修剪一些不必要的路径？或者您可能只需要一些汇总统计信息，例如预期值，这样我们就不必生成所有可能的路径？必须有比仅使用这样的蛮力方法更好的方法。

推荐阅读