В чем разница между matrixpower () и markov (), когда дело доходит до вычисления P ^ n?

Question

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

P =  0.1  0.2  0.4  0.3
     0.4  0.0  0.4  0.2
     0.3  0.3  0.0  0.4
     0.2  0.1  0.4  0.3 

И взгляните на следующий исходный код:

# markov function
markov <- function(init,mat,n,labels) 
{ 
    if (missing(labels)) 
    {
        labels <- 1:length(init) 
    }

    simlist <- numeric(n+1) 
    states <- 1:length(init) 
    simlist[1] <- sample(states,1,prob=init)  

    for (i in 2:(n+1))
    { 
        simlist[i] <- sample(states, 1, prob = mat[simlist[i-1],])  
    }

    labels[simlist] 
}

# matrixpower function
matrixpower <- function(mat,k) 
{
    if (k == 0) return (diag(dim(mat)[1])) 
    if (k == 1) return(mat)
    if (k > 1) return( mat %*% matrixpower(mat, k-1))
}

tmat = matrix(c(0.1, 0.2, 0.4, 0.3,
                0.4, 0.0, 0.4, 0.2,
                0.3, 0.3, 0.0, 0.4,
                0.2, 0.1, 0.4, 0.3), nrow=4, ncol=4, byrow=TRUE)

p10 = matrixpower(mat = tmat, k=10)  
rowMeans(p10)

nn <- 10 
alpha <- c(0.25, 0.25, 0.25, 0.25)

set.seed(1)

steps <- markov(init=alpha, mat=tmat, n=nn)
table(steps)/(nn + 1)

Выход

> rowMeans(p10)
[1] 0.25 0.25 0.25 0.25
> 
.
.
.
> table(steps)/(nn + 1)
steps
         1          2          3          4 
0.09090909 0.18181818 0.18181818 0.54545455 
> ?rowMeans

Почему результаты так отличаются?

В чем разница между использованием matrixpower() и markov(), когда дело доходит до вычисления P ^п * * 1023

Answer 1

В настоящее время вы сравниваете совершенно разные вещи.Во-первых, я сосредоточусь не на вычислениях P ⁿ , а на A * P ⁿ , где A - начальное распределение.В этом случае matrixpower выполняет работу:

p10 <- matrixpower(mat = tmat, k = 10)  
alpha <- c(0.25, 0.25, 0.25, 0.25)
alpha %*% p10
#           [,1]      [,2]      [,3]      [,4]
# [1,] 0.2376945 0.1644685 0.2857105 0.3121265

- это истинные вероятности нахождения в состояниях 1, 2, 3, 4, соответственно, после 10 шагов (после первоначального розыгрыша с использованием A).

Между тем, markov(init = alpha, mat = tmat, n = nn) является только одной реализацией длины nn + 1, и только A * P ⁿ имеет отношение только к последнему числу этой реализации.Итак, чтобы попытаться получить числа, подобные теоретическим, нам нужно много реализаций с nn <- 10, как в

table(replicate(markov(init = alpha, mat = tmat, n = nn)[nn + 1], n = 10000)) / 10000
#
#      1      2      3      4 
# 0.2346 0.1663 0.2814 0.3177

, где я моделирую 10000 реализаций и принимаю только последнее состояние каждой реализации.