data{ int N; int N_id; int id[N]; real ce[N]; real high[N]; real low[N]; real p[N]; } parameters{ real rho; real alpha; real beta; real tau; vector[3] sigma_sub; cholesky_factor_corr[3] L_Rho_sub; } transformed parameters{ vector[N_id] rho_sub; vector[N_id] alpha_sub; vector[N_id] beta_sub; matrix[3,3] Rho_sub; matrix[N_id,3] v_N_id; rho_sub = col(v_N_id,1); alpha_sub = col(v_N_id,2); beta_sub = col(v_N_id,3); Rho_sub = L_Rho_sub * L_Rho_sub'; } model{ vector[N] sigma; vector[N] wp; vector[N] u_low; vector[N] u_high; vector[N] mu; tau ~ normal( 0.2 , 0.1 ); alpha ~ normal( 0.7 , 0.25 ); beta ~ normal( 1 , 0.25 ); rho ~ normal( 0 , 0.1 ); sigma_sub ~ std_normal( ); L_Rho_sub ~ lkj_corr_cholesky( 3 ); to_vector(v_N_id) ~ multi_normal_cholesky( rep_vector(0,3) , diag_pre_multiply(sigma_sub, L_Rho_sub)); for ( i in 1:N ) { u_high[i] = (1 - exp(-(rho + rho_sub[id[i]]) * high[i]))/(rho + rho_sub[id[i]]); u_low[i] = (1 - exp(-(rho + rho_sub[id[i]]) * low[i]))/(rho + rho_sub[id[i]]); wp[i] = (exp(-(beta + beta_sub[id[i]])*(-log(p[i]))^(alpha + alpha_sub[id[i]]))); sigma[i] = (tau) * (high[i] - low[i]); mu[i] = (log(-(((wp[i] * u_high[i] + (1 - wp[i]) * u_low[i]) * (rho + rho_sub[id[i]])) - 1)))/(-(rho + rho_sub[id[i]])); } ce ~ normal( mu , sigma ); }