library(R2jags)

source("../DATA/simNmix.R")

# Bundle data
#attach(simdata)

win.data <- list(y = y, R = R, T = T, x = x, gID = gID, nG = nG)

# Inits function
Nst <- apply(y, 1, max) + 1
inits <- function(){list(N = Nst, sigma = rlnorm(1))}
  
# Parameters to be estimated
params <- c("lambda", "loglam", "p0", "p1", "sigma", "logsigma","u", "totalN")

# Could include "N" to get local population size

# MCMC settings (takes about 13 min on a slow laptop)
nc <- 3
nb <- 5000
ni <- 10000
nt <- 5

# MCMC test settings
if(F)
{
  nc <- 1 
  nb <- 2
  ni <- 10
  nt <- 5
  }



outJ <- jags(win.data, inits, params, "Nmix.txt", n.chains=nc,n.iter=ni, n.burn = nb, n.thin=nt)


if(F)
{
results <- list(obj=NULL,                     ## negative log L
                coef=coef,                   ## coefficients
                sd=sd,                       ## standard deviations 
                confint.quad=ci.q,           ## Nx2 - matrix 95% ci 
                confint.profile=NULL,        
                confint.mcmc=ci.mc,           ## 0.025 0.975 quantiles
                time=c(fit=time,profile=NULL),  
                terminfo=c(maxgrad=NULL,eratio=NULL), ## MLE termination info
                convinfo=NULL       ## MCMC convergence info
                )

save(results,mcmcall,file="fit.RData")
}