Estimando las probabilidades de transición de Markov a partir de datos de secuencia

16

Tengo un conjunto completo de secuencias (432 observaciones para ser precisos) de 4 estados : por ejemploAD

Y=(ACDDBACBAACABCADABA)

EDITAR : ¡Las secuencias de observación son de longitudes desiguales! ¿Esto cambia algo?

¿Hay alguna forma de calcular la matriz de transición en Matlab o R o similar? Creo que el paquete HMM podría ayudar. ¿Alguna idea?

Pij(Yt=j|Yt1=i)

por ejemplo: Estimando las probabilidades de la cadena de Markov

HCAI
fuente
3
Tiene estados: . Sea el número de veces que la cadena hizo una transición del estado al estado , para . Calcule los de su muestra y calcule la matriz de transición mediante la máxima probabilidad utilizando las estimaciones . 4S={1:=A,2:=B,3:=C,4:=D}nijijij,=1,2,3,4nij(pij)p^ij=nij/j=14nij
Zen
Estas notas derivan las estimaciones de MLE: stat.cmu.edu/~cshalizi/462/lectures/06/markov-mle.pdf
Zen
2
Pregunta similar: stats.stackexchange.com/questions/26722/…
B_Miner
@B_Miner, ¿podrías escribir tu código en forma de pseudocódigo para mí? O explicarlo en términos simples ... Sin embargo, veo que funciona en mi consola R.
HCAI
Tengo una pregunta: entiendo su implementación y me parece bien, pero me preguntaba ¿por qué no puedo simplemente usar la función himatestimate de Matlab para calcular la matriz T? Algo así como: estados = [1,2,3,4] [T, E] = himatestimate (x, estados); donde T es la matriz de transición en la que estoy interesado. Soy nuevo en las cadenas de Markov y HMM, así que me gustaría entender la diferencia entre las dos implementaciones (si hay alguna).
Cualquier

Respuestas:

18

Por favor, revise los comentarios anteriores. Aquí hay una implementación rápida en 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, ])

Resultados:

> 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

Una implementación (probablemente tonta) en MATLAB (que nunca he usado, así que no sé si esto va a funcionar. Acabo de buscar en Google "declarar matriz de vectores MATLAB" para obtener la sintaxis):

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
fuente
¡Se ve muy bien! Sin embargo, no estoy seguro de lo que hace la tercera línea en su código (principalmente porque estoy familiarizado con Matlab). ¿Alguna posibilidad de que puedas escribirlo en matlab o pseudocódigo? Estaría muy agradecido.
HCAI
2
La tercera línea hace esto: los valores de la cadena son . Para t = 1 , , n - 1 , incremente p x t , x t + 1 . x1,,xnt=1,,n1pxt,xt+1
Zen
La cuarta línea normaliza cada línea de la matriz . (pij)
Zen
Desnudo con mi lentitud aquí. Aprecio la traducción del código MATLAB, aunque todavía no puedo ver lo que intenta hacer en su primer forbucle. La tercera línea del código original está contando el número de veces que va de un estado x i al estado x j ? Si pudieras decirlo en palabras, te lo agradecería mucho Saludosxxixj
HCAI
1
No, es solo una fila. No concatene porque introducirá transiciones "falsas": último estado de una línea primer estado de la línea siguiente. Debe cambiar el código para recorrer las líneas de su matriz y contar las transiciones. Al final, normalice cada línea de la matriz de transición. x
Zen
9

Aquí está mi implementación en 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
fuente
1
Solicitud del usuario 32041 (publicada como una edición en lugar de un comentario ya que carece de reputación): ¿Cómo puedo forzar la transición Matriz del resultado markovchainFit a un data.frame?
chl
data.frameas(mcX,"data.frame")
@GiorgioSpedicato, ¿puede comentar sobre cómo manejar secuencias de longitudes desiguales (no puedo concatenar) en su paquete?
HCAI
@HCAI, vea la viñeta actual en la página 35-36
Giorgio Spedicato
@GiorgioSpedicato gracias por la referencia cran.r-project.org/web/packages/markovchain/vignettes/… . Todavía tengo n matrices de transición, una para cada secuencia. Lo que busco es uno general que tenga en cuenta todas las observaciones de secuencia. ¿Se me escapa algo?
HCAI
2

Aquí hay una manera de hacerlo en Matlab:

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

Reconocimiento adeudado a SomptingGuy: http://www.eng-tips.com/viewthread.cfm?qid=236532

John
fuente