2

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

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]
4

1 回答 1

0

您在 JAGS 代码中遇到了标签切换问题，也就是说，状态不限z[i]=1于较低的后验值Beta和z[i]=2较高的Beta. 因此，对于 MCMC 的每次迭代，它们都可以切换。有几种方法可以解决这个问题。其中之一是部分重新排序，即对于每次 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参数进行排序。选择用于部分重新排序的参数取决于问题。

于 2019-05-02T14:38:40.697 回答