1

私は一連のイベントを持っています:

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
4

1 に答える 1

1

あなたの質問は少しわかりにくいと思いますが、これが私の解釈です。

次の状態を考えてみましょう"LHR - BA"。州全体の確率分布を特定するにはどうすればよいですか

したがって、これを読んだ方法では、ある時点t +1 でシステムが状態にあり、時点t"LHR - BA"での分布確率を知りたいことがわかります。つまり、条件付き確率が必要です

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

ベイズの法則によれば、この確率は次のようになります。

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,], "/"))

最初の 2 行は遷移行列を取得し、行 (最初のインデックス) がターゲット状態で、列 (2 番目のインデックス) がソース状態であることを確認します。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)

その結果のテーブルの各行は、特定の後継状態の 1 つの分布になります。あなたが求めたものを含む：

> 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
于 2014-01-22T10:20:05.157 に答える