RunasSudo/hpstat
1
Fork 0
Code Issues Pull Requests Releases Activity

turnbull: Fix issue with likelihood ratio-based confidence intervals

Browse Source 
Previous logic resulted in EM-ICM algorithm terminating too early when constrained, resulting in imprecise confidence intervals
This commit is contained in:
RunasSudo 2023-12-24 22:16:40 +11:00
parent 434b432cf7
commit 606c1636e0
Signed by: RunasSudo
GPG Key ID: 7234E476BF21C61A
1 changed files with 15 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

27
src/turnbull.rs
View File

@ -490,18 +490,18 @@ fn do_icm_step(data: &TurnbullData, p: &Vec<f64>, s: &Vec<f64>, ll_tolerance: f6
let likelihood_obs_new = get_likelihood_obs(data, &s_new); let likelihood_obs_new = get_likelihood_obs(data, &s_new);
ll_model_new = likelihood_obs_new.iter().map(|l| l.ln()).sum(); 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 { 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();
}
return (p_new, s_new, ll_model_new); return (p_new, s_new, ll_model_new);
} }
@ -513,7 +513,10 @@ fn do_icm_step(data: &TurnbullData, p: &Vec<f64>, s: &Vec<f64>, ll_tolerance: f6
step_size_exponent += 1; 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"); panic!("ICM fails to increase log-likelihood");
} }
} }