Estimation des probabilités de transition de Markov à partir des données de séquence

J'ai un ensemble complet de séquences (432 observations pour être précis) de 4 états $A-D$ : par ex.

Y = (\begin{array}{ccccccc} A & C & D & D & B & A & C \\ B & A & A & C & A & - & - \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ B & C & A & D & A & B & A \end{array})

$Y=\left(\begin{array}{c c c c c c c} A& C& D&D & B & A &C\\ B& A& A&C & A&- &-\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ B& C& A&D & A & B & A\\ \end{array}\right)$

EDIT : Les séquences d'observation sont de longueurs inégales! Est-ce que cela change quelque chose?

Existe-t-il un moyen de calculer la matrice de transition dans Matlab ou R ou similaire? Je pense que le package HMM pourrait aider. Des pensées?

P_{i j} (Y_{t} = j | Y_{t - 1} = i)

$P_{ij}(Y_{t}=j|Y_{t-1}=i)$

Exemple: Estimation des probabilités de la chaîne de Markov

r matlab markov-process HCAI
la source

Vous avez

états:

. Soit

le nombre de fois que la chaîne est passée de l'état

à l'état

, pour

. Calculer le

4

$4$

S = {1 := A, 2 := B, 3 := C, 4 := D}

$S=\{1:=A,2:=B,3:=C,4:=D\}$

n_{i j}

$n_{ij}$

i

$i$

j

$j$

i j, = 1, 2, 3, 4

$ij,=1,2,3,4$

n_{i j}

$n_{ij}$ « s à partir de l' échantillon et d' estimer la matrice de transition

par maximum de vraisemblance en utilisant les estimations

(p_{i j})

$(p_{ij})$

{\hat{p}}_{i j} = n_{i j} / \sum_{j = 1}^{4} n_{i j}

$\hat{p}_{ij}=n_{ij}/\sum_{j=1}^4 n_{ij}$

Zen

Ces notes dérivent des estimations MLE: stat.cmu.edu/~cshalizi/462/lectures/06/markov-mle.pdf

Zen

Question similaire: stats.stackexchange.com/questions/26722/…

B_Miner

@B_Miner pourriez-vous écrire votre code sous forme de pseudo-code pour moi? Ou expliquez-le en termes simples ... Cependant, je vois que cela fonctionne dans ma console R.

HCAI

J'ai une question: je comprends votre implémentation et ça me va très bien, mais je me demandais pourquoi ne puis-je pas simplement utiliser la fonction hmmestimate de Matlab pour calculer la matrice T? Quelque chose comme: états = [1,2,3,4] [T, E] = hmmestimate (x, états); où T est la matrice de transition qui m'intéresse. Je suis nouveau dans les chaînes de Markov et HMM, donc j'aimerais comprendre la différence entre les deux implémentations (s'il y en a).

N'importe quel

Réponses:

Veuillez vérifier les commentaires ci-dessus. Voici une implémentation rapide dans 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, ])

Résultats:

> 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

Une implémentation (probablement stupide) dans MATLAB (que je n'ai jamais utilisée, donc je ne sais pas si cela va fonctionner. Je viens de googler "déclarer la matrice vectorielle MATLAB" pour obtenir la syntaxe):

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
la source

Ça a l'air super! Je ne sais pas ce que fait la 3e ligne dans votre code (principalement parce que je connais Matlab). Avez-vous des chances de l'écrire dans matlab ou pseudo-code? Je serais bien obligé.

HCAI

x_{1}, \dots, x_{n}

$x_1,\dots,x_n$

t = 1, \dots, n - 1

$t=1,\dots,n-1$

p_{x_{t}, x_{t + 1}}

$p_{x_t,x_{t+1}}$

(p_{i j})

$(p_{ij})$

for

x

$x$

x_{i}

$x_i$

x_{j}

$x_j$

x

$x$

\to

$\to$

Voici mon implémentation 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
la source

user32041 de (publiée en tant que modification au lieu d'un commentaire car il / elle manque de réputation): Comment puis-je contraindre la transitionMatrix du résultat markovchainFit à un data.frame?

chl

Vous pouvez convertir en

d a t a . f r a m e

$data.frame$ en utilisant

a s (m c X, " d a t a . f r a m e ")

$as(mcX,"data.frame")$

Giorgio Spedicato

@GiorgioSpedicato pouvez-vous commenter comment traiter des séquences de longueurs inégales (je ne peux pas concaténer) s'il vous plaît dans votre colis?

HCAI

@HCAI, veuillez consulter la vignette actuelle page 35-36

Giorgio Spedicato

@GiorgioSpedicato merci pour la référence cran.r-project.org/web/packages/markovchain/vignettes/… . J'ai encore n matrices de transition, une pour chaque séquence. Ce que je recherche, c'est un général qui prend en compte toutes les observations de séquence. Y a-t-il quelque chose qui me manque?

HCAI

Voici une façon de le faire dans 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))

Remerciements dus à SomptingGuy: http://www.eng-tips.com/viewthread.cfm?qid=236532

John
la source