从 depmixS4 包中获取 HMM 的估计协方差矩阵

Question

我正在为隐马尔可夫模型使用 depmixS4 包。

我将其应用于具有 m 个状态的 n 个向量的联合多元高斯分布。我的问题涉及获取估计参数 - 特别是：

我的问题
如何获得估计的m个协方差矩阵？

Answer 1

我直接问了包的作者，得到了如下答案：

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.

Answer 2

在上述示例的扩展中，以下脚本显示了如何扩展到 3 维情况。请注意，输入/输出参数使用 variance/covariance 矩阵而不是标准偏差，如上例中的名称所示

# 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

估计的 variance/covariance 具有字段的矩阵

 s_11 s_12 s_13
      s_22 s_23 
           s_33

(空字段可以对称填充) 在输出中给出为 Re1.Sigma1, ..., Re1.Sigma6

从 depmixS4 包中获取 HMM 的估计协方差矩阵

Obtaining estimated covariance matrices for HMM from depmixS4 package

r

hidden-markov-models