我正在尋找'msm'包,但是用於離散馬爾可夫鏈。例如,如果我有定義爲這樣用於離散馬爾可夫鏈模擬的R庫
Pi <- matrix(c(1/3,1/3,1/3,
0,2/3,1/6,
2/3,0,1/2))
對於狀態A,B,C的過渡矩陣。如何根據該轉換矩陣模擬馬爾可夫鏈?
感謝,
我正在尋找'msm'包,但是用於離散馬爾可夫鏈。例如,如果我有定義爲這樣用於離散馬爾可夫鏈模擬的R庫
Pi <- matrix(c(1/3,1/3,1/3,
0,2/3,1/6,
2/3,0,1/2))
對於狀態A,B,C的過渡矩陣。如何根據該轉換矩陣模擬馬爾可夫鏈?
感謝,
後來我寫了一組用於模擬和估計離散馬爾可夫鏈概率矩陣的函數:http://www.feferraz.net/files/lista/DTMC.R。你在問什麼
相關代碼:
simula <- function(trans,N) {
transita <- function(char,trans) {
sample(colnames(trans),1,prob=trans[char,])
}
sim <- character(N)
sim[1] <- sample(colnames(trans),1)
for (i in 2:N) {
sim[i] <- transita(sim[i-1],trans)
}
sim
}
#example
#Obs: works for N >= 2 only. For higher order matrices just define an
#appropriate mattrans
mattrans <- matrix(c(0.97,0.03,0.01,0.99),ncol=2,byrow=TRUE)
colnames(mattrans) <- c('0','1')
row.names(mattrans) <- c('0','1')
instancia <- simula(mattrans,255) # simulates 255 steps in the process
哎呀,你找到了解決辦法,而我在寫它給你的。這裏有一個簡單的例子,我想出了:
run = function()
{
# The probability transition matrix
trans = matrix(c(1/3,1/3,1/3,
0,2/3,1/3,
2/3,0,1/3), ncol=3, byrow=TRUE);
# The state that we're starting in
state = ceiling(3 * runif(1, 0, 1));
cat("Starting state:", state, "\n");
# Make twenty steps through the markov chain
for (i in 1:20)
{
p = 0;
u = runif(1, 0, 1);
cat("> Dist:", paste(round(c(trans[state,]), 2)), "\n");
cat("> Prob:", u, "\n");
newState = state;
for (j in 1:ncol(trans))
{
p = p + trans[state, j];
if (p >= u)
{
newState = j;
break;
}
}
cat("*", state, "->", newState, "\n");
state = newState;
}
}
run();
請注意,您的概率轉移矩陣不每一行,它應該做的總和爲1。我的例子有一個稍微改變的概率轉換矩陣,它遵守這個規則。
感謝答案。你的代碼非常易讀。對此,我真的非常感激。 – stevejb 2010-05-02 20:56:18
可讀的代碼?根據我的經驗,這個概念在大多數使用'R'的人中完全沒有了。希望能幫助到你! – icio 2010-05-02 22:09:56
爲了生成一個從1到3的隨機整數,我認爲'sample(1:3,1)'比'ceiling(3 * runif(1,0,1))'更容易理解。另外,對於最內層的for循環,您可以簡單地使用'newState < - sample(1:ncol(trans),1,prob = trans [state,])''。這更清楚地表明發生了什麼事情。然後你甚至不需要對行進行標準化。 – 2010-05-03 14:51:51
您現在可以使用在CRAN可用markovchain
包。用戶manual。很不錯,有幾個例子。
NVM:找到了答案 http://books.google.com/books?id=AALhk_mt7SYC&lpg=PA116&dq=r%20markov%20chain&pg=PA119#v=onepage&q=r%20markov%20chain&f=false – stevejb 2010-05-02 19:04:54