r语言 - 从 depmixS4 包中获取 HMM 的估计协方差矩阵

Question

>我正在使用 depmixS4 包来存储隐藏的马尔可夫模型。

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

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

Answer 1

在上述示例的扩展中，以下脚本演示如何扩展到三维案例。请注意，输入/输出参数使用方差/协方差矩阵而不是标准差，如上面示例中的名称所示

# 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

Answer 2

我直接询问了软件包的作者，得到了以下答案：

见?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

Sigma1-3 是协方差矩阵的（唯一）值。