Блог

У меня есть набор данных, подобный следующим моделируемым данным:

 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} . Чтобы оценить альфа и бета в двух разных состояниях, я могу использовать пакет ‘depmixS4’ и найти оценки максимального правдоподобия или я могу использовать MCMC в пакете ‘rjags’.

Я ожидаю, что эти две оценки будут почти одинаковыми, в то время как при запуске следующей программы для разных моделируемых данных в несколько раз ответы не совпадают и сильно отличаются!!

 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]
  

У вас возникла проблема с переключением меток в вашем коде 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 параметрам. Выбор того, какой параметр использовать для частичного изменения порядка, зависит от проблемы.