diff --git a/src/turnbull.rs b/src/turnbull.rs index 916e6e3..48cb7ed 100644 --- a/src/turnbull.rs +++ b/src/turnbull.rs @@ -490,18 +490,18 @@ fn do_icm_step(data: &TurnbullData, p: &Vec, s: &Vec, ll_tolerance: f6 let likelihood_obs_new = get_likelihood_obs(data, &s_new); ll_model_new = likelihood_obs_new.iter().map(|l| l.ln()).sum(); - if ll_model_new > ll_model { - // Constrain if required - if let Some(c) = constraint { - let cur_survival_prob = s_new[c.time_index]; - let _ = &mut p_new[0..c.time_index].iter_mut().for_each(|x| *x *= (1.0 - c.survival_prob) / (1.0 - cur_survival_prob)); // Desired failure probability over current failure probability - let _ = &mut p_new[c.time_index..].iter_mut().for_each(|x| *x *= c.survival_prob / cur_survival_prob); - - s_new = p_to_s(&p_new); - let likelihood_obs_new = get_likelihood_obs(data, &s_new); - ll_model_new = likelihood_obs_new.iter().map(|l| l.ln()).sum(); - } + // Constrain if required + if let Some(c) = constraint { + let cur_survival_prob = s_new[c.time_index]; + let _ = &mut p_new[0..c.time_index].iter_mut().for_each(|x| *x *= (1.0 - c.survival_prob) / (1.0 - cur_survival_prob)); // Desired failure probability over current failure probability + let _ = &mut p_new[c.time_index..].iter_mut().for_each(|x| *x *= c.survival_prob / cur_survival_prob); + s_new = p_to_s(&p_new); + let likelihood_obs_new = get_likelihood_obs(data, &s_new); + ll_model_new = likelihood_obs_new.iter().map(|l| l.ln()).sum(); + } + + if ll_model_new > ll_model { return (p_new, s_new, ll_model_new); } @@ -513,7 +513,10 @@ fn do_icm_step(data: &TurnbullData, p: &Vec, s: &Vec, ll_tolerance: f6 step_size_exponent += 1; - if step_size_exponent > 10 { + // FIXME: Nasty hack! + // Limit of 10 should be good enough, but sometimes ICM step does not work well when constrained so the limit is increased + // It would be greatly preferable to simply handle this case more gracefully when constrained + if step_size_exponent > 100 { panic!("ICM fails to increase log-likelihood"); } }