From 0e39402d3dc3656400fd692147031e174baf7c12 Mon Sep 17 00:00:00 2001
From: RunasSudo <runassudo@yingtongli.me>
Date: Sat, 28 Oct 2023 00:49:44 +1100
Subject: [PATCH] turnbull: Compute and display log-likelihood

---
 src/turnbull.rs | 17 +++++++++++++++++
 1 file changed, 17 insertions(+)

diff --git a/src/turnbull.rs b/src/turnbull.rs
index caac5e9..d08e827 100644
--- a/src/turnbull.rs
+++ b/src/turnbull.rs
@@ -81,6 +81,7 @@ pub fn main(args: TurnbullArgs) {
 		OutputFormat::Text => {
 			println!();
 			println!();
+			println!("LL = {:.5}", result.ll_model);
 			
 			let mut summary = Table::new();
 			let format = format::FormatBuilder::new()
@@ -216,6 +217,9 @@ pub fn fit_turnbull(data_times: MatrixXx2<f64>, progress_bar: ProgressBar, max_i
 		survival_prob.push(acc);
 	}
 	
+	// Compute log-likelihood
+	let ll = compute_log_likelihood(&data, &s);
+	
 	// --------------------------------------------------
 	// Compute standard errors for survival probabilities
 	
@@ -263,6 +267,7 @@ pub fn fit_turnbull(data_times: MatrixXx2<f64>, progress_bar: ProgressBar, max_i
 		failure_prob: s,
 		survival_prob: survival_prob,
 		survival_prob_se: survival_prob_se.data.as_vec().clone(),
+		ll_model: ll,
 	};
 }
 
@@ -361,6 +366,17 @@ fn compute_pi(data: &TurnbullData, s: &Vec<f64>, sum_fail_prob: Vec<f64>) -> Vec
 	return pi;
 }
 
+fn compute_log_likelihood(data: &TurnbullData, s: &Vec<f64>) -> f64 {
+	let mut ll = 0.0;
+	
+	for (idx_left, idx_right) in data.data_time_interval_indexes.iter() {
+		let likelihood_ob: f64 = s[*idx_left..(*idx_right + 1)].iter().sum();
+		ll += likelihood_ob.ln();
+	}
+	
+	return ll;
+}
+
 fn compute_hessian(data: &TurnbullData, s: &Vec<f64>) -> DMatrix<f64> {
 	let mut hessian: DMatrix<f64> = DMatrix::zeros(data.num_intervals() - 1, data.num_intervals() - 1);
 	
@@ -415,4 +431,5 @@ pub struct TurnbullResult {
 	pub failure_prob: Vec<f64>,
 	pub survival_prob: Vec<f64>,
 	pub survival_prob_se: Vec<f64>,
+	pub ll_model: f64,
 }