hidden markov models - Why the jags result and depmixS4 are sometimes different? -

hidden markov models - Why the jags result and depmixS4 are sometimes different? -

i have info set next simulated data:

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)))

this hmm problem 2 states meaning hidden markoff chain z_t belongs {1,2}. estimate alpha , beta in 2 different states can utilize bundle 'depmixs4' , find maximum likelihood estimates or can utilize mcmc in 'rjags' package.

i expect these 2 estimations same while when run next programme different simulated data, in several times, answers not same , different!!

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]) (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 (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', info = list('x' = simulatedbernollidata$x, 'y' = simulatedbernollidata$y, 'tt' = simulatedbernollidata$tt, t=t, n = n, = 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]

hidden-markov-models mcmc jags