diff --git a/redmapper/mask.py b/redmapper/mask.py
index f998340..dbe04f8 100644
--- a/redmapper/mask.py
+++ b/redmapper/mask.py
@@ -217,7 +217,7 @@ def gen_maskgals(self, maskgalfile):
         normmag = mstar - 2.5 * np.log10(self.config.lval_reference)
         steps = np.arange(10.0, normmag, 0.01)
         f = schechter_pdf(steps, alpha=self.config.calib_lumfunc_alpha, mstar=mstar)
-        n = scipy.integrate.simps(f, steps)
+        n = scipy.integrate.simpson(f, steps)
         maskgals.lum_pdf = schechter_pdf(maskgals.m + mstar, mstar=mstar, alpha=self.config.calib_lumfunc_alpha)
         maskgals.lumwt = maskgals.lum_pdf / n
 
diff --git a/redmapper/zlambda.py b/redmapper/zlambda.py
index e1a2dca..6157913 100644
--- a/redmapper/zlambda.py
+++ b/redmapper/zlambda.py
@@ -542,7 +542,7 @@ def _zlambda_calc_pz(self, z_lambda, wtvals, maxrad, maxmag, slow=False):
         pz = np.exp(ln_lkhd) * self.zredstr.volume_factor[self.zredstr.zindex(pzbins)]
 
         # now normalize
-        n = scipy.integrate.simps(pz, pzbins)
+        n = scipy.integrate.simpson(pz, pzbins)
         pz = pz / n
 
         self.pzbins = pzbins
@@ -723,7 +723,7 @@ def apply_correction(self, lam, zlam, zlam_e, pzbins=None, pzvals=None, noerr=Fa
 
                 pzbins = pdz * np.arange(npzbins) + zlam - pdz * (npzbins - 1)/2. + offset
 
-                n = scipy.integrate.simps(pzvals, pzbins)
+                n = scipy.integrate.simpson(pzvals, pzbins)
                 pzvals /= n
 
                 zlam_e = np.sqrt(zlam_e**2. + extra_err**2.)
diff --git a/redmapper/zred_color.py b/redmapper/zred_color.py
index 49d8071..2298570 100644
--- a/redmapper/zred_color.py
+++ b/redmapper/zred_color.py
@@ -223,7 +223,7 @@ def compute_zred(self, galaxy, no_corrections=False):
             zred_samp = np.zeros(galaxy.zred_samp.size) + zred
         else:
             pz = dist.copy()
-            n = scipy.integrate.simps(pz[gdzbins], self.zredstr.z[gdzbins])
+            n = scipy.integrate.simpson(pz[gdzbins], self.zredstr.z[gdzbins])
             pz /= n
 
             pdf = scipy.interpolate.interp1d(self.zredstr.z[gdzbins], pz[gdzbins], kind='quadratic',