diff --git a/pints/_error_measures.py b/pints/_error_measures.py
index b04d128c1..6ceae641b 100644
--- a/pints/_error_measures.py
+++ b/pints/_error_measures.py
@@ -88,7 +88,7 @@ class MeanSquaredError(ProblemErrorMeasure):
     """
     def __init__(self, problem, weights=None):
         super(MeanSquaredError, self).__init__(problem)
-        self._ninv = 1.0 / np.product(self._values.shape)
+        self._ninv = 1.0 / np.prod(self._values.shape)
 
         if weights is None:
             weights = [1] * self._n_outputs
diff --git a/pints/_log_likelihoods.py b/pints/_log_likelihoods.py
index 3c7792a7b..d0e37a2c9 100644
--- a/pints/_log_likelihoods.py
+++ b/pints/_log_likelihoods.py
@@ -1120,7 +1120,7 @@ def __init__(self, log_likelihood):
         self._log_likelihood = log_likelihood
 
         # Pre-calculate parts
-        self._f = 1.0 / np.product(self._values.shape)
+        self._f = 1.0 / np.prod(self._values.shape)
 
     def __call__(self, x):
         return self._f * self._log_likelihood(x)
diff --git a/pints/_log_priors.py b/pints/_log_priors.py
index 99311fb7c..eb669ab98 100644
--- a/pints/_log_priors.py
+++ b/pints/_log_priors.py
@@ -1245,7 +1245,7 @@ def __init__(self, lower_or_boundaries, upper=None):
         # Use normalised value (1/area) for rectangular boundaries,
         # otherwise just use 1.
         if isinstance(self._boundaries, pints.RectangularBoundaries):
-            self._value = -np.log(np.product(self._boundaries.range()))
+            self._value = -np.log(np.prod(self._boundaries.range()))
         else:
             self._value = 1
 
diff --git a/pints/_mcmc/__init__.py b/pints/_mcmc/__init__.py
index f65140363..f0ec64baf 100644
--- a/pints/_mcmc/__init__.py
+++ b/pints/_mcmc/__init__.py
@@ -90,7 +90,7 @@ def __init__(self, x0, sigma0=None):
             self._sigma0 = np.diag(0.01 * self._sigma0)
         else:
             self._sigma0 = np.array(sigma0, copy=True)
-            if np.product(self._sigma0.shape) == self._n_parameters:
+            if np.prod(self._sigma0.shape) == self._n_parameters:
                 # Convert from 1d array
                 self._sigma0 = self._sigma0.reshape((self._n_parameters,))
                 self._sigma0 = np.diag(self._sigma0)
@@ -192,7 +192,7 @@ def __init__(self, chains, x0, sigma0=None):
             self._sigma0 = np.diag(0.01 * self._sigma0)
         else:
             self._sigma0 = np.array(sigma0, copy=True)
-            if np.product(self._sigma0.shape) == self._n_parameters:
+            if np.prod(self._sigma0.shape) == self._n_parameters:
                 # Convert from 1d array
                 self._sigma0 = self._sigma0.reshape((self._n_parameters,))
                 self._sigma0 = np.diag(self._sigma0)
@@ -352,7 +352,7 @@ def __init__(
                 else:
                     n_parameters = log_pdf[0].n_parameters()
                 # Make sure sigma0 is a (covariance) matrix
-                if np.product(sigma0.shape) == n_parameters:
+                if np.prod(sigma0.shape) == n_parameters:
                     # Convert from 1d array
                     sigma0 = sigma0.reshape((n_parameters,))
                     sigma0 = np.diag(sigma0)
diff --git a/pints/_optimisers/__init__.py b/pints/_optimisers/__init__.py
index e95c31a64..5887051d3 100644
--- a/pints/_optimisers/__init__.py
+++ b/pints/_optimisers/__init__.py
@@ -1125,7 +1125,7 @@ def __init__(self, function, dimension, x, y):
         self.d = dimension
         self.x = x
         self.y = y
-        self.n = 1 / np.product(y.shape)    # Total number of points in data
+        self.n = 1 / np.prod(y.shape)    # Total number of points in data
 
     def n_parameters(self):
         return self.d
diff --git a/pints/tests/test_log_priors.py b/pints/tests/test_log_priors.py
index a77c9fca8..499ac341d 100755
--- a/pints/tests/test_log_priors.py
+++ b/pints/tests/test_log_priors.py
@@ -1081,7 +1081,7 @@ def test_uniform_prior(self):
         self.assertEqual(p([10, 10]), m)
         self.assertEqual(p([5, 20]), m)
 
-        w = -np.log(np.product(upper - lower))
+        w = -np.log(np.prod(upper - lower))
         self.assertEqual(p([1, 2]), w)
         self.assertEqual(p([1, 5]), w)
         self.assertEqual(p([1, 20 - 1e-14]), w)
@@ -1110,7 +1110,7 @@ def test_uniform_prior(self):
         self.assertEqual(p([10, 10]), m)
         self.assertEqual(p([5, 20]), m)
 
-        w = -np.log(np.product(upper - lower))
+        w = -np.log(np.prod(upper - lower))
         self.assertEqual(p([1, 2]), w)
         self.assertEqual(p([1, 5]), w)
         self.assertEqual(p([1, 20 - 1e-14]), w)