# 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: example332
# File: example332
# Dated: Thu Aug 04 19:11:58 CEST 2016

model {
  for (ind_subj in 1:N_subj) {
    for (ind_t in 1:N_t[ind_subj]) {
      C[ind_subj,ind_t] <- ((DOSE[ind_subj,1] / V[ind_subj]) * (exp(((-(K[ind_subj])) * grid[ind_subj,ind_t]))))
      Y_mean[ind_subj,ind_t] <- (C[ind_subj,ind_t] + 0)
      Y[ind_subj,ind_t] ~ dnorm(Y_mean[ind_subj,ind_t], Y_prec[ind_subj])
      Y_res[ind_subj,ind_t] <- Y[ind_subj,ind_t]-Y_mean[ind_subj,ind_t]
    }
    Y_prec[ind_subj] <- 1 / ((pow((SD_ADD[ind_subj]),2))*1)
    K[ind_subj] <- (exp(DIST_K_V[1]))
    V[ind_subj] <- (exp(DIST_K_V[2]))
    SD_ADD[ind_subj] <- (1 / (sqrt(TAU[ind_subj])))
    eta_T[ind_subj] ~ dnorm(0, eta_T_prec)
    log(TAU[ind_subj])<- log(POP_T) + eta_T[ind_subj]
  }
  DIST_K_V[1:2] ~ dmnorm(DIST_K_V_mean[],DIST_K_V_prec[,])
  DIST_K_V_mean[1] <- (log(POP_K))
  DIST_K_V_mean[2] <- (log(POP_V))
  DIST_K_V_prec[1:2,1:2] <- inverse(OMEGA_P[,])
  eta_T_prec <- (1/(OMEGA_T))
  MU_POP_K <- 0.1
  MU_POP_V <- 8
  SIGMA_POP_P[1,1] <- 0.1
  SIGMA_POP_P[1,2] <- 5.0E-5
  SIGMA_POP_P[2,1] <- 5.0E-5
  SIGMA_POP_P[2,2] <- 0.1
  R[1,1] <- 1
  R[1,2] <- 0.1
  R[2,1] <- 0.1
  R[2,2] <- 1
  rho <- 5
  a_OMEGA_T <- 1
  b_OMEGA_T <- 0.01
  a_POP_T <- 1
  b_POP_T <- 0.01
  lMU_POP_V <- (log(MU_POP_V))
  POP_K <- (exp(lPOP_P[1]))
  POP_V <- (exp(lPOP_P[2]))
  OMEGA_T <- (1 / TAU_T)
  OMEGA_P[1:2,1:2] <- inverse(invOMEGA_P[,])
  # PRIORS
  TAU_T ~ dgamma(a_OMEGA_T, TAU_T_rate)
  TAU_T_rate <- b_OMEGA_T
  lPOP_P[1:2] ~ dmnorm(lPOP_P_mean[],lPOP_P_prec[,])
  lPOP_P_mean[1] <- (log(MU_POP_K))
  lPOP_P_mean[2] <- lMU_POP_V
  lPOP_P_prec[1:2,1:2] <- inverse(SIGMA_POP_P[,])
  POP_T ~ dgamma(a_POP_T, POP_T_rate)
  POP_T_rate <- b_POP_T
  invOMEGA_P[1:2,1:2] ~ dwish(invOMEGA_P_inverseScaleMatrix[,],rho)
  invOMEGA_P_inverseScaleMatrix[1:2,1:2] <- inverse(R[,])
}
