在 R 中使用 markovchainFit 获得马尔可夫链后如何识别状态转移概率？

Question

我有一系列事件：

library(markovchain)
sequence<-c("LHR - BA","BOS - BA","BOS - ZE","IAD - ZE","BOS - BA","LHR - BA",
            "LGW - BA","TPA - BA","TPA - BA","LGW - BA","LHR - BA","BOS - BA",
            "BOS - ZE","BOS - ZE","BOS - BA","LHR - BA","LHR - BA","BOS - BA",
            "BOS - ZE","BOS - ZE","BOS - DL","ATL - DL","LHR - BA","BRU - BA")

使用这个序列，我得到一个使用以下函数的马尔可夫链：

sequenceMatr <- createSequenceMatrix(sequence, sanitize=FALSE)
mcFit        <- markovchainFit(data=sequence, method="mle")

考虑下一个状态是"LHR - BA"；那么如何以以下格式识别跨州的概率分布：

           "LHR - BA", "BOS - ZE", "IAD - ZE", "BOS - BA", "LHR - BA", "LGW - BA", "TPA - BA"
"LHR - BA"     0.1   ,     0.2   ,    0.2    ,     0.3   ,     0.1   ,     0.1   ,     0.1

Answer 1

我觉得你的问题有点难以理解，但这是我的解释。

考虑下一个状态是"LHR - BA"；那么我如何确定跨州的概率分布

所以按照我的阅读方式，您知道在某个时间 t+1 您的系统处于状态 "LHR - BA" 并且您想知道时间 t 的分布概率我>。换句话说，你想要条件概率

P(S(t)=x | S(t+1)="LHR - BA")

根据Bayes' law这个概率等于

P(S(t)=x) * P(S(t+1)="LHR - BA" | S(t)=x) / P(S(t+1)="LHR - BA")

为此，您需要对无条件分布进行一些估计。对于大 t 和合理的（不知道这里的正确术语）马尔可夫链，您可以在这里简单地采用（希望是唯一的）稳定状态。通过这种解释，您可以以非常直接的方式将上述公式转换为 R：

mcEst <- mcFit$estimate
mcSteady <- steadyStates(mcEst)
sapply(states(mcEst), function(x) transitionProbability(mcEst, x, "LHR - BA")*mcSteady[1,x]/mcSteady[1,"LHR - BA"])

但也许您想以更好的方式编写此代码，并为每个可能的S(t+1) 提供行，而不仅仅是"LHR - BA"。为此，您必须更直接地处理转换矩阵，而不是调用 transitionProbability，因为该方法不适用于矢量化参数。

tm <- mcEst@transitionMatrix
if (mcEst@byrow) tm <- t(tm)
res <- tm * t(outer(mcSteady[1,], mcSteady[1,], "/"))

前两行获取转移矩阵并确保行（第一个索引）是目标状态，列（第二个索引）是源状态。请参阅getMethods("transitionProbability") 了解执行此类操作的已发布功能。所以那时你有

tm[i,j] = P(S(t+1)=i | S(t)=j)

然后你采取稳定状态，并做所有可能的组合。你得到

outer(mcSteady[1,], mcSteady[1,], "/")[i,j] = mcSteady[1,i]/mcSteady[1,j]

这是错误的方式。所以你转置它并得到

t(outer(mcSteady[1,], mcSteady[1,], "/"))[i,j] = mcSteady[1,j]/mcSteady[1,i]

你乘以tm得到最终结果：

res[i,j] = P(S(t+1)=i | S(t)=j) * P(S(t)=j) / P(S(t+1)=i)

结果表的每一行将是给定后继状态的一个分布。包括你要求的那个：

> res["LHR - BA",]
  ATL - DL   BOS - BA   BOS - DL   BOS - ZE   BRU - BA   IAD - ZE   LGW - BA 
0.23379630 0.37037037 0.00000000 0.00000000 0.02083333 0.00000000 0.20833333 
  LHR - BA   TPA - BA 
0.16666667 0.00000000