From 307aff6f14bb5c42911f501da6263818264374fa Mon Sep 17 00:00:00 2001 From: RunasSudo Date: Mon, 25 Dec 2023 22:21:50 +1100 Subject: [PATCH] turnbull: Disregard ICM step when computing likelihood-ratio confidence intervals ICM step performance is heavily degraded when constraints are required It is much faster to rely on the EM step alone 1275% speedup! --- src/turnbull.rs | 34 +++++++++++++++++++++------------- 1 file changed, 21 insertions(+), 13 deletions(-) diff --git a/src/turnbull.rs b/src/turnbull.rs index cf7b7bf..3c956e9 100644 --- a/src/turnbull.rs +++ b/src/turnbull.rs @@ -312,6 +312,8 @@ fn fit_turnbull_estimator(data: &TurnbullData, progress_bar: ProgressBar, max_it // ------- // EM step + // TODO: Do EM step multiple times per ICM step? + let p_after_em = do_em_step(data, &p, &s, &constraint); let s_after_em = p_to_s(&p_after_em); @@ -324,13 +326,18 @@ fn fit_turnbull_estimator(data: &TurnbullData, progress_bar: ProgressBar, max_it // -------- // ICM step - let (p_new, s_new, ll_model_new) = do_icm_step(data, &p, &s, ll_tolerance, &constraint, ll_model_after_em); + let ll_model_new; + + if constraint.is_none() { + (p, s, ll_model_new) = do_icm_step(data, &p, &s, ll_tolerance, ll_model_after_em); + } else { + // ICM step is very slow with constraints, so skip it and just do EM + ll_model_new = ll_model_after_em; + } let ll_change = ll_model_new - ll_model; let converged = ll_change <= ll_tolerance; - p = p_new; - s = s_new; ll_model = ll_model_new; // Estimate progress bar according to either the order of magnitude of the ll_change relative to tolerance, or iteration/max_iterations @@ -421,7 +428,7 @@ fn do_em_step(data: &TurnbullData, p: &Vec, s: &Vec, constraint: &Opti return p_new; } -fn do_icm_step(data: &TurnbullData, p: &Vec, s: &Vec, ll_tolerance: f64, constraint: &Option, ll_model: f64) -> (Vec, Vec, f64) { +fn do_icm_step(data: &TurnbullData, p: &Vec, s: &Vec, ll_tolerance: f64, /* constraint: &Option, */ ll_model: f64) -> (Vec, Vec, f64) { // Compute Λ, the cumulative hazard // Since Λ = -inf when survival is 1, and Λ = inf when survival is 0, these are omitted // The entry at lambda[j] corresponds to the survival immediately before time point j + 1 @@ -495,15 +502,16 @@ fn do_icm_step(data: &TurnbullData, p: &Vec, s: &Vec, ll_tolerance: f6 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(); - } + // This is very slow, so support constraints only in the EM step + //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);