2017-10-12 9 views
0

チェーンが状態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 
} 

答えて

0

私はライン

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 
} 
0

あなたはXをk-1に設定しています。 Rでは長さ1のベクトルとして扱われます.Xが2になると、Xには2番目の要素がないため、X [i]はインデックスエラーを返します。

注:2つの異なる入れ子レベルで同じインデックスを使用することは非常に悪い形式です。また、if-then-elseステートメントの大量のリストが始まると、コードを再考するときです。この場合、p1 [i]> = Yに1:17をサブセット化し、最小値をとり、Xをそれに設定することができます。

関連する問題