Szacowanie prawdopodobieństwa przejścia Markowa na podstawie danych sekwencji

16

Mam pełny zestaw sekwencji (dokładnie 432 obserwacje) z 4 stanów AD : np

Y=(ACDDBACBAACABCADABA)

EDYCJA : Sekwencje obserwacji mają nierówne długości! Czy to coś zmienia?

Czy istnieje sposób obliczania macierzy przejścia

Pij(Yt=j|Yt1=i)
w Matlabie lub R lub podobnym? Myślę, że pakiet HMM może pomóc. jakieś pomysły?

np .: Szacowanie prawdopodobieństwa łańcucha Markowa

HCAI
źródło
3
Masz 4 stany: S={1:=A,2:=B,3:=C,4:=D} . Niech nij będzie liczbą razy, gdy łańcuch dokonał przejścia ze stanu i do stanu j , dla ij,=1,2,3,4 . Oblicz nij„S z próbki i oszacować macierz transformacji (pij) przez maksimum prawdopodobieństwo użyciu szacunkowych P i J = n i j / Σ 4 J = 1 n i j . p^ij=nij/j=14nij
Zen.
Te notatki pochodzą z szacunków MLE: stat.cmu.edu/~cshalizi/462/lectures/06/markov-mle.pdf
Zen.
2
Podobne pytanie: stats.stackexchange.com/questions/26722/…
B_Miner
@B_Miner czy mógłbyś dla mnie napisać kod w formie pseudokodu? Lub wyjaśnij to w kategoriach laikatu ... Widzę jednak, że działa na mojej konsoli R.
HCAI,
Mam pytanie: rozumiem twoją implementację i jest ona dla mnie w porządku, ale zastanawiałem się, dlaczego nie mogę po prostu użyć funkcji hmmestimate Matlaba do obliczenia macierzy T? Coś w stylu: stany = [1,2,3,4] [T, E] = hmmestimate (x, stany); gdzie T jest interesującą mnie macierzą przejścia. Jestem nowy w łańcuchach Markowa i HMM, więc chciałbym zrozumieć różnicę między tymi dwoma implementacjami (jeśli takie istnieją).
Dowolny

Odpowiedzi:

18

Sprawdź komentarze powyżej. Oto szybkie wdrożenie w R.

x <- c(1,2,1,1,3,4,4,1,2,4,1,4,3,4,4,4,3,1,3,2,3,3,3,4,2,2,3)
p <- matrix(nrow = 4, ncol = 4, 0)
for (t in 1:(length(x) - 1)) p[x[t], x[t + 1]] <- p[x[t], x[t + 1]] + 1
for (i in 1:4) p[i, ] <- p[i, ] / sum(p[i, ])

Wyniki:

> p
          [,1]      [,2]      [,3]      [,4]
[1,] 0.1666667 0.3333333 0.3333333 0.1666667
[2,] 0.2000000 0.2000000 0.4000000 0.2000000
[3,] 0.1428571 0.1428571 0.2857143 0.4285714
[4,] 0.2500000 0.1250000 0.2500000 0.3750000

(Prawdopodobnie głupia) implementacja w MATLAB (z której nigdy nie korzystałem, więc nie wiem, czy to zadziała. Właśnie przejrzałem „zadeklaruj macierz wektorową MATLAB”, aby uzyskać składnię):

x = [ 1, 2, 1, 1, 3, 4, 4, 1, 2, 4, 1, 4, 3, 4, 4, 4, 3, 1, 3, 2, 3, 3, 3, 4, 2, 2, 3 ]
n = length(x) - 1
p = zeros(4,4)
for t = 1:n
  p(x(t), x(t + 1)) = p(x(t), x(t + 1)) + 1
end
for i = 1:4
  p(i, :) = p(i, :) / sum(p(i, :))
end
Zen
źródło
Wygląda świetnie! Nie jestem jednak pewien, co robi trzecia linia w twoim kodzie (głównie dlatego, że znam Matlaba). Czy jest szansa, że ​​napiszesz go w matlabie lub pseudokodzie? Byłbym bardzo zobowiązany.
HCAI
2
Trzeci wiersz robi to: wartości łańcucha to . Dla t = 1 , , n - 1 , przyrost p x t , x t + 1 . x1,,xnt=1,,n1pxt,xt+1
Zen.
Czwarta linia normalizuje każdą linię macierzy . (pij)
Zen.
Naga tutaj moja powolność. Doceniam tłumaczenie kodu MATLAB, chociaż nadal nie widzę, co próbuje zrobić w pierwszej forpętli. Trzecia linia oryginalnego kodu liczy, ile razy przechodzi ze stanu x i do stanu x j ? Jeśli mógłbyś to powiedzieć słowami, bardzo bym to docenił. Pozdrawiamxxixj
HCAI
1
Nie, to tylko jeden wiersz. Nie łącz konkatenacji, ponieważ wprowadzisz „fałszywe” przejścia: ostatni stan jednej linii pierwszy stan następnej linii. Musisz zmienić kod, aby przejść przez linie macierzy i policzyć przejścia. Na koniec znormalizuj każdą linię macierzy przejścia. x
Zen.
9

Oto moja implementacja w języku R.

x <- c(1,2,1,1,3,4,4,1,2,4,1,4,3,4,4,4,3,1,3,2,3,3,3,4,2,2,3)
xChar<-as.character(x)
library(markovchain)
mcX<-markovchainFit(xChar)$estimate
mcX
Giorgio Spedicato
źródło
1
user32041 żądanie „s (wysłane jako edit zamiast komentarza, ponieważ on / ona pozbawiona reputacji): W jaki sposób można zmusić transitionMatrix rezultatu markovchainFit do data.frame?
chl
Możesz przekonwertować na za pomocą s ( m c X , d a t a . f r a m e )data.frameas(mcX,"data.frame")
Giorgio Spedicato,
@GiorgioSpedicato, czy możesz skomentować, jak radzić sobie z sekwencjami o nierównej długości (nie mogę połączyć), proszę w swoim pakiecie?
HCAI,
@HCAI, proszę zobaczyć aktualną stronę winiety 35-36
Giorgio Spedicato,
@GiorgioSpedicato dziękuję za odniesienie cran.r-project.org/web/packages/markovchain/vignettes/… . Nadal mam n macierzy przejściowych, po jednej dla każdej sekwencji. Poszukuję jednego ogólnego, który bierze pod uwagę wszystkie obserwacje sekwencji. Czy czegoś brakuje?
HCAI,
2

Oto sposób, aby to zrobić w Matlabie:

x = [1,2,1,1,3,4,4,1,2,4,1,4,3,4,4,4,3,1,3,2,3,3,3,4,2,2,3];
counts_mat = full(sparse(x(1:end-1),x(2:end),1));
trans_mat = bsxfun(@rdivide,counts_mat,sum(counts_mat,2))

Acknowledgement owed to SomptingGuy: http://www.eng-tips.com/viewthread.cfm?qid=236532

John
źródło