diff --git a/mapie/conformity_scores/residual_conformity_scores.py b/mapie/conformity_scores/residual_conformity_scores.py
index f1252fbef..1010aee9f 100644
--- a/mapie/conformity_scores/residual_conformity_scores.py
+++ b/mapie/conformity_scores/residual_conformity_scores.py
@@ -433,7 +433,7 @@ def get_estimation_distribution(
         ``conformity_scores`` can be either the conformity scores or
         the quantile of the conformity scores.
         """
-        r_pred = self._predict_residual_estimator(X)
+        r_pred = self._predict_residual_estimator(X).reshape((-1, 1))
         if not self.prefit:
             return np.add(
                 y_pred,
diff --git a/mapie/regression/regression.py b/mapie/regression/regression.py
index 0524665e8..805b9afd2 100644
--- a/mapie/regression/regression.py
+++ b/mapie/regression/regression.py
@@ -71,8 +71,9 @@ class MapieRegressor(BaseEstimator, RegressorMixin):
         - ``"split"``, does not involve cross-validation but a division
           of the data into training and calibration subsets. The splitter
           used is the following: ``sklearn.model_selection.ShuffleSplit``.
+          ``method`` parameter is set to ``"base"``.
         - ``"prefit"``, assumes that ``estimator`` has been fitted already,
-          and the ``method`` parameter is ignored.
+          and the ``method`` parameter is set to ``"base"``.
           All data provided in the ``fit`` method is then used
           for computing conformity scores only.
           At prediction time, quantiles of these conformity scores are used
@@ -407,13 +408,19 @@ def _check_fit_parameters(
         ------
         ValueError
             If conformity score is FittedResidualNormalizing score and method
-            is neither ``"prefit"`` or ``"split"``
+            is neither ``"prefit"`` or ``"split"``.
+
+        ValueError
+            If ``cv`` is `"prefit"`` or ``"split"`` and ``method`` is not
+            ``"base"``.
         """
         # Checking
         self._check_parameters()
         cv = check_cv(
             self.cv, test_size=self.test_size, random_state=self.random_state
         )
+        if self.cv in ["split", "prefit"] and self.method != "base":
+            self.method = "base"
         estimator = self._check_estimator(self.estimator)
         agg_function = self._check_agg_function(self.agg_function)
         cs_estimator = check_conformity_score(
@@ -422,8 +429,8 @@ def _check_fit_parameters(
         if isinstance(cs_estimator, ConformalResidualFittingScore) and \
            self.cv not in ["split", "prefit"]:
             raise ValueError(
-                "The FittedResisualNormalizingScore can be used only with "
-                "``split`` and ``prefit`` methods"
+                "The ConformalResidualFittingScore can be used only with "
+                "``cv='split'`` and ``cv='prefit'``"
             )
 
         X, y = indexable(X, y)
diff --git a/mapie/tests/test_conformity_scores.py b/mapie/tests/test_conformity_scores.py
index 1fd4481db..9129963a6 100644
--- a/mapie/tests/test_conformity_scores.py
+++ b/mapie/tests/test_conformity_scores.py
@@ -1,3 +1,5 @@
+from typing import Any
+
 import numpy as np
 import pytest
 
@@ -372,3 +374,23 @@ def test_crf_prefit_get_estimation_distribution() -> None:
     crf_conf_score.get_estimation_distribution(
         X_toy, y_pred_list, conf_scores
     )
+
+
+@pytest.mark.parametrize("score", [AbsoluteConformityScore(),
+                                   GammaConformityScore(),
+                                   ConformalResidualFittingScore()])
+@pytest.mark.parametrize("alpha", [[0.3], [0.5, 0.4]])
+def test_intervals_shape_with_every_score(
+    score: ConformityScore,
+    alpha: Any
+) -> None:
+    mapie_reg = MapieRegressor(
+        method="base", cv="split", conformity_score=score
+    )
+    X = np.concatenate((X_toy, X_toy))
+    y = np.concatenate((y_toy, y_toy))
+    mapie_reg = mapie_reg.fit(X, y)
+    y_pred, intervals = mapie_reg.predict(X, alpha=alpha)
+    n_samples = X.shape[0]
+    assert y_pred.shape[0] == n_samples
+    assert intervals.shape == (n_samples, 2, len(alpha))