diff --git a/src/turnbull.rs b/src/turnbull.rs
index 7713aaf..ffc2c0c 100644
--- a/src/turnbull.rs
+++ b/src/turnbull.rs
@@ -328,13 +328,17 @@ fn get_sum_fail_prob(data: &TurnbullData, s: &DVector<f64>) -> DVector<f64> {
 }
 
 fn compute_pi(data: &TurnbullData, s: &DVector<f64>, sum_fail_prob: DVector<f64>) -> DVector<f64> {
-	let mut pi: DVector<f64> = DVector::zeros(data.num_intervals());
-	for (i, (idx_left, idx_right)) in data.data_time_interval_indexes.iter().enumerate() {
+	// Faster to repeatedly index Vec than DVector, so first work on Vec then convert to DVector at the end
+	let mut pi: Vec<f64> = vec![0.0; data.num_intervals()];
+	let s_vec = s.data.as_vec();
+	
+	for ((idx_left, idx_right), sum_fail_prob_i) in data.data_time_interval_indexes.iter().zip(sum_fail_prob.iter()) {
 		for j in *idx_left..(*idx_right + 1) {
-			pi[j] += s[j] / sum_fail_prob[i] / data.num_obs() as f64;
+			pi[j] += s_vec[j] / sum_fail_prob_i / data.num_obs() as f64;
 		}
 	}
-	return pi;
+	
+	return DVector::from_vec(pi);
 }
 
 fn compute_hessian(data: &TurnbullData, s: &DVector<f64>) -> DMatrix<f64> {