diff --git a/.github/workflows/python-tests.yml b/.github/workflows/python-tests.yml
index a0a064ac..012b281a 100644
--- a/.github/workflows/python-tests.yml
+++ b/.github/workflows/python-tests.yml
@@ -97,7 +97,7 @@ jobs:
         flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
 
     - name: Setup pip dependencies
-      run: pip install pytest-cov coveralls coveragepy-lcov
+      run: pip install pytest-cov coveralls coveragepy-lcov 'coverage<7'
 
     - name: Print conda environment (for debugging)
       run: |
diff --git a/tests/test_coreg/test_biascorr.py b/tests/test_coreg/test_biascorr.py
index 42eea966..fda108a0 100644
--- a/tests/test_coreg/test_biascorr.py
+++ b/tests/test_coreg/test_biascorr.py
@@ -349,7 +349,7 @@ def test_biascorr__bin_and_fit_1d(self, fit_args, fit_func, fit_optimizer, bin_s
             elev_fit_args.update({"niter": 1})
 
         # Run with input parameter, and using only 100 subsamples for speed
-        bcorr.fit(**elev_fit_args, subsample=100, random_state=42)
+        bcorr.fit(**elev_fit_args, subsample=1000, random_state=42)
 
         # Check that variable names are defined during fit
         assert bcorr._meta["bias_var_names"] == ["elevation"]
@@ -521,10 +521,10 @@ def test_deramp__synthetic(self, fit_args, order: int) -> None:
         deramp = biascorr.Deramp(poly_order=order)
         elev_fit_args = fit_args.copy()
         if isinstance(elev_fit_args["to_be_aligned_elev"], gpd.GeoDataFrame):
-            bias_elev = bias_dem.to_pointcloud(data_column_name="z", subsample=30000).ds
+            bias_elev = bias_dem.to_pointcloud(data_column_name="z", subsample=50000).ds
         else:
             bias_elev = bias_dem
-        deramp.fit(elev_fit_args["reference_elev"], to_be_aligned_elev=bias_elev, subsample=20000, random_state=42)
+        deramp.fit(elev_fit_args["reference_elev"], to_be_aligned_elev=bias_elev, subsample=40000, random_state=42)
 
         # Check high-order fit parameters are the same within 10%
         fit_params = deramp._meta["fit_params"]
diff --git a/xdem/fit.py b/xdem/fit.py
index a7ce993e..01b083ac 100644
--- a/xdem/fit.py
+++ b/xdem/fit.py
@@ -535,7 +535,7 @@ def wrapper_cost_sumofsin(p: NDArrayf, x: NDArrayf, y: NDArrayf) -> float:
         # Insert in a scipy bounds object
         scipy_bounds = scipy.optimize.Bounds(lb, ub)
         # First guess for the mean parameters
-        p0 = np.divide(lb + ub, 2).squeeze()
+        p0 = ((lb + ub) / 2).squeeze()
 
         if verbose:
             print("Bounds")