如何找到马尔可夫链概率？答案

【问题标题】：How to find the Markov Chain Probability?如何找到马尔可夫链概率？
【发布时间】：2018-03-24 17:44:28
【问题描述】：

我试图找出链在到达状态 k 之前从状态 k-1 跳转到状态 1 的概率。谁能发现我的错误？

我试图模拟马尔可夫链，但我想编写一个代码，让我找到k ={1, 2, 3, ........17} 的概率。但是我真的拿不到代码。

这是我经常收到的错误消息

Error in while (X[i] > 1 && X[i] < k) { : 
  missing value where TRUE/FALSE needed

这是我的代码：

k <- 17
{   p <- 0.5
q <- 0.1
P <- matrix (0, nrow = k, ncol = k, byrow = TRUE)
for (i in 1:k)
{   for (j in 1:k)
    {   if (i == 1 && i == j)
        {   P[i,j] <- 1
        }
        else if (i == k && i == j)
        {   P[i,j] <- 1
        }
        else if (i == j)
        {   P[i,j] <- p*(1-q)
        }
        else if (j == k && i != 1)
        {   P[i,j] <- q
        }   
        else if (i == j+1 && i != k)
        {   P[i,j] <- (1-p)*(1-q)
        }
    }
}
P
X <- (k-1) 
trials <- 1000
hits <- 0 #counter for no. of hits 
for (i in 1:trials)
{   i <- 1 #no. of steps
    while(X[i] > 1 && X[i] < k)
    {   Y <- runif(1) #uniform samples
        p1 <- P[X[i],] #calculating the p-value
        p1 <- cumsum(p1)
        # changes in the chain
        if(Y <= p1[1])
        {   X[i+1] = 1}
        else if(Y <= p1[2])
        {   X[i+1] = 2}
        else if(Y <= p1[3])
        {   X[i+1] = 3}
        else if(Y <= p1[4])
        {   X[i+1] = 4}
        else if(Y <= p1[5])
        {   X[i+1] = 5}
        else if(Y <= p1[6])
        {   X[i+1] = 6}
        else if(Y <= p1[7])
        {   X[i+1] = 7}
        else if(Y <= p1[8])
        {   X[i+1] = 8}
        else if(Y <= p1[9])
        {   X[i+1] = 9}
        else if(Y <= p1[10])
        {   X[i+1] = 10}
        else if(Y <= p1[11])
        {   X[i+1] = 11}
        else if(Y <= p1[12])
        {   X[i+1] = 12}
        else if(Y <= p1[13])
        {   X[i+1] = 13}
        else if(Y <= p1[14])
        {   X[i+1] = 14}
        else if(Y <= p1[15])
        {   X[i+1] = 15}
        else if(Y <= p1[16])
        {   X[i+1] = 16}
        else if(Y <= p1[17])
        {   X[i+1] <= 17}
        i <- i+1
    }
    if(X[i]==1)
    {   hits <- hits+1}
    else
    {   hits <- hits+0}
}

Probability <- hits/trials
Probability
}

【问题讨论】：

标签： r simulation probability markov-chains

【解决方案1】：

我觉得行

i <- 1 #no. of steps

不应该在那里。试试这个：

k <- 17
{   p <- 0.5
q <- 0.1
P <- matrix (0, nrow = k, ncol = k, byrow = TRUE)
for (i in 1:k)
{   for (j in 1:k)
    {   if (i == 1 && i == j)
        {   P[i,j] <- 1
        }
        else if (i == k && i == j)
        {   P[i,j] <- 1
        }
        else if (i == j)
        {   P[i,j] <- p*(1-q)
        }
        else if (j == k && i != 1)
        {   P[i,j] <- q
        }   
        else if (i == j+1 && i != k)
        {   P[i,j] <- (1-p)*(1-q)
        }
    }
}
P
X <- (k-1) 
trials <- 1000
hits <- 0 #counter for no. of hits 
for (i in 1:trials)
{
    while(X[i] > 1 && X[i] < k)
    {   Y <- runif(1) #uniform samples
        p1 <- P[X[i],] #calculating the p-value
        p1 <- cumsum(p1)
        # changes in the chain
        if(Y <= p1[1])
        {   X[i+1] = 1}
        else if(Y <= p1[2])
        {   X[i+1] = 2}
        else if(Y <= p1[3])
        {   X[i+1] = 3}
        else if(Y <= p1[4])
        {   X[i+1] = 4}
        else if(Y <= p1[5])
        {   X[i+1] = 5}
        else if(Y <= p1[6])
        {   X[i+1] = 6}
        else if(Y <= p1[7])
        {   X[i+1] = 7}
        else if(Y <= p1[8])
        {   X[i+1] = 8}
        else if(Y <= p1[9])
        {   X[i+1] = 9}
        else if(Y <= p1[10])
        {   X[i+1] = 10}
        else if(Y <= p1[11])
        {   X[i+1] = 11}
        else if(Y <= p1[12])
        {   X[i+1] = 12}
        else if(Y <= p1[13])
        {   X[i+1] = 13}
        else if(Y <= p1[14])
        {   X[i+1] = 14}
        else if(Y <= p1[15])
        {   X[i+1] = 15}
        else if(Y <= p1[16])
        {   X[i+1] = 16}
        else if(Y <= p1[17])
        {   X[i+1] <= 17}
        i <- i+1
    }
    if(X[i]==1)
    {   hits <- hits+1}
    else
    {   hits <- hits+0}
}

Probability <- hits/trials
Probability
}

【讨论】：

【解决方案2】：

您将 X 设置为 k-1。在 R 中，它被视为长度为 1 的向量。一旦 i 达到 2，X[i] 返回索引错误，因为 X 没有第二个元素。

进一步说明：在两个不同的嵌套级别中使用相同的索引是非常糟糕的形式。此外，当您开始拥有大量 if-then-else 语句时，是时候重新考虑您的代码了。在这种情况下，您可以只在 p1[i] >=Y 上设置 1:17 的子集，取最小值，然后将 X 设置为那个值。

【讨论】：