Нахождение стационарного распределения Марковского процесса при заданной матрице вероятностей перехода

Существует два потока, связанных с этой проблемой переполнения стека:

Все вышесказанное просто, но очень дорого. Если у нас есть матрица перехода порядка n, то на каждой итерации мы вычисляем матрицу-матрицу умножения по затратам O(n ^ 3).

Существует ли более эффективный способ сделать это? Одна вещь, которая приходит мне в голову, - это использование собственного разложения. Известно, что матрица Маркова:

  • быть диагонализуемым в комплексной области: A = E * D * E^{-1};
  • имеют действительное собственное значение 1, а другие (комплексные) собственные значения с длиной меньше 1.
Стационарное распределение - это собственный вектор, связанный с собственным значением 1, то есть первый собственный вектор. Что ж, теория хороша, но я не могу заставить ее работать. Принимая матрицу P в первом связанном вопросе: 
P <- structure(c(0, 0.1, 0, 0, 0, 0, 0, 0.1, 0.2, 0, 0, 0, 0, 0, 0.2, 
0.3, 0, 0, 0.5, 0.4, 0.3, 0.5, 0.4, 0, 0, 0, 0, 0, 0.6, 0.4, 
0.5, 0.4, 0.3, 0.2, 0, 0.6), .Dim = c(6L, 6L))

Если я это сделаю: 

Re(eigen(P)$vectors[, 1])
# [1] 0.4082483 0.4082483 0.4082483 0.4082483 0.4082483 0.4082483

Что происходит? Согласно предыдущим вопросам, стационарное распределение имеет вид: 

# [1] 0.002590673 0.025906737 0.116580322 0.310880848 0.272020713 0.272020708
2 5
matrix statistics markov-chains probability

2 ответа:

Ваш вектор y = Re(eigen(P)$vectors[, 1]) не является распределением (поскольку он не складывается в единицу) и решает P'y = y, а не x'P = x. Решение из вашего связанного Q&A приблизительно решает последнее:

x = c(0.00259067357512953, 0.0259067357512953, 0.116580310880829, 
0.310880829015544, 0.272020725388601, 0.272020725388601)
all(abs(x %*% P - x) < 1e-10) # TRUE

Транспонируя P, вы можете использовать свой подход к собственным значениям: 

x2 = Re(eigen(t(P))$vectors[, 1])
x2 <- x2/sum(x2) 
(function(x) all(abs(x %*% P - x) < 1e-10))(
  x2
) # TRUE
Однако в данном случае он находит другой стационарный вектор.
2

Ну, чтобы использовать собственное разложение, нам нужно работать с t(P).

Определение матрицы вероятности перехода отличается между вероятностью / статистикой и линейной алгеброй. В статистике все строки P суммируются в 1, в то время как в линейной алгебре все столбцы P суммируются в 1. Поэтому вместо eigen(P) нам нужно eigen(t(P)): 

e <- Re(eigen(t(P))$vectors[, 1])
e / sum(e)
# [1] 0.002590673 0.025906737 0.116580322 0.310880848 0.272020713 0.272020708
Как мы видим, мы использовали только первый собственный вектор, то есть собственный вектор наибольшего собственного значения. Поэтому нет необходимости вычислять все Собственные значения / векторы с использованием eigen. Степенной метод может быть использован для нахождения собственного вектора наибольшего собственного значения. Реализуем это в R: 
stydis1 <- function (A) {
  n <- dim(A)[1L]
  ## checking
  if (any(.rowSums(A, n, n) != 1)) 
    stop (" 'A' is not a Markov matrix")
  ## implement power method
  e <- runif (n)
  oldnorm <- sqrt(c(crossprod(e)))
  repeat {
    e <- crossprod(A, e)
    newnorm <- sqrt(c(crossprod(e)))
    if (abs(newnorm / oldnorm - 1) < 1e-8) break
    e <- e / newnorm
    oldnorm <- newnorm
    }
  ## rescale `e` so that it sums up to 1
  c(e / sum(e))
  }

stydis1 (P)
# [1] 0.002590673 0.025906737 0.116580322 0.310880848 0.272020713 0.272020708

И результат правильный.

На самом деле, нам не нужно использовать собственное разложение. Мы можем скорректировать метод, используемый в вашем втором связанном вопросе. Там мы взяли матричную мощность, которая, как вы заметили, стоит дорого; но почему бы не перелить ее в матрицу-векторное умножение?

stydis2 <- function (A) {
  n <- dim(A)[1L]
  ## checking
  if (any(.rowSums(A, n, n) != 1)) 
    stop (" 'A' is not a Markov matrix")
  ## direct computation
  b <- A[1, ]
  oldnorm <- sqrt(c(crossprod(b)))
  repeat {
    b <- crossprod(A, b)
    newnorm <- sqrt(c(crossprod(b)))
    if (abs(newnorm / oldnorm - 1) < 1e-8) break
    oldnorm <- newnorm
    }
  ## return stationary distribution
  c(b)
  }

stydis2 (P)
# [1] 0.002590673 0.025906737 0.116580322 0.310880848 0.272020713 0.272020708

Мы начинаем из произвольного начального распределения, скажем A[1, ], и итеративно применять матрицу перехода, пока распределение не сойдет. Опять же, результат правильный.

5