model { for(j in 1:M){ sigmay[j] <- Sigma[1] + Sigma[2]*sst[j] + Sigma[3]*nao[j] xiy[j] <- Xi[1] + Xi[2]*soi[j] log(rate[j]) <- Rate[1] + Rate[2]*sst[j] + Rate[3]*nao[j] + Rate[4]*soi[j] + Rate[5]*ssn[j] L[j] ~ dpois(rate[j]) } for(i in 1:N){ sigma1[i] <- sigmay[offset[i]] # sigma[i] <- step(sigma1[i] - minval)*sigma1[i] + (1 - step(sigma1[i] - minval))*minval sigma[i] <- exp(sigma1[i]) xi[i] <- xiy[offset[i]] ones[i] ~ dbern(p[i]) p[i] <- C/(sigma[i]*pow(w[i], (1/xi[i] + 1)))*step(z[i] - epsilon) #likelihood z[i] <- (1 + xi[i]/sigma[i]*(y[i] - u)) w[i] <- max(z[i], epsilon) } Xi[1] ~ dunif(-1, 1) Xi[2] ~ dunif(-100, 100) Sigma[1] ~ dunif(-100, 100) Sigma[2] ~ dunif(-100, 100) Sigma[3] ~ dunif(-100, 100) Rate[1] ~ dunif(-100, 100) Rate[2] ~ dunif(-100, 100) Rate[3] ~ dunif(-100, 100) Rate[4] ~ dunif(-100, 100) Rate[5] ~ dunif(-100, 100) } data { for(i in 1:N){ ones[i] <- 1 offset[i] <- Yr[i] - Yr0 + 1 } # C <- 1/1000 # epsilon <- 1.0E-12 }