Ces logiciels permettent le calcul de probabilités, la restauration des états cachés et l'estimation des paramètres pour les modèles ci-dessus. Il est également possible de simuler des processus. L'estimation peut se faire par l'algorithme EM, CEM, SEM et EM à la Gibbs ou suivant toute combinaison séquentielle de ces algorithmes. L'utilisateur a la possibilité de choisir les lois d'émission parmi les 14 modèles gaussiens, la loi exponentielle ou les lois non paramétriques pour des observations discrètes. L'utilisateur peut également imposer des contraintes sur la matrice de transition (matrice de type diagonale par bandes) ou sur la loi de l'état initial (loi stationnaire).
Les logiciels développés reprennent les mêmes fonctionnalités que XEMGaus pour l'estimation, à savoir :
Les paramètres du modèle sont alors estimés par
l'algorithme
EM à la Gibbs suivi de l'algorithme EM, en utilisant
simultanément
les deux séquences simulées.
% simulation d'une chaine de Markov cachee a deux etats
» model.fam = 'normal';
% lois d'emissions gaussiennes
» par.A = [ 0.8000 0.2000; 0.4000 0.6000 ];
% matrice de transition
» par.p = [ 0.6667 0.3333];
% distribution stationnaire de A
» par.mu = [-3; 3];
% les deux parametres d'emission : moyenne des lois gaussiennes
» par.S(:,:,1) = 1;
» par.S(:,:,2) = 0.5;
% les deux parametres d'emission : variance des lois gaussiennes
» [x,z] = chainxrnd([500 600], par,model);
» size(x{1})
ans =
500 1
% x{1} est une sequence simulee de longueur 500
% z{1} est la chaine de Markov correspondante
% estimation des parametres
» model.k=2;
% nombre d'etats caches du modele
» model.model='plkI';
% modele stationnaire, gaussien, avec une variance
% dependante de l'etat cache
» algo = {'gibbs','em'};
% algorithmes utilises : EM a la Gibbs puis EM
» cvg = 'xmlORmaxit';
% critere d'arret : stabilisation de la log vraisemblance ou
% depassement du nombre maximal d'iterations
» maxit = 300;
% nombre maximal d'iterations pour l'algorithme EM a la Gibbs
% puis pour EM
» parinit.A = 'random';
» parinit.p = 'random';
» parinit.mu = 'random';
» parinit.S = 'rand. var.';
% initialisation aleatoire de l'algorithme
» nbxem = 3;
% nombre de valeurs initiales pour l'algorithme
» bestpar = chainxem('x',x,'model',model,'algo',algo,'cvg',cvg,...
'maxit',maxit,'nbxem',nbxem,'parinit',parinit);
|------------------------------|
| 2 state(s) - model plkI |
|------------------------------|
| xem 1: [gibbs.9][em.3]
| xem 2: [gibbs.100.147][em.3]
| xem 3: [gibbs.5][em.3]
% estimation des parametres par l'algorithme EM a la Gibbs puis EM
% le nombre d'iterations effectuees est indiquee pour chaque execution
% de EM avec un parametre initial different
» bestpar.A
ans =
0.6173 0.3827
0.1939 0.8061
% estimateur de la matrice de transition
» bestpar.p
ans =
0.3363 0.6637
% estimateur de la matrice de la loi stationnaire
» bestpar.mu
ans =
2.9874
-3.0489
% estimateur des moyennes des parametres d'emission
» bestpar.S
ans(:,:,1) =
0.4586
ans(:,:,2) =
0.9289
% estimateur des variances des parametres d'emission
% modele a lois d'emission "non parametriques" (autrement dit
multinomiales)
% lignes : etats, colonnes : valeurs
par.O = [0.1 0.1 0.4 0.4; 0.3 0.4 0.2 0.1];
model.fam = 'np';
[x,z] = chainxrnd([500 600], par,model);
The computation of probabilities (typically the likelihood), the hidden state restoration and the parameter estimation can be performed for the models above. These models can also be simulated. The estimation is achieved by the EM algorithm or any combination of its variants: CEM (Classification EM), SEM (Stochastic EM) and EM à la Gibbs. The following observation distributions can be chosen: the 14 Gaussian models, the exponential distribution or nonparametric (i.e. multinomial) distributions for observed processes with finite values. The user can also define some constraints on the transition probability matrix (band-diagonal matrix) or on the initial state distribution (stationary distribution).
These routines have the same features than XEMGaus concerning the estimation, namely:
Then the model parameter are estimated using the EM algorithm à
la Gibbs followed by the classical EM algorithm from the two
simulated sequences:
% simulation of a two-state hidden Markov chain
» model.fam = 'normal';
% Gaussian emission distributions
» par.A = [ 0.8000 0.2000; 0.4000 0.6000 ];
% transition probability matrix
» par.p = [ 0.6667 0.3333];
% stationary distribution associated with A
» par.mu = [-3; 3];
% parameters of the emission distributions: mean of the Gaussian distributions
» par.S(:,:,1) = 1;
» par.S(:,:,2) = 0.5;
% parameters of the emission distributions: variance of the Gaussian distributions
» [x,z] = chainxrnd([500 600], par,model);
» size(x{1})
ans =
500 1
% x{1} is a simulated sequence with length 500
% z{1} is the associated Markov chain
% parameter estimation
» model.k=2;
% number of hidden states
» model.model='plkI';
% stationary Gaussian model with state-dependent variance
» algo = {'gibbs','em'};
% chosen algorithms: EM a la Gibbs and EM
» cvg = 'xmlORmaxit';
% stopping criterion: stabilization of the log-likelihood or
% maximal number of iterations
» maxit = 300;
% maximal number of iterations concerning the EM algorithm a la Gibbs
% and then the EM algorithm
» parinit.A = 'random';
» parinit.p = 'random';
» parinit.mu = 'random';
» parinit.S = 'rand. var.';
% random determination of the initial parameters
» nbxem = 3;
% number of initial values of the parameter
» bestpar = chainxem('x',x,'model',model,'algo',algo,'cvg',cvg,...
'maxit',maxit,'nbxem',nbxem,'parinit',parinit);
|------------------------------|
| 2 state(s) - model plkI |
|------------------------------|
| xem 1: [gibbs.9][em.3]
| xem 2: [gibbs.100.147][em.3]
| xem 3: [gibbs.5][em.3]
% parameter estimation using EM a la Gibbs followed by EM
% for each run of the algorithm, the effective number of iterations is printed
% (each run corresponds to a different initial parameter)
» bestpar.A
ans =
0.6173 0.3827
0.1939 0.8061
% estimate of the transition probability matrix
» bestpar.p
ans =
0.3363 0.6637
% estimate of the stationary distribution
» bestpar.mu
ans =
2.9874
-3.0489
% estimate of the emission parameters: means of the Gaussian distributions
» bestpar.S
ans(:,:,1) =
0.4586
ans(:,:,2) =
0.9289
% estimate of the emission parameters: variances of the Gaussian
distributions
% nonparametric (i.e. multinomial) emission
distributions
% rows: states, columns: values
par.O = [0.1 0.1 0.4 0.4; 0.3 0.4 0.2 0.1];
model.fam = 'np';
[x,z] = chainxrnd([500 600], par,model);
Back to my homepage (in French).
Email : |
Laboratoire Jean Kuntzmann |
|