Generating annual estimates of N of a binomial distribution

I'm trying to estimate N from a binomial distribution as in https://stats.stackexchange.com/questions/113851/bayesian-estimation-of-n-of-a-binomial-distribution. I have 40 years of counts and an estimated success rate. I'd like to generate an estimate of N for each year. I'm new to JAGS and thought I could get estimates of N by adding an index (n[i]), but I get the following error: "Index out of range taking subset of n." Am I completely off base in what I'm trying to do here? Below is a model and random dataset that I thought would generate annual estimates of N. I do have a model working that generates a single estimate of N, which is just removing the index ([i]) from n in the likelihood (dbin(theta, n)). Thanks in advance for your help.

sink("file.jags")
cat("
    model {

    ## Likelihood    
    for (i in 1:nyear) {
      
      x[i] ~ dbin(theta, n[i])
     
    }

    ## Priors    
    lambda ~ dgamma(0.005, 0.005)    
    theta ~ dbeta(22, 37)
    mu <- lambda/theta
    n ~ dpois(mu)
    
}    
", fill = TRUE)

sink()

# Initial values
inits <- function() { 
  list(
    n = sample(1:100
               , size = 1) 
    , lambda = runif(1, 1, 100)
  )
}

# Parameters monitored
params <- c('theta'
            , 'n'
)

# MCMC settings
ni <- 2000
nb <- 200
nc <- 3

# Bundle data
jags.data <- list(n = c(round(runif(40, 1, 100)))
                  , nyear = 40)

# Fit model
out <- jagsUI(data = jags.data
              , inits = inits
              , parameters.to.save = params
              , model.file = "file.jags"
              , n.chains = nc
              , n.iter = ni
              , n.burnin = nb
              , verbose = TRUE)

# Summarize posteriors
print(out$summary, dig = 2)

EDIT 1: As suggested by @DaveArmstrong, I needed to include the prior for n in the for loop to generate annual priors. Additionally, I removed n from the initial values list because I was receiving a dimension mismatch error in 'setParameters'. Once I removed n from the 'inits' the model ran.

I also switched the prior on lambda from a gamma to a uniform based on my interpretation of Dave's point re: the model not being able to initialize. After I got the model running, I checked to see if the different prior on lambda had an effect and it didn't. The model still ran using the 'dgamma(0.005, 0.005)' but I'm leaving the dunif prior in.

"model {

    # Likelihood    
    for (i in 1:nyear) {
      
      reported_morts[i] ~ dbin(report_rate, total_morts[i])
      
      # Annual prior for total mortality. Necessary to get annual estimates.
      total_morts[i] ~ dpois(mu)
    }

    # Priors
    lambda ~ dunif(1, 100)
    report_rate ~ dbeta(22, 37)
    mu <- round(lambda/report_rate)
    
}"  

# Initial values. Removed "n = sample(1:100, size = 1)".
inits <- function() { 
  list(
    lambda = runif(1, 1, 100)
  )
}
  

You should be able to do it by changing the model to something like this:

mod <- "model {

    ## Likelihood    
    for (i in 1:nyear) {
      
      x[i] ~ dbin(theta, n[i])
      n[i] ~ dpois(mu)
     
    }

    ## Priors    
    lambda ~ dgamma(0.005, 0.005)    
    theta ~ dbeta(22, 37)
    mu <- lambda/theta
}    
"

The problem before was that n was a scalar in the prior, but a vector in the model. You've got to put a prior on each value of n. That model is syntactically correct, but it doesn't initialize because the priors you've set are very far away from the data. If you just get the prior predictive distribution of x, their posterior means are all less than 0.1, but the data for x is integer values between 1 and 100. This generates a 'node inconsistent with parent' error. I assume these are just dummy data, so perhaps your real problem will work alright. Otherwise, the priors will have to be a bit more realistic to work.