Why the parameter estimation of my beta-binomial model using jags differ from maximum likelihood estimations

Question

I have a beta-binomial model like this

where $B$ is the beta function.

I want to estimate the parameters $\theta_1,\theta_2,\ldots,\theta_5$.

I used a Maximum likelihood method:

BBlikelihood = function(theta){
    k = theta[1]/(1+exp(theta[2]*x+theta[3]))+theta[4]
    mu = exp(k)/(1+exp(k))
    if(theta[5]<0) theta[5]=0
    return(-sum(lchoose(n,y)+
                   lbeta(y+mu*theta[5],n-y+(1-mu)*theta[5]) -
                   lbeta(mu*theta[5],(1-mu)*theta[5])))
}

theta.init = c(8,100,-3,-6,1)
BBtheta = optim(theta.init,BBlikelihood,control=list(maxit=1000))  
BBtheta$par

Then, I used MCMC as below:

HastingsBBlikelihood = function(theta,theta.prime){
    k = theta[1]/(1+exp(theta[2]*x+theta[3]))+theta[4]
    mu = exp(k)/(1+exp(k))
    k.prime = theta.prime[1]/(1+exp(theta.prime[2]*x+theta.prime[3]))+theta.prime[4]
    mu.prime = exp(k.prime)/(1+exp(k.prime))
    if(theta.prime[5]<0) 
        return(-Inf) 
    sum( ( lbeta(y+mu.prime*theta.prime[5],n-y+(1-mu.prime)*theta.prime[5])
          - lbeta(mu.prime*theta.prime[5],(1-mu.prime)*theta.prime[5]))
        - ( lbeta(y+mu*theta[5],n-y+(1-mu)*theta[5])
           lbeta(mu*theta[5],(1-mu)*theta[5])))
}

SigmaBB = list(0.1*diag(5))
Sigma.i = 1
thetaBB.mcmc = matrix(c(8,100,-3,-6,1),1,5)
for(i in 2:25000){
    if((i %% 1000)==0){ 
        Sigma.i <- Sigma.i + 1; 
        SigmaBB[[Sigma.i]] <- 2.38^2*var(thetaBB.mcmc)/5
    }
    thetaBB.mcmc = rbind(thetaBB.mcmc,
                            mvrnorm(n=1,
                                       thetaBB.mcmc[i-1,],
                                       Sigma=SigmaBB[[Sigma.i]]))
    if(log(runif(1))>HastingsBBlikelihood(thetaBB.mcmc[i-1,],thetaBB.mcmc[i,]))
        thetaBB.mcmc[i,] <- thetaBB.mcmc[i-1,]
}

the MCMC estimations agree with maximum likelihood estimations.

Finally, I used jags to estimate the parameters as the following:

  cat("model {

    # Parameters
    theta[1] ~ dunif(-1.e10, 1.e10)   
    theta[2] ~ dunif(-1.e10, 1.e10)
    theta[3] ~ dunif(-1.e10, 1.e10)
    theta[4] ~ dunif(-1.e10, 1.e10)
    theta[5] ~ dunif(0, 1.e10)

    # Observations
    for (i in 1:TT){

    k[i] <- theta[1]/(1+exp(theta[2]*x[i]+theta[3]))+theta[4]
    mu[i] <- exp(k[i])/(1+exp(k[i]))

    p[i] ~ dbeta(mu[i]*theta[5],(1-mu[i])*theta[5])     

    y[i] ~ dbin(p[i],n[i])

  }
    }",
    file="BetaBinom.bug")    

library("rjags")
TT<-length(y)
jags <- jags.model('BetaBinom.bug', data = list('x'=x, 'y'=y, 'n'=n , TT=TT ))

res <- coda.samples(jags,c('theta'),1000)
res.median = apply(res[[1]],2,median)
res.median
res.mean = apply(res[[1]],2,mean)
res.mean
res.sd = apply(res[[1]],2,sd)
res.sd

The estimations in this way are very strange and not as the previous values. I wonder what is wrong in the jags function!? Thanks for any comment or suggestion in advance.

https://ehc.ac/p/mcmc-jags/discussion/610037/thread/dc35eac1/#4823