### coefficients 
coef_intercept = Beta ('coef_intercept', 0.0, None, None, 0)
coef_gender = Beta ('coef_gender', 0.0, None, None, 0)
coef_age_resp = Beta ('coef_age_resp', 0.0, None, None, 0)
coef_cattle_vac = Beta ('coef_cattle_vac', 0.0, None, None, 0)
coef_move_herd = Beta ('coef_move_herd', 0.0, None, None, 0)
coef_hh_size = Beta ('coef_hh_size', 0.0, None, None, 0)
coef_tarmac_access = Beta ('coef_tarmac_access', 0.0, None, None, 0)
coef_vet_shop_access = Beta ('coef_vet_shop_access', 0.0, None, None, 0)
coef_pdn_sys = Beta ('coef_pdn_sys', 0.0, None, None, 0)


beta_gender_REF = Beta ('beta_gender_REF', 1, None, None, 0)
beta_age_resp_REF = Beta ('beta_age_resp_REF', 1, None, None, 0)
beta_cattle_vac_REF = Beta ('beta_cattle_vac_REF', 1, None, None, 0)
beta_move_herd_REF = Beta ('beta_move_herd_REF', 1, None, None, 0)
beta_hh_size_REF = Beta ('beta_hh_size_REF', 1, None, None, 0)
beta_tarmac_access_REF = Beta ('beta_tarmac_access_REF', 1, None, None, 0)
beta_vet_shop_access_REF = Beta ('beta_vet_shop_access_REF', 1, None, None, 0)
beta_pdn_sys_REF = Beta ('beta_pdn_sys_REF', 1, None, None, 0)

beta_gender_cl = coef_gender + beta_gender_REF * bioDraws('beta_gender_cl','UNIFORM')
beta_age_resp_cl = coef_age_resp + beta_age_resp_REF * bioDraws('beta_age_resp_cl','NORMAL_MLHS')
beta_cattle_vac_cl = coef_cattle_vac + beta_cattle_vac_REF * bioDraws('beta_cattle_vac_cl','UNIFORM')
beta_move_herd_cl = coef_move_herd + beta_move_herd_REF * bioDraws('beta_move_herd_cl','UNIFORM')
beta_hh_size_cl = coef_hh_size + beta_hh_size_REF * bioDraws('beta_hh_size_cl','NORMAL_MLHS')
beta_tarmac_access_cl = coef_tarmac_access + beta_tarmac_access_REF * bioDraws('beta_tarmac_access_cl','UNIFORM')
beta_vet_shop_access_cl = coef_vet_shop_access + beta_vet_shop_access_REF * bioDraws('beta_vet_shop_access_cl','UNIFORM')
beta_pdn_sys_cl = coef_pdn_sys + beta_pdn_sys_REF * bioDraws('beta_pdn_sys_cl','UNIFORM')

omega = bioDraws('omega', 'NORMAL')
density = dist.normalpdf(omega)
sigma_s1 = Beta('sigma_s1', 0, None, None, 0)
errorComponent = bioDraws('errorComponent', 'NORMAL_HALTON5')
ec_sigma = Beta('ec_sigma', 1, None, None, 0)

LV = (
    coef_intercept 
    + coef_gender * gender 
    + coef_age_resp * age_resp 
    + coef_cattle_vac * cattle_vac 
    + coef_move_herd * move_herd 
    + coef_hh_size * hh_size 
    + coef_tarmac_access * tarmac_access 
    + coef_vet_shop_access * vet_shop_access 
    + coef_pdn_sys * pdn_sys
    + sigma_s1 * omega
    + ec_sigma * errorComponent
)

### Measurement equations 
INTER_VRP1 = Beta('INTER_VRP1', 0, None, None, 0)
INTER_VRP2 = Beta('INTER_VRP2', 0, None, None, 0)
INTER_VRP3 = Beta('INTER_VRP3', 0, None, None, 1)
INTER_VRP4 = Beta('INTER_VRP4', 0, None, None, 0)
INTER_VRP5 = Beta('INTER_VRP5', 0, None, None, 0,)
INTER_VRP6 = Beta('INTER_VRP6', 0, None, None, 0,)

B_VRP1_F1 = Beta('B_VRP1_F1', 0, None, None, 0)
B_VRP2_F1 = Beta('B_VRP2_F1', 0, None, None, 0)
B_VRP3_F1 = Beta('B_VRP3_F1', 1, None, None, 1)
B_VRP4_F1 = Beta('B_VRP4_F1', 0, None, None, 0)
B_VRP5_F1 = Beta('B_VRP5_F1', 0, None, None, 0)
B_VRP6_F1 = Beta('B_VRP6_F1', 0, None, None, 0) 

MODEL_VRP1 = INTER_VRP1 + B_VRP1_F1 * LV
MODEL_VRP2 = INTER_VRP2 + B_VRP2_F1 * LV
MODEL_VRP3 = INTER_VRP3 + B_VRP3_F1 * LV
MODEL_VRP4 = INTER_VRP4 + B_VRP4_F1 * LV
MODEL_VRP5 = INTER_VRP5 + B_VRP5_F1 * LV 
MODEL_VRP6 = INTER_VRP6 + B_VRP6_F1 * LV 

SIGMA_STAR_VRP1 = Beta('SIGMA_STAR_VRP1', 1, 1.0e-5, None, 0)
SIGMA_STAR_VRP2 = Beta('SIGMA_STAR_VRP2', 1, 1.0e-5, None, 0)
SIGMA_STAR_VRP3 = Beta('SIGMA_STAR_VRP3', 1, 1.0e-5, None, 1)
SIGMA_STAR_VRP4 = Beta('SIGMA_STAR_VRP4', 1, 1.0e-5, None, 0)
SIGMA_STAR_VRP5 = Beta('SIGMA_STAR_VRP5', 1, 1.0e-5, None, 0)
SIGMA_STAR_VRP6 = Beta('SIGMA_STAR_VRP6', 1, 1.0e-5, None, 0)

delta_1 = Beta('delta_1', 0.1, 1.0e-5, None, 0)
delta_2 = Beta('delta_2', 0.2, 1.0e-5, None, 0)
tau_1 = -delta_1 - delta_2
tau_2 = -delta_1
tau_3 = delta_1
tau_4 = delta_1 + delta_2

VRP1_tau_1 = (tau_1 - MODEL_VRP1) / SIGMA_STAR_VRP1
VRP1_tau_2 = (tau_2 - MODEL_VRP1) / SIGMA_STAR_VRP1
VRP1_tau_3 = (tau_3 - MODEL_VRP1) / SIGMA_STAR_VRP1
VRP1_tau_4 = (tau_4 - MODEL_VRP1) / SIGMA_STAR_VRP1
IndVRP1 = {
    1: bioNormalCdf(VRP1_tau_1),
    2: bioNormalCdf(VRP1_tau_2) - bioNormalCdf(VRP1_tau_1),
    3: bioNormalCdf(VRP1_tau_3) - bioNormalCdf(VRP1_tau_2),
    4: bioNormalCdf(VRP1_tau_4) - bioNormalCdf(VRP1_tau_3),
    5: 1 - bioNormalCdf(VRP1_tau_4),
} 

P_VRP1 = Elem(IndVRP1, VRP1)

VRP2_tau_1 = (tau_1 - MODEL_VRP2) / SIGMA_STAR_VRP2
VRP2_tau_2 = (tau_2 - MODEL_VRP2) / SIGMA_STAR_VRP2
VRP2_tau_3 = (tau_3 - MODEL_VRP2) / SIGMA_STAR_VRP2
VRP2_tau_4 = (tau_4 - MODEL_VRP2) / SIGMA_STAR_VRP2
IndVRP2 = {
    1: bioNormalCdf(VRP2_tau_1),
    2: bioNormalCdf(VRP2_tau_2) - bioNormalCdf(VRP2_tau_1),
    3: bioNormalCdf(VRP2_tau_3) - bioNormalCdf(VRP2_tau_2),
    4: bioNormalCdf(VRP2_tau_4) - bioNormalCdf(VRP2_tau_3),
    5: 1 - bioNormalCdf(VRP2_tau_4),
}

P_VRP2 = Elem(IndVRP2, VRP2)

VRP3_tau_1 = (tau_1 - MODEL_VRP3) / SIGMA_STAR_VRP3
VRP3_tau_2 = (tau_2 - MODEL_VRP3) / SIGMA_STAR_VRP3
VRP3_tau_3 = (tau_3 - MODEL_VRP3) / SIGMA_STAR_VRP3
VRP3_tau_4 = (tau_4 - MODEL_VRP3) / SIGMA_STAR_VRP3
IndVRP3 = {
    1: bioNormalCdf(VRP3_tau_1),
    2: bioNormalCdf(VRP3_tau_2) - bioNormalCdf(VRP3_tau_1),
    3: bioNormalCdf(VRP3_tau_3) - bioNormalCdf(VRP3_tau_2),
    4: bioNormalCdf(VRP3_tau_4) - bioNormalCdf(VRP3_tau_3),
    5: 1 - bioNormalCdf(VRP3_tau_4),
} 

P_VRP3 = Elem(IndVRP3, VRP3)

VRP4_tau_1 = (tau_1 - MODEL_VRP4) / SIGMA_STAR_VRP4
VRP4_tau_2 = (tau_2 - MODEL_VRP4) / SIGMA_STAR_VRP4
VRP4_tau_3 = (tau_3 - MODEL_VRP4) / SIGMA_STAR_VRP4
VRP4_tau_4 = (tau_4 - MODEL_VRP4) / SIGMA_STAR_VRP4
IndVRP4 = {
    1: bioNormalCdf(VRP4_tau_1),
    2: bioNormalCdf(VRP4_tau_2) - bioNormalCdf(VRP4_tau_1),
    3: bioNormalCdf(VRP4_tau_3) - bioNormalCdf(VRP4_tau_2),
    4: bioNormalCdf(VRP4_tau_4) - bioNormalCdf(VRP4_tau_3),
    5: 1 - bioNormalCdf(VRP4_tau_4),
} 

P_VRP4 = Elem(IndVRP4, VRP4)

VRP5_tau_1 = (tau_1 - MODEL_VRP5) / SIGMA_STAR_VRP5
VRP5_tau_2 = (tau_2 - MODEL_VRP5) / SIGMA_STAR_VRP5
VRP5_tau_3 = (tau_3 - MODEL_VRP5) / SIGMA_STAR_VRP5
VRP5_tau_4 = (tau_4 - MODEL_VRP5) / SIGMA_STAR_VRP5
IndVRP5 = {
    1: bioNormalCdf(VRP5_tau_1),
    2: bioNormalCdf(VRP5_tau_2) - bioNormalCdf(VRP5_tau_1),
    3: bioNormalCdf(VRP5_tau_3) - bioNormalCdf(VRP5_tau_2),
    4: bioNormalCdf(VRP5_tau_4) - bioNormalCdf(VRP5_tau_3),
    5: 1 - bioNormalCdf(VRP5_tau_4),
} 

P_VRP5 = Elem(IndVRP5, VRP5)

VRP6_tau_1 = (tau_1 - MODEL_VRP6) / SIGMA_STAR_VRP6
VRP6_tau_2 = (tau_2 - MODEL_VRP6) / SIGMA_STAR_VRP6
VRP6_tau_3 = (tau_3 - MODEL_VRP6) / SIGMA_STAR_VRP6
VRP6_tau_4 = (tau_4 - MODEL_VRP6) / SIGMA_STAR_VRP6
IndVRP6 = {
    1: bioNormalCdf(VRP6_tau_1),
    2: bioNormalCdf(VRP6_tau_2) - bioNormalCdf(VRP6_tau_1),
    3: bioNormalCdf(VRP6_tau_3) - bioNormalCdf(VRP6_tau_2),
    4: bioNormalCdf(VRP6_tau_4) - bioNormalCdf(VRP6_tau_3),
    5: 1 - bioNormalCdf(VRP6_tau_4),
} 

P_VRP6 = Elem(IndVRP6, VRP6)

INTER_VPC1 = Beta('INTER_VPC1', 0, None, None, 1)
INTER_VPC2 = Beta('INTER_VPC2', 0, None, None, 0)
INTER_VPC3 = Beta('INTER_VPC3', 0, None, None, 0)
INTER_VPC4 = Beta('INTER_VPC4', 0, None, None, 0)
INTER_VPC5 = Beta('INTER_VPC5', 0, None, None, 0)
INTER_VPC6 = Beta('INTER_VPC6', 0, None, None, 0)

B_VPC1_F1 = Beta('B_VPC1_F1', 1, None, None, 1)
B_VPC2_F1 = Beta('B_VPC2_F1', 0, None, None, 0)
B_VPC3_F1 = Beta('B_VPC3_F1', 0, None, None, 0)
B_VPC4_F1 = Beta('B_VPC4_F1', 0, None, None, 0)
B_VPC5_F1 = Beta('B_VPC5_F1', 0, None, None, 0)
B_VPC6_F1 = Beta('B_VPC6_F1', 0, None, None, 0)

MODEL_VPC1 = INTER_VPC1 + B_VPC1_F1 * LV
MODEL_VPC2 = INTER_VPC2 + B_VPC2_F1 * LV
MODEL_VPC3 = INTER_VPC3 + B_VPC3_F1 * LV
MODEL_VPC4 = INTER_VPC4 + B_VPC4_F1 * LV
MODEL_VPC5 = INTER_VPC5 + B_VPC5_F1 * LV 
MODEL_VPC6 = INTER_VPC6 + B_VPC6_F1 * LV

SIGMA_STAR_VPC1 = Beta('SIGMA_STAR_VPC1', 1, 1.0e-5, None, 1)
SIGMA_STAR_VPC2 = Beta('SIGMA_STAR_VPC2', 1, 1.0e-5, None, 0)
SIGMA_STAR_VPC3 = Beta('SIGMA_STAR_VPC3', 1, 1.0e-5, None, 0)
SIGMA_STAR_VPC4 = Beta('SIGMA_STAR_VPC4', 1, 1.0e-5, None, 0)
SIGMA_STAR_VPC5 = Beta('SIGMA_STAR_VPC5', 1, 1.0e-5, None, 0)
SIGMA_STAR_VPC6 = Beta('SIGMA_STAR_VPC6', 1, 1.0e-5, None, 0)

VPC1_tau_1 = (tau_1 - MODEL_VPC1) / SIGMA_STAR_VPC1
VPC1_tau_2 = (tau_2 - MODEL_VPC1) / SIGMA_STAR_VPC1
VPC1_tau_3 = (tau_3 - MODEL_VPC1) / SIGMA_STAR_VPC1
VPC1_tau_4 = (tau_4 - MODEL_VPC1) / SIGMA_STAR_VPC1
IndVPC1 = {
    1: bioNormalCdf(VPC1_tau_1),
    2: bioNormalCdf(VPC1_tau_2) - bioNormalCdf(VPC1_tau_1),
    3: bioNormalCdf(VPC1_tau_3) - bioNormalCdf(VPC1_tau_2),
    4: bioNormalCdf(VPC1_tau_4) - bioNormalCdf(VPC1_tau_3),
    5: 1 - bioNormalCdf(VPC1_tau_4),
} 

P_VPC1 = Elem(IndVPC1, VPC1)

VPC2_tau_1 = (tau_1 - MODEL_VPC2) / SIGMA_STAR_VPC2
VPC2_tau_2 = (tau_2 - MODEL_VPC2) / SIGMA_STAR_VPC2
VPC2_tau_3 = (tau_3 - MODEL_VPC2) / SIGMA_STAR_VPC2
VPC2_tau_4 = (tau_4 - MODEL_VPC2) / SIGMA_STAR_VPC2
IndVPC2 = {
    1: bioNormalCdf(VPC2_tau_1),
    2: bioNormalCdf(VPC2_tau_2) - bioNormalCdf(VPC2_tau_1),
    3: bioNormalCdf(VPC2_tau_3) - bioNormalCdf(VPC2_tau_2),
    4: bioNormalCdf(VPC2_tau_4) - bioNormalCdf(VPC2_tau_3),
    5: 1 - bioNormalCdf(VPC2_tau_4),
} 

P_VPC2 = Elem(IndVPC2, VPC2)

VPC3_tau_1 = (tau_1 - MODEL_VPC3) / SIGMA_STAR_VPC3
VPC3_tau_2 = (tau_2 - MODEL_VPC3) / SIGMA_STAR_VPC3
VPC3_tau_3 = (tau_3 - MODEL_VPC3) / SIGMA_STAR_VPC3
VPC3_tau_4 = (tau_4 - MODEL_VPC3) / SIGMA_STAR_VPC3
IndVPC3 = {
    1: bioNormalCdf(VPC3_tau_1),
    2: bioNormalCdf(VPC3_tau_2) - bioNormalCdf(VPC3_tau_1),
    3: bioNormalCdf(VPC3_tau_3) - bioNormalCdf(VPC3_tau_2),
    4: bioNormalCdf(VPC3_tau_4) - bioNormalCdf(VPC3_tau_3),
    5: 1 - bioNormalCdf(VPC3_tau_4),
} 

P_VPC3 = Elem(IndVPC3, VPC3)

VPC4_tau_1 = (tau_1 - MODEL_VPC4) / SIGMA_STAR_VPC4
VPC4_tau_2 = (tau_2 - MODEL_VPC4) / SIGMA_STAR_VPC4
VPC4_tau_3 = (tau_3 - MODEL_VPC4) / SIGMA_STAR_VPC4
VPC4_tau_4 = (tau_4 - MODEL_VPC4) / SIGMA_STAR_VPC4
IndVPC4 = {
    1: bioNormalCdf(VPC4_tau_1),
    2: bioNormalCdf(VPC4_tau_2) - bioNormalCdf(VPC4_tau_1),
    3: bioNormalCdf(VPC4_tau_3) - bioNormalCdf(VPC4_tau_2),
    4: bioNormalCdf(VPC4_tau_4) - bioNormalCdf(VPC4_tau_3),
    5: 1 - bioNormalCdf(VPC4_tau_4),
} 

P_VPC4 = Elem(IndVPC4, VPC4)

VPC5_tau_1 = (tau_1 - MODEL_VPC5) / SIGMA_STAR_VPC5
VPC5_tau_2 = (tau_2 - MODEL_VPC5) / SIGMA_STAR_VPC5
VPC5_tau_3 = (tau_3 - MODEL_VPC5) / SIGMA_STAR_VPC5
VPC5_tau_4 = (tau_4 - MODEL_VPC5) / SIGMA_STAR_VPC5
IndVPC5 = {
    1: bioNormalCdf(VPC5_tau_1),
    2: bioNormalCdf(VPC5_tau_2) - bioNormalCdf(VPC5_tau_1),
    3: bioNormalCdf(VPC5_tau_3) - bioNormalCdf(VPC5_tau_2),
    4: bioNormalCdf(VPC5_tau_4) - bioNormalCdf(VPC5_tau_3),
    5: 1 - bioNormalCdf(VPC5_tau_4),
} 

P_VPC5 = Elem(IndVPC5, VPC5)

VPC6_tau_1 = (tau_1 - MODEL_VPC6) / SIGMA_STAR_VPC6
VPC6_tau_2 = (tau_2 - MODEL_VPC6) / SIGMA_STAR_VPC6
VPC6_tau_3 = (tau_3 - MODEL_VPC6) / SIGMA_STAR_VPC6
VPC6_tau_4 = (tau_4 - MODEL_VPC6) / SIGMA_STAR_VPC6
IndVPC6 = {
    1: bioNormalCdf(VPC6_tau_1),
    2: bioNormalCdf(VPC6_tau_2) - bioNormalCdf(VPC6_tau_1),
    3: bioNormalCdf(VPC6_tau_3) - bioNormalCdf(VPC6_tau_2),
    4: bioNormalCdf(VPC6_tau_4) - bioNormalCdf(VPC6_tau_3),
    5: 1 - bioNormalCdf(VPC6_tau_4),
} 

P_VPC6 = Elem(IndVPC6, VPC6)

### Choice Model
ASC1 = Beta('ASC1', 0, None, None,1) 
ASC2 = Beta('ASC2', 0, None, None,1)  
ASC3 = Beta('ASC3', 0, None, None,0) 
BETA_PRICE = Beta('BETA_PRICE', 0, None, None, 0)
BETA_SINGDOSE = Beta('BETA_SINGDOSE', 0, None, None, 0) 
BETA_INIT1ANNUAL = Beta('BETA_INITANNUAL', 0, None, None, 0) 
BETA_INIT2ANNUAL = Beta('BETA_INIT2ANNUAL', 0, None, None, 0) 
BETA_POOLED = Beta('BETA_POOLED', 0, None, None, 0) 
BETA_5KM = Beta('BETA_5KM_class', 0, None, None, 0) 
BETA_MORE5KM = Beta('BETA_MORE5KM', 0, None, None, 0)
BETA_NORXN = Beta('BETA_NORXN', 0, None, None, 0) 
BETA_FEVER = Beta('BETA_FEVER', 0, None, None, 0) 
BETA_ABORT = Beta('BETA_ABORT', 0, None, None, 0)
BETA_DAY7 = Beta('BETA_DAY7', 0, None, None, 0) 
BETA_DAY14 = Beta('BETA_DAY14', 0, None, None, 0)
BETA_DAY21 = Beta('BETA_DAY21', 0, None, None, 0)
BETA_NOEXEMPT = Beta('BETA_NOEXEMPT', 0, None, None, 0) 
BETA_CALF6MNTH = Beta('BETA_CALF6MNTH', 0, None, None, 0) 
BETA_OUTBREAK = Beta('BETA_OUTBREAK', 0, None, None, 0)

BETA_PRICE_REF = Beta('BETA_PRICE_REF', 1, None, None, 0)
BETA_SINGDOSE_REF = Beta('BETA_SINGDOSE_REF', 1, None, None, 0)
BETA_INIT1ANNUAL_REF = Beta('BETA_INITANNUAL_REF', 1, None, None, 0)
BETA_INIT2ANNUAL_REF = Beta('BETA_INIT2ANNUAL_REF', 1, None, None, 0)
BETA_POOLED_REF = Beta('BETA_POOLED_REF', 1, None, None, 0)
BETA_5KM_REF = Beta('BETA_5KM_REF', 1, None, None, 0)
BETA_MORE5KM_REF = Beta('BETA_MORE5KM_REF', 1, None, None, 0)
BETA_NORXN_REF = Beta('BETA_NORXN_REF', 1, None, None, 0)
BETA_FEVER_REF = Beta('BETA_FEVER_REF', 1, None, None, 0)
BETA_ABORT_REF = Beta('BETA_ABORT_REF', 1, None, None, 0)
BETA_DAY7_REF = Beta('BETA_DAY7_REF', 1, None, None, 0)
BETA_DAY14_REF = Beta('BETA_DAY14_REF', 1, None, None, 0)
BETA_DAY21_REF = Beta('BETA_DAY21_REF', 1, None, None, 0)
BETA_NOEXEMPT_REF = Beta('BETA_NOEXEMPT_REF', 1, None, None, 0)
BETA_CALF6MNTH_REF = Beta('BETA_CALF6MNTH_REF', 1, None, None, 0)
BETA_OUTBREAK_REF = Beta('BETA_OUTBREAK_REF', 1, None, None, 0)

BETA_PRICE_CL = BETA_PRICE + BETA_PRICE_REF * bioDraws('BETA_PRICE_CL', 'UNIFORM')
BETA_SINGDOSE_CL = BETA_SINGDOSE + BETA_SINGDOSE_REF * bioDraws('BETA_SINGDOSE_CL', 'UNIFORM') 
BETA_INIT1ANNUAL_CL = BETA_INIT1ANNUAL + BETA_INIT1ANNUAL_REF * bioDraws('BETA_INIT1ANNUAL_CL', 'UNIFORM') 
BETA_INIT2ANNUAL_CL = BETA_INIT2ANNUAL + BETA_INIT2ANNUAL_REF * bioDraws('BETA_INIT2ANNUAL_CL', 'UNIFORM') 
BETA_POOLED_CL = BETA_POOLED + BETA_POOLED_REF * bioDraws('BETA_POOLED_CL', 'UNIFORM') 
BETA_5KM_CL = BETA_5KM + BETA_5KM_REF * bioDraws('BETA_5KM_CL', 'UNIFORM') 
BETA_MORE5KM_CL = BETA_MORE5KM + BETA_MORE5KM_REF * bioDraws('BETA_MORE5KM_CL', 'UNIFORM') 
BETA_NORXN_CL = BETA_NORXN + BETA_NORXN_REF * bioDraws('BETA_NORXN_CL', 'UNIFORM') 
BETA_FEVER_CL = BETA_FEVER + BETA_FEVER_REF * bioDraws('BETA_FEVER_CL', 'UNIFORM') 
BETA_ABORT_CL = BETA_ABORT + BETA_ABORT_REF * bioDraws('BETA_ABORT_CL', 'UNIFORM') 
BETA_DAY7_CL = BETA_DAY7 + BETA_DAY7_REF * bioDraws('BETA_DAY7_CL', 'UNIFORM') 
BETA_DAY14_CL = BETA_DAY14 + BETA_DAY14_REF * bioDraws('BETA_DAY14_CL', 'UNIFORM')  
BETA_DAY21_CL = BETA_DAY21 + BETA_DAY21_REF * bioDraws('BETA_DAY21_CL', 'UNIFORM') 
BETA_NOEXEMPT_CL = BETA_NOEXEMPT + BETA_NOEXEMPT_REF * bioDraws('BETA_NOEXEMPT_CL', 'UNIFORM') 
BETA_CALF6MNTH_CL = BETA_CALF6MNTH + BETA_CALF6MNTH_REF * bioDraws('BETA_CALF6MNTH_CL', 'UNIFORM') 
BETA_OUTBREAK_CL = BETA_OUTBREAK + BETA_OUTBREAK_REF * bioDraws('BETA_OUTBREAK_CL', 'UNIFORM')

V1 = (
    ASC1 
    + BETA_PRICE * price_a 
    + BETA_SINGDOSE * single_a 
    + BETA_INIT1ANNUAL * initial_annual_a 
    + BETA_INIT2ANNUAL * intitial2_annual_a 
    + BETA_POOLED * vac_pool_a 
    + BETA_5KM * vac5km_a 
    + BETA_MORE5KM * vacmore5km_a 
    + BETA_NORXN * norxn_a 
    + BETA_FEVER * fever_a 
    + BETA_ABORT * likabort_a 
    + BETA_DAY7 * days7_a 
    + BETA_DAY14 * days14_a 
    + BETA_DAY21 * days21_a 
    + BETA_NOEXEMPT * noexmpt_a 
    + BETA_CALF6MNTH * calve6mth_a 
    + BETA_OUTBREAK * outbreak_a 
    + ec_sigma * errorComponent
)

V2 = (
    ASC2
    + BETA_PRICE * price_b 
    + BETA_SINGDOSE * single_b 
    + BETA_INIT1ANNUAL * initial_annual_b 
    + BETA_INIT2ANNUAL * intitial2_annual_b 
    + BETA_POOLED * vac_pool_b 
    + BETA_5KM * vac5km_b 
    + BETA_MORE5KM * vacmore5km_b 
    + BETA_NORXN * norxn_b 
    + BETA_FEVER * fever_b 
    + BETA_ABORT * likabort_b 
    + BETA_DAY7 * days7_b 
    + BETA_DAY14 * days14_b 
    + BETA_DAY21 * days21_b 
    + BETA_NOEXEMPT * noexmpt_b 
    + BETA_CALF6MNTH * calve6mth_b 
    + BETA_OUTBREAK * outbreak_b
    + ec_sigma * errorComponent
) 

V3 = ASC3 

V = {1: V1, 2: V2, 3: V3} 

condprob =  PanelLikelihoodTrajectory(models.logit(V, None, choice))

loglike = log(MonteCarlo(condlike))

biogeme = bio.BIOGEME(database, loglike, numberOfDraws=250)
biogeme.modelName = 'LatentFullmd01_Ver2'