Блог

Чтобы проиллюстрировать проблему, давайте определим следующую матрицу (где 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 мы вычисляем 2 вероятности из 1 группировки, в t=3 мы вычисляем 4 вероятности из 2 группировок (обратите внимание, как один вариант исключается из сравнения из-за последовательной зависимости и неотъемлемого предположения, что он не был самым высоким в t-1 ) и в t=4 , мы вычисляем 8 вероятностей из 4 группировок. Затем конечные вероятности являются произведением вероятностей в каждом 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)]))
  

который выглядит следующим образом (1 из 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
  

И мы можем вычислить вероятности в каждом периоде следующим образом (все строки суммируются до единицы, как и должно быть):

 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 роста? Или хороший способ заполнения набора маскирующих матриц для перебора? Я пытаюсь избежать проблемы, выходящей из-под контроля. Возможно ли, что полное перечисление и исключение путей является лучшим вариантом, т. Е. Более быстрым и надежным? Любая помощь, идеи и указатели полезны.

Это то, что вы хотите?

 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 просто не сможет обработать такое количество вычислений. Я предлагаю вам пересмотреть свои ожидаемые результаты. Выполняете ли вы симуляцию с некоторыми другими ограничениями, кроме зависимости от времени, чтобы мы могли дополнительно обрезать некоторые ненужные пути? Или, может быть, вам нужна только некоторая сводная статистика, такая как ожидаемое значение, чтобы нам не приходилось генерировать все возможные пути? Должен быть лучший способ, чем просто использовать подход грубой силы, подобный этому.