为啥 jags 结果和 depmixS4 有时不同？

Posted 2023-02-19

【中文标题】为啥 jags 结果和 depmixS4 有时不同？

【英文标题】：Why the jags result and depmixS4 are sometimes different?为什么 jags 结果和 depmixS4 有时不同？

【发布时间】：2014-08-15 06:05:58

【问题描述】：

我有一个类似以下模拟数据的数据集：

Pi = matrix(c(0.9,0.1,0.3,0.7),2,2,byrow=TRUE)
delta = c(.5,.5)

z = sample(c(1,2),1,prob=delta)
T = 365
for( t in 2:T)
  z[t] = sample(x=c(1,2),1,prob=Pi[z[t-1],])


x <- sample(x=seq(-1, 1.5, length.out=T),T,replace=TRUE)

alpha = c(-1, -3.2)
Beta = c(-4,3)

y<-NA
for(i in 1:T)
  y[i] = rbinom(1,size=10,prob=1/(1+exp(-Beta[z[i]]*x[i]-alpha[z[i]])))


SimulatedBinomData <- data.frame('y' = y, 'x' = x , size=rep(10,T), 'z' = z)

yy<-NA
xx<-NA
for(i in 1:dim(SimulatedBinomData)[1])
  yy<-c(yy,c(rep(1,SimulatedBinomData$y[i]),rep(0,(SimulatedBinomData$size[i]-SimulatedBinomData$y[i]))))
  xx<-c(xx,rep(SimulatedBinomData$x[i],SimulatedBinomData$size[i]))

yy<-yy[-1]
xx<-xx[-1]

SimulatedBernolliData<-data.frame(y=yy,x=xx, tt=rep(c(1:T),rep(10,T)))

这是一个具有两个状态的 HMM 问题，这意味着隐马尔可夫链 z_t 属于 1,2。要估计两种不同状态下的 alpha 和 Beta，我可以使用包“depmixS4”并找到最大似然估计值，或者我可以使用“rjags”包中的 MCMC。

我希望这两个估计值几乎相同，而当我针对不同的模拟数据运行以下程序时，有好几次，答案都不相同而且非常不同！

library("rjags")
library("depmixS4")

mod <- depmix(cbind(y,(size-y))~x, data=SimulatedBinomData, nstates=2, family=binomial(logit))
fm <- fit(mod)
getpars(fm)

n<-length(SimulatedBernolliData$y)
T<-max(SimulatedBernolliData$tt)                 

cat("model 

    # Transition Probability
    Ptrans[1,1:2] ~ ddirch(a)
    Ptrans[2,1:2] ~ ddirch(a)    

    # States
    Pinit[1] <- 0.5 #failor
    Pinit[2] <- 0.5 #success
    state[1] ~ dbern(Pinit[2])

    for (t in 2:T) 
    state[t] ~ dbern(Ptrans[(state[t-1]+1),2])
    

    # Parameters
    alpha[1] ~ dunif(-1.e10, 1.e10)
    alpha[2] ~ dunif(-1.e10, 1.e10)

    Beta[1] ~ dunif(-1.e10, 1.e10)
    Beta[2] ~ dunif(-1.e10, 1.e10)

    # Observations
    for (i in 1:n)
    z[i] <- state[tt[i]] 
    y[i] ~ dbern(1/(1+exp(-(alpha[(z[i]+1)]+Beta[(z[i]+1)]*x[i]))))
    
    ",
file="LeftBehindHiddenMarkov.bug")

jags <- jags.model('LeftBehindHiddenMarkov.bug', data = list('x' =  SimulatedBernolliData$x, 'y' = SimulatedBernolliData$y, 'tt' = SimulatedBernolliData$tt, T=T, n = n, a = c(1,1) ))
res <- coda.samples(jags,c('alpha', 'Beta', 'Ptrans','state'),1000)
res.median = apply(res[[1]],2,median)
res.median[1:8]
res.mean = apply(res[[1]],2,mean)
res.mean[1:8]
res.sd = apply(res[[1]],2,sd)
res.sd[1:8]
res.mode = apply(res[[1]],2,function(x)as.numeric(names(table(x))
[which.max(table(x))]) )
res.mode[1:8]

【参考方案1】：

您在 JAGS 代码中遇到了标签切换问题，即状态 z[i]=1 没有限制到 Beta 和 z[i]=2 的较低后验值到较高的 Beta。因此，对于 MCMC 的每次迭代，它们都可以切换。有several ways to solve this problem。其中之一是部分重排序，即对于每次 MCMC 迭代，为 Beta 绘制两个独立的值并对其进行排序，使得 Beta[1] Beta[2]。

你可以通过替换来做到这一点

Beta[1] ~ dunif(-1.e10, 1.e10)
Beta[2] ~ dunif(-1.e10, 1.e10)

为

Beta[1:2] <- sort(Betaaux)
Betaaux[1] ~ dunif(-1.e10, 1.e10)
Betaaux[2] ~ dunif(-1.e10, 1.e10)

当然，也可以在alpha 参数上进行排序。选择用于部分重新排序的参数取决于问题。