Obtaining estimated covariance matrices for HMM from depmixS4 package

Question

I am using the depmixS4 package for hidden markov models.

I am applying it on a joint multivariate gaussian distribution for n vectors with m states. My question concerns obtaining the estimated parameters - in particular:

My question
How can I obtain the estimated m covariance matrices?

Answer 1

In extension of the above example, the following script shows how to extend to the 3-dimensional case. Beware that the input / output parameters use the variance/covariance matrix and not standard deviations, as suggested by the names in the above example

# DepMixS4- Multivariate Normal
# Exmaple from the Help Page extended to three dimensions        

library(depmixS4)        

# generate data from two different 3-dim normal distributions
#mean 
m1 <- c(0,0,1) 
# the first one has equal covariance of size 0.3
# Sigma1 denotes the covariance matrix
Sigma1 <- matrix(c(2,0.3,0.3,0.3,1,0.3,0.3,0.3,1),nrow=3)
#mean
m2 <- c(1,0,0)
# the second one has equal covariance of size -0.3
Sigma2 <- matrix(c(2,-0.3,-0.3,-0.3,1,-0.3,-0.3,-0.3,1),nrow=3)        

y1 <- mvrnorm(1000,m1,Sigma1)
y2 <- mvrnorm(1000,m2,Sigma2)
# Check that Sigma1_11 is indeed variance:
# The following gives approx 2, as it should be
var(y1[,1])    


# this creates data with a single change point
y <- rbind(y1,y2)        

# now use makeDepmix to create a depmix model for this 3-dim normal timeseries
# response is a 2-dim list of response models. 
rModels <-  list()
rModels[[1]] <- list(MVNresponse(y~1))
rModels[[2]] <- list(MVNresponse(y~1))        

trstart=c(0.9,0.1,0.1,0.9)        

transition <- list()
transition[[1]] <- transInit(~1,nstates=2,data=data.frame(1),pstart=c(trstart[1:2]))
transition[[2]] <- transInit(~1,nstates=2,data=data.frame(1),pstart=c(trstart[3:4]))        

instart=runif(2)
inMod <- transInit(~1,ns=2,ps=instart,data=data.frame(1))        

mod <- makeDepmix(response=rModels,transition=transition,prior=inMod)        

fm3 <- fit(mod,emc=em.control(random=FALSE))      

summary(fm3)    

The last command gives output (only the last part)

Response parameters 
    Re1.coefficients1 Re1.coefficients2 Re1.coefficients3 Re1.Sigma1 Re1.Sigma2 Re1.Sigma3 Re1.Sigma4 Re1.Sigma5 Re1.Sigma6
St1        0.92688362        0.04410173      -0.004027074   2.042186 -0.3347858 -0.2467676   1.066815 -0.3113216  0.9507479
St2        0.03676585        0.05259029       1.022748735   2.037637  0.4235656  0.3856682   1.146025  0.3401500  0.9763637

The estimated variance/covariancematrix with fields

 s_11 s_12 s_13
      s_22 s_23 
           s_33

(Empty fields can be filled in by symmetry) are given in the output as Re1.Sigma1, ... , Re1.Sigma6

Answer 2

I asked the authour of the package directly and got the following answer:

See ?makeDepmix

if you run the example with multivariate data like so:

# generate data from two different multivariate normal distributions
m1 <- c(0,1)
sd1 <- matrix(c(1,0.7,.7,1),2,2)
m2 <- c(1,0)
sd2 <- matrix(c(2,.1,.1,1),2,2)
set.seed(2)
y1 <- mvrnorm(50,m1,sd1)
y2 <- mvrnorm(50,m2,sd2)
# this creates data with a single change point
y <- rbind(y1,y2)

# now use makeDepmix to create a depmix model for this bivariate normal timeseries
rModels <-  list()
rModels[[1]] <- list(MVNresponse(y~1))
rModels[[2]] <- list(MVNresponse(y~1))

trstart=c(0.9,0.1,0.1,0.9)

transition <- list()
transition[[1]] <- transInit(~1,nstates=2,data=data.frame(1),pstart=c(trstart[1:2]))
transition[[2]] <- transInit(~1,nstates=2,data=data.frame(1),pstart=c(trstart[3:4]))

instart=runif(2)
inMod <- transInit(~1,ns=2,ps=instart,data=data.frame(1))

mod <- makeDepmix(response=rModels,transition=transition,prior=inMod)

fm2 <- fit(mod,emc=em.control(random=FALSE))
The output gives this: 

> summary(fm2)
Initial state probabilties model 
            St1     St2
(Intercept)   0 -10.036

Transition matrix 
       toS1 toS2
fromS1 0.98 0.02
fromS2 0.00 1.00

Response parameters 
    Re1.coefficients1 Re1.coefficients2 Re1.Sigma1 Re1.Sigma2 Re1.Sigma3
St1             0.125             1.024      1.346      0.873      1.272
St2             0.716             0.100      2.068      0.285      0.888

Sigma1-3 are the (unique) values of the covariance matrix.