如何找到马氏链概率？

我试图找到链状态从k-1跳转到状态1，然后命中状态k的概率。任何人都可以发现我的错误吗？如何找到马氏链概率？

我试图模拟马尔科夫链，但我想做一个代码，让我找到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 
}

来源

2017-10-12 Ahmed Jyad

我觉得行

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 
}

来源

2017-10-12 17:29:02

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

其他注意事项：在两个不同的嵌套级别使用相同的索引是非常糟糕的形式。另外，当你开始有大量的if-then-else语句列表时，是时候重新考虑你的代码了。在这种情况下，您可以在p1 [i]> = Y上仅1:17进行子集化，取最小值，然后将X设置为该值。

来源

2017-10-12 17:46:43 Acccumulation

如何找到马氏链概率？

回答

相关问题