用于离散马尔可夫链模拟的 R 库

发布于 2024-08-30 22:05:48 字数 168 浏览 6 评论 0原文

我正在寻找类似“msm”包的东西，但适用于离散马尔可夫链。例如，如果我有一个

Pi <- matrix(c(1/3,1/3,1/3,
0,2/3,1/6,
2/3,0,1/2))

针对状态 A、B、C 定义的转移矩阵。如何根据该转移矩阵模拟马尔可夫链？

原文

I am looking for something like the 'msm' package, but for discrete Markov chains. For example, if I had a transition matrix defined as such

Pi <- matrix(c(1/3,1/3,1/3,
0,2/3,1/6,
2/3,0,1/2))

for states A,B,C. How can I simulate a Markov chain according to that transition matrix?

分享到QQ

分享到微博

如果你对这篇内容有疑问，欢迎到本站社区发帖提问参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

发布评论

需要登录才能够评论，你可以免费注册一个本站的账号。

锦爱 2024-09-06 22:05:48

不久前，我编写了一组用于模拟和估计离散马尔可夫链概率矩阵的函数：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

A while back I wrote a set of functions for simulation and estimation of Discrete Markov Chain probability matrices: http://www.feferraz.net/files/lista/DTMC.R.

Relevant code for what you're asking:

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

回复收藏 0 原文

提赋 2024-09-06 22:05:48

啊，当我为你写下它时，你找到了解决方案。这是我想出的一个简单的例子：

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，而它应该是这样。我的例子有一个稍微改变的概率转移矩阵，它遵循这个规则。

Argh, you found the solution whilst I was writing it up for you. Here's a simple example that I came up with:

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();

Note that your probability transition matrix doesn't sum to 1 in each row, which it should do. My example has a slightly altered probability transition matrix which adheres to this rule.

回复收藏 0 原文