model{
  
  # prior distributions
  p ~ dunif(0, 1)
  alpha ~ dunif(low_a, high_a)
  lambda ~ dunif(low_l, high_l)
  
  ### CALCULATING STABLE STAGE DISTRIBUTION
  w[1] <- 1  # proportion to stable stage
  g[1] <- p  # survive to next stage class
  z[1] <- 0  # remain in current stage class
  
  for(s in 2:(S - 1)){  # loop through all stage classes
    g[s] <- p
    z[s] <- 0
    w[s] <- g[s - 1] / (lambda - z[s]) * w[s - 1]
  }
  
  g[5] <- 0
  z[5] <- p
  w[5] <- g[4] / (lambda - z[5]) * w[4]
  
  # stable stage distribution
  C <- w / sum(w)
  
  ### EXPECTED COUNT IN EACH STAGE CLASS
  ey[2] <- p * C[1]
  ey[3] <- p * C[2]
  ey[4] <- p * C[3]
  ey[5] <- p * (C[4] + C[5])
  ey[1] <- alpha * (lambda - sum(ey[2:5]))
  
  ### LIKELIHOOD
  y ~ dmulti(ey, N)

}