# Script generated by the UniPV WinBUGS converter
# 'WinBUGSConverter' copyright (c) BMS - University of Pavia (2016)
# Converter Version: 1.0.0
# Source: PharmML 0.8.1
# Target: WINBUGS 1.4.3
# Run ID: run_id_0
# Model: Generated from MDL. MOG ID: simeoni2004
# File: Simeoni_PAGE_Estimation_BUGS
# Dated: Mon Aug 01 23:35:08 CEST 2016

model {
  for (ind_subj in 1:N_subj) {
    for (ind_t in 1:N_t[ind_subj]) {
      Y_mean[ind_subj,ind_t] <- (WTOT[ind_subj,ind_t] + 0)
      Y_prec[ind_subj,ind_t] <- 1 / ((pow(((CV * WTOT[ind_subj,ind_t])),2))*1)
      Y[ind_subj,ind_t] ~ dnorm(Y_mean[ind_subj,ind_t], Y_prec[ind_subj,ind_t])
      Y_res[ind_subj,ind_t] <- Y[ind_subj,ind_t]-Y_mean[ind_subj,ind_t]
    }
    LAMBDA0[ind_subj] <- LAMBDA0_POP
    LAMBDA1[ind_subj] <- LAMBDA1_POP
    K1[ind_subj] <- K1_POP
    K2[ind_subj] <- K2_POP
    K10[ind_subj] <- K10_POP
    K12[ind_subj] <- K12_POP
    K21[ind_subj] <- K21_POP
    V1[ind_subj] <- V1_POP
    PSI[ind_subj] <- PSI_POP
    W0[ind_subj] <- W0_POP
    der99wb_unipv[ind_subj,1:N_t[ind_subj], 1:6]<- function.model1(grid[ind_subj,1:N_t[ind_subj]], amt[ind_subj,1:N_t[ind_subj]],  rate[ind_subj,1:N_t[ind_subj]], ii[ind_subj,1:N_t[ind_subj]], evid[ind_subj,1:N_t[ind_subj]], cmt[ind_subj,1:N_t[ind_subj]], addl[ind_subj,1:N_t[ind_subj]], ss[ind_subj,1:N_t[ind_subj]], theta[ind_subj,])
    theta[ind_subj,1] <- K21_POP
    theta[ind_subj,2] <- K10_POP
    theta[ind_subj,3] <- K12_POP
    theta[ind_subj,4] <- PSI_POP
    theta[ind_subj,5] <- logLAMBDA0_POP
    theta[ind_subj,6] <- logLAMBDA1_POP
    theta[ind_subj,7] <- logK2_POP
    theta[ind_subj,8] <- V1_POP
    theta[ind_subj,9] <- logK1_POP
    theta[ind_subj,10] <- 0 #IC Q1
    theta[ind_subj,11] <- 0 #IC Q2
    theta[ind_subj,12] <- W0[ind_subj] #IC X1
    theta[ind_subj,13] <- 0 #IC X2
    theta[ind_subj,14] <- 0 #IC X3
    theta[ind_subj,15] <- 0 #IC X4
    theta[ind_subj,16] <- 1 #F Q1
    theta[ind_subj,17] <- 1 #F Q2
    theta[ind_subj,18] <- 1 #F X1
    theta[ind_subj,19] <- 1 #F X2
    theta[ind_subj,20] <- 1 #F X3
    theta[ind_subj,21] <- 1 #F X4
    theta[ind_subj,22] <- 0 #TLAG Q1
    theta[ind_subj,23] <- 0 #TLAG Q2
    theta[ind_subj,24] <- 0 #TLAG X1
    theta[ind_subj,25] <- 0 #TLAG X2
    theta[ind_subj,26] <- 0 #TLAG X3
    theta[ind_subj,27] <- 0 #TLAG X4
    for (ind_q in 1:N_t[ind_subj]) {
      Q1[ind_subj,ind_q] <- der99wb_unipv[ind_subj,ind_q,1] + 0
      Q2[ind_subj,ind_q] <- der99wb_unipv[ind_subj,ind_q,2] + 0
      X1[ind_subj,ind_q] <- der99wb_unipv[ind_subj,ind_q,3] + W0[ind_subj]
      X2[ind_subj,ind_q] <- der99wb_unipv[ind_subj,ind_q,4] + 0
      X3[ind_subj,ind_q] <- der99wb_unipv[ind_subj,ind_q,5] + 0
      X4[ind_subj,ind_q] <- der99wb_unipv[ind_subj,ind_q,6] + 0
      C[ind_subj,ind_q] <- (Q1[ind_subj,ind_q] / V1[ind_subj])
      WTOT[ind_subj,ind_q] <- (((X1[ind_subj,ind_q] + X2[ind_subj,ind_q]) + X3[ind_subj,ind_q]) + X4[ind_subj,ind_q])
    }
  }
  MU_LAMBDA0_POP <- (-(1.4))
  MU_LAMBDA1_POP <- (-(1.3))
  MU_K1_POP <- (-(1.2))
  MU_K2_POP <- (-(3.2))
  MU_W0_POP <- (-(3.9))
  SIGMA_LAMBDA0_POP <- 0.8
  SIGMA_LAMBDA1_POP <- 1.2
  SIGMA_K1_POP <- 1.8
  SIGMA_K2_POP <- 1.4
  SIGMA_W0_POP <- 3.2
  a_CV2 <- 195
  b_CV2 <- 1.4
  K10_POP <- (0.868 * 24)
  K12_POP <- (0.006 * 24)
  K21_POP <- (0.0838 * 24)
  V1_POP <- 0.81
  PSI_POP <- 20
  LAMBDA0_POP <- (exp(logLAMBDA0_POP))
  LAMBDA1_POP <- (exp(logLAMBDA1_POP))
  K1_POP <- (exp(logK1_POP))
  K2_POP <- (exp(logK2_POP))
  W0_POP <- (exp(logW0_POP))
  CV <- (sqrt((1 / invCV2)))
  # PRIORS
  invCV2 ~ dgamma(a_CV2, invCV2_rate)
  invCV2_rate <- 1/b_CV2
  logLAMBDA1_POP ~ dnorm(MU_LAMBDA1_POP, logLAMBDA1_POP_prec)
  logLAMBDA1_POP_prec <- (1/(SIGMA_LAMBDA1_POP))
  logW0_POP ~ dnorm(MU_W0_POP, logW0_POP_prec)
  logW0_POP_prec <- (1/(SIGMA_W0_POP))
  logLAMBDA0_POP ~ dnorm(MU_LAMBDA0_POP, logLAMBDA0_POP_prec)
  logLAMBDA0_POP_prec <- (1/(SIGMA_LAMBDA0_POP))
  logK2_POP ~ dnorm(MU_K2_POP, logK2_POP_prec)
  logK2_POP_prec <- (1/(SIGMA_K2_POP))
  logK1_POP ~ dnorm(MU_K1_POP, logK1_POP_prec)
  logK1_POP_prec <- (1/(SIGMA_K1_POP))
}
