diff --git a/lifelines/fitters/__init__.py b/lifelines/fitters/__init__.py
index 964101ec2..686e44f69 100644
--- a/lifelines/fitters/__init__.py
+++ b/lifelines/fitters/__init__.py
@@ -18,7 +18,7 @@
 import autograd.numpy as anp
 
 from scipy.optimize import minimize, root_scalar
-from scipy.integrate import trapz
+from scipy.integrate import trapezoid
 from scipy import stats
 
 import pandas as pd
@@ -2507,7 +2507,7 @@ def predict_expectation(self, X, conditional_after=None) -> pd.Series:
         warnings.warn("""Approximating the expected value using trapezoid rule.\n""", exceptions.ApproximationWarning)
         subjects = utils._get_index(X)
         v = self.predict_survival_function(X, conditional_after=conditional_after)[subjects]
-        return pd.Series(trapz(v.values.T, v.index), index=subjects).squeeze()
+        return pd.Series(trapezoid(v.values.T, v.index), index=subjects).squeeze()
 
     @property
     def median_survival_time_(self):
diff --git a/lifelines/fitters/aalen_additive_fitter.py b/lifelines/fitters/aalen_additive_fitter.py
index a836db3bd..1f39f2c4c 100644
--- a/lifelines/fitters/aalen_additive_fitter.py
+++ b/lifelines/fitters/aalen_additive_fitter.py
@@ -6,7 +6,7 @@
 import numpy as np
 import pandas as pd
 from numpy.linalg import LinAlgError
-from scipy.integrate import trapz
+from scipy.integrate import trapezoid
 
 from lifelines.fitters import RegressionFitter
 from lifelines.utils.printer import Printer
@@ -396,7 +396,7 @@ def predict_expectation(self, X) -> pd.Series:
         """
         index = _get_index(X)
         t = self._index
-        return pd.Series(trapz(self.predict_survival_function(X)[index].values.T, t), index=index)
+        return pd.Series(trapezoid(self.predict_survival_function(X)[index].values.T, t), index=index)
 
     def _compute_confidence_intervals(self):
         ci = 100 * (1 - self.alpha)
diff --git a/lifelines/fitters/coxph_fitter.py b/lifelines/fitters/coxph_fitter.py
index 8da6658d4..d6972a78e 100644
--- a/lifelines/fitters/coxph_fitter.py
+++ b/lifelines/fitters/coxph_fitter.py
@@ -9,7 +9,7 @@
 from numpy import dot, einsum, log, exp, zeros, arange, multiply, ndarray
 import numpy as np
 from scipy.linalg import solve as spsolve, LinAlgError, norm, inv
-from scipy.integrate import trapz
+from scipy.integrate import trapezoid
 from scipy import stats
 from pandas import DataFrame, Series, Index
 import pandas as pd
@@ -2514,7 +2514,7 @@ def predict_expectation(self, X: DataFrame, conditional_after: Optional[ndarray]
         """
         subjects = utils._get_index(X)
         v = self.predict_survival_function(X, conditional_after=conditional_after)[subjects]
-        return pd.Series(trapz(v.values.T, v.index), index=subjects)
+        return pd.Series(trapezoid(v.values.T, v.index), index=subjects)
 
     def _compute_baseline_hazard(self, partial_hazards: DataFrame, name: Any) -> pd.DataFrame:
         # https://stats.stackexchange.com/questions/46532/cox-baseline-hazard
diff --git a/lifelines/utils/__init__.py b/lifelines/utils/__init__.py
index 7fbf5cbf8..9b1767f31 100644
--- a/lifelines/utils/__init__.py
+++ b/lifelines/utils/__init__.py
@@ -11,7 +11,7 @@
 from numpy import ndarray
 import numpy as np
 
-from scipy.integrate import quad, trapz
+from scipy.integrate import quad, trapezoid
 from scipy.linalg import solve
 from scipy import stats
 
@@ -266,7 +266,7 @@ def _expected_value_of_survival_up_to_t(model_or_survival_function, t: float = n
         )
         sf = model_or_survival_function.loc[:t]
         sf = pd.concat((sf, pd.DataFrame([1], index=[0], columns=sf.columns))).sort_index()
-        return trapz(y=sf.values[:, 0], x=sf.index)
+        return trapezoid(y=sf.values[:, 0], x=sf.index)
     elif isinstance(model_or_survival_function, lifelines.fitters.UnivariateFitter):
         # lifelines model
         model = model_or_survival_function
@@ -313,7 +313,7 @@ def _expected_value_of_survival_squared_up_to_t(
         sf = model_or_survival_function.loc[:t]
         sf = pd.concat((sf, pd.DataFrame([1], index=[0], columns=sf.columns))).sort_index()
         sf_tau = sf * sf.index.values[:, None]
-        return 2 * trapz(y=sf_tau.values[:, 0], x=sf_tau.index)
+        return 2 * trapezoid(y=sf_tau.values[:, 0], x=sf_tau.index)
     elif isinstance(model_or_survival_function, lifelines.fitters.UnivariateFitter):
         # lifelines model
         model = model_or_survival_function