diff --git a/doc/api.rst b/doc/api.rst
index 5300bb7a..c2d1736c 100644
--- a/doc/api.rst
+++ b/doc/api.rst
@@ -72,7 +72,7 @@ Conformity scores
 
    conformity_scores.AbsoluteConformityScore
    conformity_scores.GammaConformityScore
-
+   conformity_scores.ConformalResidualFittingScore
 
 Resampling
 ==========
diff --git a/doc/index.rst b/doc/index.rst
index b3da3ed9..53038bf4 100644
--- a/doc/index.rst
+++ b/doc/index.rst
@@ -13,6 +13,7 @@
    :caption: REGRESSION
 
    theoretical_description_regression
+   theoretical_description_conformity_scores
    examples_regression/4-tutorials/plot_main-tutorial-regression
    examples_regression/4-tutorials/plot_cqr_tutorial
    examples_regression/4-tutorials/plot_ts-tutorial
diff --git a/doc/theoretical_description_conformity_scores.rst b/doc/theoretical_description_conformity_scores.rst
new file mode 100644
index 00000000..7955db47
--- /dev/null
+++ b/doc/theoretical_description_conformity_scores.rst
@@ -0,0 +1,116 @@
+.. title:: Theoretical Description : contents
+
+.. _theoretical_description_conformity_scores:
+
+=============================================
+Theoretical Description for conformity scores
+=============================================
+
+The :class:`mapie.conformity_scores.ConformityScore` class implements various
+methods to compute conformity scores for regression.
+We give here a brief theoretical description of the scores included in the module.
+Note that it is possible for the user to create any conformal scores that are not 
+already included in MAPIE by inheriting this class.
+
+Before describing the methods, let's briefly present the mathematical setting.
+With conformal predictions we want to transform an heuristic notion of uncertainty
+from a model into a rigorous one, and the first step to do it is to choose a conformal score.
+The only requirement for the score function :math:`s(X, Y) \in \mathbb{R}` is
+that larger scores should encode worse agreement between :math:`X` and :math:`Y`. [1]
+
+There are two types of scores : the symmetric and asymmetric ones.
+The symmetric property defines the way of computing the quantile of the conformity
+scores when calculating the interval's bounds. If a score is symmetrical two
+quantiles will be computed : one on the right side of the distribution
+and the other on the left side.
+
+1. The absolute residual score
+==============================
+
+The absolute residual score (:class:`mapie.conformity_scores.AbsoluteConformityScore`)
+is the simplest and most commonly used conformal score, it translates the error
+of the model : in regression it is called the residual.
+
+.. math:: |Y-\hat{\mu}(X)|
+
+The intervals of prediction's bounds are then computed from the following formula :
+
+.. math:: [\hat{\mu}(X) - q(s), \hat{\mu}(X) + q(s)]
+
+Where :math:`q(s)` is the :math:`(1-\alpha)` quantile of the conformity scores.
+(see :doc:`theoretical_description_regression` for more details).
+
+With this score the intervals of predictions will be constant over the whole dataset.
+This score is by default symmetric (*see above for definition*).
+
+2. The gamma score
+==================
+
+The gamma score [2] (:class:`mapie.conformity_scores.GammaConformityScore`) adds a
+notion of adaptivity with the normalization of the residuals by the predictions.
+
+.. math:: \frac{|Y-\hat{\mu}(X)|}{\hat{\mu}(X)}
+
+It computes adaptive intervals : intervals of different size on each example, with
+the following formula  :
+
+.. math:: [\hat{\mu}(X) * (1 - q(s)), \hat{\mu}(X) * (1 + q(s))]
+
+Where :math:`q(s)` is the :math:`(1-\alpha)` quantile of the conformity scores.
+(see :doc:`theoretical_description_regression` for more details).
+
+This score is by default asymmetric (*see definition above*).
+
+Compared to the absolute residual score, it allows to see regions with smaller intervals
+than others which are interpreted as regions with more certainty than others.
+It is important to note that, this conformity score is inversely proportional to the
+order of magnitude of the predictions. Therefore, the uncertainty is proportional to
+the order of magitude of the predictions, implying that this score should be used
+in use cases where we want greater uncertainty when the prediction is high.
+
+3. The conformal residual fitting score
+=======================================
+
+The conformal residual fitting score [1] (:class:`mapie.conformity_scores.ConformalResidualFittingScore`)
+(CRF) is slightly more complex than the previous scores.
+The normalization of the residual is now done by the predictions of an additional model
+:math:`\hat\sigma` which learns to predict the base model residuals from :math:`X`.
+:math:`\hat\sigma` is trained on :math:`(X, |Y-\hat{\mu}(X)|)` and the formula of the score is:
+
+.. math:: \frac{|Y-\hat{\mu}(X)|}{\hat{\sigma}(X)}
+
+This score provides adaptive intervals : intervals of different sizes in each point
+with the following formula :
+
+.. math:: [\hat{\mu}(X) - q(s) * \hat{\sigma}(X), \hat{\mu}(X) + q(s) * \hat{\sigma}(X)]
+
+Where :math:`q(s)` is the :math:`(1-\alpha)` quantile of the conformity scores.
+(see :doc:`theoretical_description_regression` for more details).
+
+This score is by default symmetric (*see definition above*). Unlike the scores above,
+and due to the additionnal model required this score can only be used with split methods.
+
+Normalisation by the learned residuals from :math:`X` adds to the score a knowledge of
+:math:`X` and its similarity to the other examples in the dataset.
+Compared to the gamma score, the other adaptive score implemented in MAPIE,
+it is not proportional to the uncertainty.
+
+
+Key takeaways
+=============
+
+- The absolute residual score is the basic conformity score and gives constant intervals. It is the one used by default by :class:`mapie.regression.MapieRegressor`.
+- The gamma conformity score adds a notion of adaptivity by giving intervals of different sizes,
+  and is proportional to the uncertainty.
+- The conformal residual fitting score is a conformity score that requires an additional model
+  to learn the residuals of the model from :math:`X`. It gives very adaptive intervals
+  without specific asumptions on the data.
+
+References
+==========
+
+[1] Angelopoulos, A. N., & Bates, S. (2021). A gentle introduction to conformal
+prediction and distribution-free uncertainty quantification. arXiv preprint arXiv:2107.07511.
+[2] Cordier, T., Blot, V., Lacombe, L., Morzadec, T., Capitaine, A. & Brunel, N.. (2023).
+Flexible and Systematic Uncertainty Estimation with Conformal Prediction via the MAPIE library.
+Available from https://proceedings.mlr.press/v204/cordier23a.html.
diff --git a/doc/theoretical_description_regression.rst b/doc/theoretical_description_regression.rst
index d4f37b19..45f86b1a 100644
--- a/doc/theoretical_description_regression.rst
+++ b/doc/theoretical_description_regression.rst
@@ -26,6 +26,9 @@ feature vector :math:`X_{n+1}` such that
 .. math::
     P \{Y_{n+1} \in \hat{C}_{n, \alpha}(X_{n+1}) \} \geq 1 - \alpha
 
+All the methods below are described with the absolute residual conformity score for simplicity
+but other conformity scores are implemented in MAPIE (see :doc:`theoretical_description_conformity_scores`).
+
 1. The "Naive" method
 =====================
 
diff --git a/examples/regression/2-advanced-analysis/plot_conditional_coverage.py b/examples/regression/2-advanced-analysis/plot_conditional_coverage.py
new file mode 100644
index 00000000..17fc06a7
--- /dev/null
+++ b/examples/regression/2-advanced-analysis/plot_conditional_coverage.py
@@ -0,0 +1,336 @@
+"""
+===============================
+Estimating conditional coverage
+===============================
+This example uses :func:`~mapie.regression.MapieRegressor` with conformal
+scores that returns adaptive intervals i.e.
+(:class:`~mapie.conformity_scores.GammaConformityScore` and
+:class:`~mapie.conformity_scores.ConformalResidualFittingScore`) as well as
+:func:`~mapie.regression.MapieQuantileRegressor`.
+The conditional coverage is computed with the three
+functions that allows to estimate the conditional coverage in regression
+:func:`~mapie.metrics.regression_ssc`,
+:func:`~mapie.metrics.regression_ssc_score` and :func:`~mapie.metrics.hsic`.
+"""
+import warnings
+from typing import Tuple, Union
+
+import pandas as pd
+import numpy as np
+import matplotlib.pyplot as plt
+
+from lightgbm import LGBMRegressor
+
+from mapie._typing import NDArray
+from mapie.regression import MapieQuantileRegressor, MapieRegressor
+from mapie.conformity_scores import (GammaConformityScore,
+                                     ConformalResidualFittingScore)
+from mapie.metrics import (regression_coverage_score_v2,
+                           regression_ssc_score,
+                           hsic, regression_ssc)
+from mapie.subsample import Subsample
+
+warnings.filterwarnings("ignore")
+
+random_state = 42
+split_size = 0.20
+alpha = 0.05
+rng = np.random.default_rng(random_state)
+
+
+# Functions for generating our dataset
+def sin_with_controlled_noise(
+    min_x: Union[int, float],
+    max_x: Union[int, float],
+    n_samples: int,
+) -> Tuple[NDArray, NDArray]:
+    """
+    Generate a dataset following sinx except that one interval over two 1 or -1
+    (0.5 probability) is added to X
+
+    Parameters
+    ----------
+    min_x: Union[int, float]
+        The minimum value for X.
+
+    max_x: Union[int, float]
+        The maximum value for X.
+
+    n_samples: int
+        The number of samples wanted in the dataset.
+
+    Returns
+    -------
+    Tuple[NDArray, NDArray]
+        - X: feature data.
+        - y: target data
+    """
+    X = rng.uniform(min_x, max_x, size=(n_samples, 1)).astype(np.float32)
+    y = np.zeros(shape=(n_samples,))
+
+    for i in range(int(max_x) + 1):
+        indexes = np.argwhere(np.greater_equal(i, X) * np.greater(X, i - 1))
+        if i % 2 == 0:
+            for index in indexes:
+                noise = rng.choice([-1, 1])
+                y[index] = np.sin(X[index][0]) + noise * rng.random() + 2
+        else:
+            for index in indexes:
+                y[index] = np.sin(X[index][0]) + rng.random() / 5 + 2
+
+    return X, y
+
+
+# Data generation
+min_x, max_x, n_samples = 0, 10, 3000
+X_train, y_train = sin_with_controlled_noise(min_x, max_x, n_samples)
+X_test, y_test = sin_with_controlled_noise(min_x, max_x,
+                                           int(n_samples * split_size))
+
+# Definition of our base models
+model = LGBMRegressor(random_state=random_state, alpha=0.5)
+model_quant = LGBMRegressor(
+            objective="quantile",
+            alpha=0.5,
+            random_state=random_state
+)
+
+
+# Definition of the experimental set up
+STRATEGIES = {
+    "CV+": {
+        "method": "plus",
+        "cv": 10,
+    },
+    "JK+ab_Gamma": {
+        "method": "plus",
+        "cv": Subsample(n_resamplings=100),
+        "conformity_score": GammaConformityScore()
+    },
+    "CRF": {
+        "cv": "split",
+        "conformity_score": ConformalResidualFittingScore(
+            residual_estimator=LGBMRegressor(
+                alpha=0.5,
+                random_state=random_state),
+            split_size=0.7,
+            random_state=random_state
+        )
+    },
+    "CQR": {
+        "method": "quantile", "cv": "split", "alpha": alpha
+    },
+}
+
+y_pred, intervals, coverage, cond_coverage, coef_corr = {}, {}, {}, {}, {}
+num_bins = 10
+for strategy, params in STRATEGIES.items():
+    # computing predictions
+    if strategy == "CQR":
+        mapie = MapieQuantileRegressor(
+            model_quant,
+            **params
+        )
+        mapie.fit(X_train, y_train, random_state=random_state)
+        y_pred[strategy], intervals[strategy] = mapie.predict(X_test)
+    else:
+        mapie = MapieRegressor(model, **params, random_state=random_state)
+        mapie.fit(X_train, y_train)
+        y_pred[strategy], intervals[strategy] = mapie.predict(
+            X_test, alpha=alpha
+        )
+
+    # computing metrics
+    coverage[strategy] = regression_coverage_score_v2(
+        y_test, intervals[strategy]
+    )
+    cond_coverage[strategy] = regression_ssc_score(
+        y_test, intervals[strategy], num_bins=num_bins
+    )
+    coef_corr[strategy] = hsic(y_test, intervals[strategy])
+
+
+# Visualisation of the estimated conditional coverage
+estimated_cond_cov = pd.DataFrame(
+    columns=["global coverage", "max coverage violation", "hsic"],
+    index=STRATEGIES.keys())
+for m, cov, ssc, coef in zip(
+    STRATEGIES.keys(),
+    coverage.values(),
+    cond_coverage.values(),
+    coef_corr.values()
+):
+    estimated_cond_cov.loc[m] = [
+        round(cov[0], 2), round(ssc[0], 2), round(coef[0], 2)
+    ]
+
+with pd.option_context('display.max_rows', None, 'display.max_columns', None):
+    print(estimated_cond_cov)
+
+##############################################################################
+# We can see here that the global coverage is approximately the same for
+# all methods. What we want to understand is : "Are these methods good
+# adaptive conformal methods ?". For this we have the two metrics
+# :func:`~mapie.metrics.regression_ssc_score` and :func:`~mapie.metrics.hsic`.
+# - SSC (Size Stratified Coverage) is the maximum violation of the coverage :
+#   the intervals are grouped by width and the coverage is computed for each
+#   group. The lower coverage is the maximum coverage violation. An adaptive
+#   method is one where this maximum violation is as close as possible to the
+#   global coverage. If we interpret the result for the four methods here :
+#   CV+ seems to be the better one.
+# - And with the hsic correlation coefficient, we have the
+#   same interpretation : :func:`~mapie.metrics.hsic` computes the correlation
+#   between the coverage indicator and the interval size, a value of 0
+#   translates an independence between the two.
+#
+# We would like to highlight here the misinterpretation that can be made
+# with these metrics. In fact, here CV+ with the absolute residual score
+# calculates constant intervals which, by definition, are not adaptive.
+# Therefore, it is very important to check that the intervals widths are well
+# spread before drawing conclusions (with a plot of the distribution of
+# interval widths or a visualisation of the data for example).
+#
+# In this example, with the hsic correlation coefficient, none of the methods
+# stand out from the others. However, the SSC score for the method using the
+# gamma score is significantly worse than for CQR and CRF, even though their
+# global coverage is similar. CRF and CQR are very close here, with CRF being
+# slightly more conservative.
+
+
+# Visualition of the data and predictions
+def plot_intervals(X, y, y_pred, intervals, title="", ax=None):
+    """
+    Plots the data X, y with associated intervals and predictions points.
+
+    Parameters
+    ----------
+    X: ArrayLike
+        Observed features
+    y: ArrayLike
+        Observed targets
+    y_pred: ArrayLike
+        Predictions
+    intervals: ArrayLike
+        Prediction intervals
+    title: str
+        Title of the plot
+    ax: matplotlib axes
+        An ax can be provided to include this plot in a subplot
+    """
+    if ax is None:
+        fig, ax = plt.subplots(figsize=(6, 5))
+
+    y_pred = y_pred.reshape((X.shape[0], 1))
+    order = np.argsort(X[:, 0])
+
+    # data
+    ax.scatter(X.ravel(), y, color="#1f77b4", alpha=0.3, label="data")
+    # predictions
+    ax.scatter(
+        X.ravel(), y_pred,
+        color="#ff7f0e", marker="+", label="predictions", alpha=0.5
+    )
+    # intervals
+    for i in range(intervals.shape[-1]):
+        ax.fill_between(
+            X[order].ravel(),
+            intervals[:, 0, i][order],
+            intervals[:, 1, i][order],
+            color="#ff7f0e",
+            alpha=0.3
+        )
+
+    ax.set_xlabel("x")
+    ax.set_ylabel("y")
+    ax.set_title(title)
+    ax.legend()
+
+
+def plot_coverage_by_width(y, intervals, num_bins, alpha, title="", ax=None):
+    """
+    PLots a bar diagram of coverages by groups of interval widths.
+
+    Parameters
+    ----------
+    y: ArrayLike
+        Observed targets.
+    intervals: ArrayLike
+        Intervals of prediction
+    num_bins: int
+        Number of groups of interval widths
+    alpha: float
+        The risk level
+    title: str
+        Title of the plot
+    ax: matplotlib axes
+        An ax can be provided to include this plot in a subplot
+    """
+    if ax is None:
+        fig, ax = plt.subplots(figsize=(6, 5))
+
+    ax.bar(
+        np.arange(num_bins),
+        regression_ssc(y, intervals, num_bins=num_bins)[0]
+    )
+    ax.axhline(y=1 - alpha, color='r', linestyle='-')
+    ax.set_title(title)
+    ax.set_xlabel("intervals grouped by size")
+    ax.set_ylabel("coverage")
+    ax.tick_params(
+        axis='x', which='both', bottom=False, top=False, labelbottom=False
+    )
+
+
+max_width = np.max([
+    np.abs(intervals[strategy][:, 0, 0] - intervals[strategy][:, 1, 0])
+    for strategy in STRATEGIES.keys()])
+
+fig_distr, axs_distr = plt.subplots(nrows=2, ncols=2, figsize=(12, 10))
+fig_viz, axs_viz = plt.subplots(nrows=2, ncols=2, figsize=(12, 10))
+fig_hist, axs_hist = plt.subplots(nrows=2, ncols=2, figsize=(12, 10))
+
+for ax_viz, ax_hist, ax_distr, strategy in zip(
+    axs_viz.flat, axs_hist.flat, axs_distr.flat, STRATEGIES.keys()
+):
+    plot_intervals(
+        X_test, y_test, y_pred[strategy], intervals[strategy],
+        title=strategy, ax=ax_viz
+    )
+    plot_coverage_by_width(
+        y_test, intervals[strategy],
+        num_bins=num_bins, alpha=alpha, title=strategy, ax=ax_hist
+    )
+
+    ax_distr.hist(
+        np.abs(intervals[strategy][:, 0, 0] - intervals[strategy][:, 1, 0]),
+        bins=num_bins
+    )
+    ax_distr.set_xlabel("Interval width")
+    ax_distr.set_ylabel("Occurences")
+    ax_distr.set_title(strategy)
+    ax_distr.set_xlim([0, max_width])
+
+fig_viz.suptitle("Predicted points and intervals on test data")
+fig_distr.suptitle("Distribution of intervals widths")
+fig_hist.suptitle("Coverage by bins of intervals grouped by widths (ssc)")
+plt.tight_layout()
+plt.show()
+
+##############################################################################
+# With toy datasets like this, it is easy to compare visually the methods
+# with a plot of the data and predictions.
+# As mentionned above, a histogram of the ditribution of the interval widths is
+# important to accompany the metrics. It is clear from this histogram
+# that CV+ is not adaptive, the metrics presented here should not be used
+# to evaluate its adaptivity. A wider spread of intervals indicates a more
+# adaptive method.
+# Finally, with the plot of coverage by bins of intervals grouped by widths
+# (which is the output of :func:`~mapie.metrics.regression_ssc`), we want
+# the bins to be as constant as possible around the global coverage (here 0.9).
+
+# As the previous metrics show, gamma score does not perform well in terms of
+# size stratified coverage. It either over-covers or under-covers too much.
+# For CRF and CQR, while the first one has several bins with over-coverage,
+# the second one has more under-coverage. These results are confirmed by the
+# visualisation of the data: CQR is better when the data are more spread out,
+# whereas CRF is better with small intervals.
diff --git a/examples/regression/4-tutorials/plot_cqr_tutorial.py b/examples/regression/4-tutorials/plot_cqr_tutorial.py
index d4f31d34..f370fa78 100644
--- a/examples/regression/4-tutorials/plot_cqr_tutorial.py
+++ b/examples/regression/4-tutorials/plot_cqr_tutorial.py
@@ -46,7 +46,7 @@ class :class:`~mapie.subsample.Subsample` (note that the `alpha` parameter is
 from mapie.regression import MapieQuantileRegressor, MapieRegressor
 from mapie.subsample import Subsample
 
-random_state = 23
+random_state = 18
 rng = np.random.default_rng(random_state)
 round_to = 3
 
@@ -135,7 +135,8 @@ class :class:`~mapie.subsample.Subsample` (note that the `alpha` parameter is
     n_jobs=-1,
     n_iter=10,
     cv=KFold(n_splits=5, shuffle=True),
-    verbose=0
+    verbose=0,
+    random_state=random_state
 )
 optim_model.fit(X_train, y_train)
 estimator = optim_model.best_estimator_
@@ -197,14 +198,14 @@ def plot_prediction_intervals(
     axs.errorbar(
         y_test_sorted_[~warnings],
         y_pred_sorted_[~warnings],
-        yerr=error[~warnings],
+        yerr=np.abs(error[~warnings]),
         capsize=5, marker="o", elinewidth=2, linewidth=0,
         label="Inside prediction interval"
         )
     axs.errorbar(
         y_test_sorted_[warnings],
         y_pred_sorted_[warnings],
-        yerr=error[warnings],
+        yerr=np.abs(error[warnings]),
         capsize=5, marker="o", elinewidth=2, linewidth=0, color="red",
         label="Outside prediction interval"
         )
@@ -265,10 +266,14 @@ def plot_prediction_intervals(
 for strategy, params in STRATEGIES.items():
     if strategy == "cqr":
         mapie = MapieQuantileRegressor(estimator, **params)
-        mapie.fit(X_train, y_train, X_calib=X_calib, y_calib=y_calib)
+        mapie.fit(
+            X_train, y_train,
+            X_calib=X_calib, y_calib=y_calib,
+            random_state=random_state
+        )
         y_pred[strategy], y_pis[strategy] = mapie.predict(X_test)
     else:
-        mapie = MapieRegressor(estimator, **params)
+        mapie = MapieRegressor(estimator, **params, random_state=random_state)
         mapie.fit(X_train, y_train)
         y_pred[strategy], y_pis[strategy] = mapie.predict(X_test, alpha=0.2)
     (
diff --git a/examples/regression/4-tutorials/plot_crf_tutorial.py b/examples/regression/4-tutorials/plot_crf_tutorial.py
new file mode 100644
index 00000000..819b6ae1
--- /dev/null
+++ b/examples/regression/4-tutorials/plot_crf_tutorial.py
@@ -0,0 +1,291 @@
+"""
+=======================================================
+Tutorial for conformal residual fitting score (CRF)
+=======================================================
+We will use the sklearn california housing dataset to understand how the
+conformal residual fitting score works and show the multiple ways of using it.
+
+We will explicit the experimental setup below.
+"""
+import warnings
+
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.ticker import FormatStrFormatter
+from numpy.typing import ArrayLike
+
+from sklearn.datasets import fetch_california_housing
+from sklearn.model_selection import train_test_split
+from sklearn.ensemble import RandomForestRegressor
+from sklearn.linear_model import LinearRegression
+
+from mapie.conformity_scores import ConformalResidualFittingScore
+from mapie.metrics import regression_coverage_score_v2, regression_ssc_score
+from mapie.regression import MapieRegressor
+
+warnings.filterwarnings("ignore")
+
+random_state = 23
+rng = np.random.default_rng(random_state)
+round_to = 3
+
+
+##############################################################################
+# 1. Data
+# --------------------------------------------------------------------------
+# The target variable of this dataset is the median house value for the
+# California districts. This dataset is composed of 8 features, including
+# variables such as the age of the house, the median income of the
+# neighborhood, the average number of rooms or bedrooms or even the location in
+# latitude and longitude. In total there are around 20k observations.
+
+
+data = fetch_california_housing()
+X = data.data
+y = data.target
+
+
+##############################################################################
+# Now let's visualize a histogram of the price of the houses.
+
+fig, axs = plt.subplots(1, 1, figsize=(5, 5))
+axs.hist(y, bins=50)
+axs.set_xlabel("Median price of houses")
+axs.set_title("Histogram of house prices")
+axs.xaxis.set_major_formatter(FormatStrFormatter('%.0f' + "k"))
+plt.show()
+
+
+##############################################################################
+# Let's now create the different splits for the dataset, with a training,
+# calibration, residual and test set. Recall that the calibration set is used
+# for calibrating the prediction intervals and the residual set is used to fit
+# the residual estimator used by the
+# :class:`~mapie.conformity_scores.ConformalResidualFittingScore`.
+
+np.array(X)
+np.array(y)
+X_train, X_test, y_train, y_test = train_test_split(
+    X,
+    y,
+    random_state=random_state,
+    test_size=0.02
+)
+X_train, X_calib, y_train, y_calib = train_test_split(
+    X_train,
+    y_train,
+    random_state=random_state
+)
+X_calib_prefit, X_res, y_calib_prefit, y_res = train_test_split(
+    X_calib,
+    y_calib,
+    random_state=random_state,
+    test_size=0.5
+)
+
+
+##############################################################################
+# 2. Models
+# --------------------------------------------------------------------------
+# We will now define 4 different ways of using the CRF score. Remember that the
+# CRF score is only available in the split setup. First, the simplest one
+# with all the default parameters :
+# a :class:`~sklearn.linear_model.LinearRegression` is used for the residual
+# estimator. (Note that to avoid negative values it is trained with the log
+# of the features and the exponential of the predictions are used).
+# It is also possible to use it with ``cv="prefit"`` i.e. with
+# the base model trained beforehand. The third setup that we illustrate here
+# is with the residual model prefitted : we can set the estimator in parameters
+# of the class, not forgetting to specify ``prefit="True"``. Finally, as an
+# example of the exotic parameterisation we can do : we use as a residual
+# estimator a :class:`~sklearn.linear_model.LinearRegression` wrapped to avoid
+# negative values like it is done by default in the class.
+
+class PosEstim(LinearRegression):
+    def __init__(self):
+        super().__init__()
+
+    def fit(self, X, y):
+        super().fit(
+            X, np.log(np.maximum(y, np.full(y.shape, np.float64(1e-8))))
+        )
+        return self
+
+    def predict(self, X):
+        y_pred = super().predict(X)
+        return np.exp(y_pred)
+
+
+base_model = RandomForestRegressor(n_estimators=10, random_state=random_state)
+base_model = base_model.fit(X_train, y_train)
+
+residual_estimator = RandomForestRegressor(
+    n_estimators=20,
+    max_leaf_nodes=70,
+    min_samples_leaf=7,
+    random_state=random_state
+)
+residual_estimator = residual_estimator.fit(
+    X_res, np.abs(np.subtract(y_res, base_model.predict(X_res)))
+)
+wrapped_residual_estimator = PosEstim().fit(
+    X_res, np.abs(np.subtract(y_res, base_model.predict(X_res)))
+)
+# Estimating prediction intervals
+STRATEGIES = {
+    "Default": {
+        "cv": "split",
+        "conformity_score": ConformalResidualFittingScore()
+    },
+    "Base model prefit": {
+        "cv": "prefit",
+        "estimator": base_model,
+        "conformity_score": ConformalResidualFittingScore(
+            split_size=0.5, random_state=random_state
+        )
+    },
+    "Base and residual model prefit": {
+        "cv": "prefit",
+        "estimator": base_model,
+        "conformity_score": ConformalResidualFittingScore(
+            residual_estimator=residual_estimator,
+            random_state=random_state,
+            prefit=True
+        )
+    },
+    "Wrapped residual model": {
+        "cv": "prefit",
+        "estimator": base_model,
+        "conformity_score": ConformalResidualFittingScore(
+            residual_estimator=wrapped_residual_estimator,
+            random_state=random_state,
+            prefit=True
+        )
+    },
+}
+
+y_pred, intervals, coverage, cond_coverage = {}, {}, {}, {}
+num_bins = 10
+alpha = 0.1
+for strategy, params in STRATEGIES.items():
+    mapie = MapieRegressor(**params, random_state=random_state)
+    if mapie.conformity_score.prefit:
+        mapie.fit(X_calib_prefit, y_calib_prefit)
+    else:
+        mapie.fit(X_calib, y_calib)
+    y_pred[strategy], intervals[strategy] = mapie.predict(X_test, alpha=alpha)
+
+    coverage[strategy] = regression_coverage_score_v2(
+        y_test, intervals[strategy]
+    )
+    cond_coverage[strategy] = regression_ssc_score(
+        y_test, intervals[strategy], num_bins=num_bins
+    )
+
+
+def yerr(y_pred, intervals) -> ArrayLike:
+    """
+    Returns the error bars with the point prediction and its interval
+
+    Parameters
+    ----------
+    y_pred: ArrayLike
+        Point predictions.
+    intervals: ArrayLike
+        Predictions intervals.
+
+    Returns
+    -------
+    ArrayLike
+        Error bars.
+    """
+    return np.abs(np.concatenate(
+        [
+            np.expand_dims(y_pred, 0) - intervals[:, 0, 0].T,
+            intervals[:, 1, 0].T - np.expand_dims(y_pred, 0),
+        ],
+        axis=0,
+    ))
+
+
+def plot_predictions(y, y_pred, intervals, coverage, cond_coverage, ax=None):
+    """
+    Plots y_true against y_pred with the associated interval.
+
+    Parameters
+    ----------
+    y: ArrayLike
+       Observed targets
+    y_pred: ArrayLike
+       Predictions
+    intervals: ArrayLike
+       Prediction intervals
+    coverage: float
+       Global coverage
+    cond_coverage: float
+        Maximum violation coverage
+    ax: matplotlib axes
+       An ax can be provided to include this plot in a subplot
+    """
+    if ax is None:
+        fig, ax = plt.subplots(figsize=(6, 5))
+
+    ax.set_xlim([-0.5, 5.5])
+    ax.set_ylim([-0.5, 5.5])
+
+    error = y_pred - intervals[:, 0, 0]
+    warning1 = y_test > y_pred + error
+    warning2 = y_test < y_pred - error
+    warnings = warning1 + warning2
+    ax.errorbar(
+        y[~warnings],
+        y_pred[~warnings],
+        yerr=np.abs(error[~warnings]),
+        color="g",
+        alpha=0.2,
+        linestyle="None",
+        label="Inside prediction interval"
+    )
+    ax.errorbar(
+        y[warnings],
+        y_pred[warnings],
+        yerr=np.abs(error[warnings]),
+        color="r",
+        alpha=0.3,
+        linestyle="None",
+        label="Outside prediction interval"
+    )
+
+    ax.scatter(y, y_pred, s=3, color="black")
+    ax.plot([0, max(max(y), max(y_pred))], [0, max(max(y), max(y_pred))], "-r")
+    ax.set_title(
+        f"{strategy} - coverage={coverage:.0%} " +
+        f"- max violation={cond_coverage:.0%}"
+    )
+    ax.set_xlabel("y true")
+    ax.set_ylabel("y pred")
+    ax.legend()
+    ax.grid()
+
+
+fig, axs = plt.subplots(nrows=2, ncols=2, figsize=(12, 10))
+for ax, strategy in zip(axs.flat, STRATEGIES.keys()):
+    plot_predictions(
+        y_test,
+        y_pred[strategy],
+        intervals[strategy],
+        coverage[strategy][0],
+        cond_coverage[strategy][0],
+        ax=ax
+    )
+
+fig.suptitle(f"Predicted values and intervals of level {alpha}")
+plt.tight_layout()
+plt.show()
+
+##############################################################################
+# The results show that all the setups reach the global coverage guaranteed of
+# 1-alpha.
+# It is interesting to note that the "base model prefit" and the "wrapped
+# residual model" give exactly the same results. And this is because they are
+# the same models : one prefitted and one fitted directly in the class.
diff --git a/mapie/conformity_scores/residual_conformity_scores.py b/mapie/conformity_scores/residual_conformity_scores.py
index 9c1921e3..bfa95fb6 100644
--- a/mapie/conformity_scores/residual_conformity_scores.py
+++ b/mapie/conformity_scores/residual_conformity_scores.py
@@ -146,7 +146,7 @@ def get_estimation_distribution(
 
 class ConformalResidualFittingScore(ConformityScore):
     """
-    ConformalResidualFittingScore (CRF) score.
+    Conformal residual fitting score (CRF).
 
     The signed conformity score = (|y - y_pred|) / r_pred. r_pred being the
     predicted residual (|y - y_pred|) of the base estimator.
diff --git a/mapie/metrics.py b/mapie/metrics.py
index d4a6cf5f..f50699ac 100644
--- a/mapie/metrics.py
+++ b/mapie/metrics.py
@@ -692,7 +692,7 @@ def hsic(
         and if it has negative or null values.
 
     Examples
-    -------
+    --------
     >>> from mapie.metrics import hsic
     >>> import numpy as np
     >>> y_true = np.array([9.5, 10.5, 12.5])