Ручное моделирование цепи Маркова в R - PullRequest
Новая
       86

Ручное моделирование цепи Маркова в R

0 голосов
/ 24 апреля 2019

Рассмотрим цепочку Маркова с пространством состояний S = {1, 2} , матрица переходов

enter image description here

и начальное распределение α = (1/2, 1/2) .

  1. Имитировать 5 ступеней цепи Маркова (то есть имитировать X 0 , X 1 , ..., X 5 ).Повторите симуляцию 100 раз.Используйте результаты моделирования для решения следующих задач.

    • Оценка P (X 1 = 1 | X 0 = 1) .Сравните ваш результат с точной вероятностью.

Мое решение: 

# returns Xn 
func2 <- function(alpha1, mat1, n1) 
{
  xn <- alpha1 %*% matrixpower(mat1, n1+1)

  return (xn)
}

alpha <- c(0.5, 0.5)
mat <- matrix(c(0.5, 0.5, 0, 1), nrow=2, ncol=2)
n <- 10


for (variable in 1:100) 
{
   print(func2(alpha, mat, n))
}

Какая разница, если я запускаю этот кододин или 100 раз (как сказано в постановке задачи)?

Как с этого момента найти условную вероятность?

1 Ответ

3 голосов
/ 24 апреля 2019

Пусть 

alpha <- c(1, 1) / 2
mat <- matrix(c(1 / 2, 0, 1 / 2, 1), nrow = 2, ncol = 2) # Different than yours

- начальное распределение и матрица перехода. Ваш func2 находит только n-й шаг распределения, который не нужен и ничего не моделирует. Вместо этого мы можем использовать 

chainSim <- function(alpha, mat, n) {
  out <- numeric(n)
  out[1] <- sample(1:2, 1, prob = alpha)
  for(i in 2:n)
    out[i] <- sample(1:2, 1, prob = mat[out[i - 1], ])
  out
}

где out[1] генерируется с использованием только начального распределения, а затем для последующих членов мы используем матрицу перехода.

Тогда у нас есть 

set.seed(1)
# Doing once
chainSim(alpha, mat, 1 + 5)
# [1] 2 2 2 2 2 2

так что цепь инициировалась в 2 и застряла там из-за указанных вероятностей перехода.

Делая это 100 раз, мы имеем 

# Doing 100 times
sim <- replicate(chainSim(alpha, mat, 1 + 5), n = 100)
rowMeans(sim - 1)
# [1] 0.52 0.78 0.87 0.94 0.99 1.00

где последняя строка показывает, как часто мы оказывались в состоянии 2, а не в 1. Это дает одну (из многих) причин, почему 100 повторений более информативны: мы застряли в состоянии 2, выполняя только одну симуляцию, повторяя при этом 100 раз мы исследовали больше возможных путей.

Тогда условная вероятность может быть найдена с помощью 

mean(sim[2, sim[1, ] == 1] == 1)
# [1] 0.4583333

в то время как истинная вероятность равна 0,5 (определяется верхней левой записью матрицы перехода).

...