Получение оценочных ковариационных матриц для HMM из пакета depmixS4

В расширении приведенного выше примера следующий сценарий показывает, как расширить его до трехмерного случая. Помните, что входные / выходные параметры используют матрицу дисперсии / ковариации, а не стандартные отклонения, как указано в названиях в приведенном выше примере.

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

Последняя команда дает результат (только последняя часть)

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

Расчетная матрица дисперсии / ковариации с полями

 s_11 s_12 s_13
      s_22 s_23 
           s_33

(Пустые поля могут быть заполнены симметрично) отображаются в выводе как Re1.Sigma1, ..., Re1.Sigma6

Я спросил напрямую у автора пакета и получил следующий ответ:

См. ?makeDepmix

если вы запустите пример с многомерными данными, например:

# 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

Сигма1-3 - (уникальные) значения ковариационной матрицы.