diff --git a/doc/source/api-reference/qibo.rst b/doc/source/api-reference/qibo.rst
index 451569479c..c2dc1358a5 100644
--- a/doc/source/api-reference/qibo.rst
+++ b/doc/source/api-reference/qibo.rst
@@ -2269,6 +2269,12 @@ Hellinger fidelity
 .. autofunction:: qibo.quantum_info.hellinger_fidelity
 
 
+Hellinger shot error
+""""""""""""""""""""
+
+.. autofunction:: qibo.quantum_info.hellinger_fidelity
+
+
 Haar integral
 """""""""""""
 
diff --git a/src/qibo/noise_model.py b/src/qibo/noise_model.py
index 44311c2661..b3be7e84b8 100644
--- a/src/qibo/noise_model.py
+++ b/src/qibo/noise_model.py
@@ -1,7 +1,7 @@
 import numpy as np
 
 from qibo import gates, models
-from qibo.quantum_info import hellinger_fidelity
+from qibo.quantum_info.utils import hellinger_fidelity, hellinger_shot_error
 
 
 def noisy_circuit(circuit, params):
@@ -200,28 +200,6 @@ def freq_to_prob(freq):
     return prob
 
 
-def hellinger_shot_error(p, q, nshots):
-    """Hellinger fidelity error caused by using two probability distributions estimated using a finite number of shots.
-    It is calculated propagating the probability error of each state of the system. The complete formula is:
-    :math:`(1 - H^{2}(p, q))/\\sqrt{nshots} * \\sum_{i=1}^{n}(\\sqrt{p_i(1-q_i)}+\\sqrt{q_i(1-p_i)})`
-    where the sum is made all over the possible states and :math:`H(p, q)` is the Hellinger distance.
-
-       Args:
-           p (numpy.ndarray): (discrete) probability distribution :math:`p`.
-           q (numpy.ndarray): (discrete) probability distribution :math:`q`.
-           nshots (int): the number of shots we used to run the circuit to obtain :math:`p` and :math:`q`.
-
-       Returns:
-           (float): The Hellinger fidelity error.
-
-    """
-    hellinger_fid = hellinger_fidelity(p, q)
-    hellinger_fid_e = np.sqrt(hellinger_fid / nshots) * np.sum(
-        np.sqrt(q * (1 - p)) + np.sqrt(p * (1 - q))
-    )
-    return hellinger_fid_e
-
-
 def loss(parameters, *args):
     """The loss function used to be maximized in the fit method of the :class:`qibo.noise_model.CompositeNoiseModel`.
     It is the hellinger fidelity calculated between the probability distribution of the noise model and the experimental target distribution using the :func:`qibo.quantum_info.hellinger_fidelity`.
diff --git a/src/qibo/quantum_info/utils.py b/src/qibo/quantum_info/utils.py
index ec772168e0..a162561fd2 100644
--- a/src/qibo/quantum_info/utils.py
+++ b/src/qibo/quantum_info/utils.py
@@ -258,14 +258,13 @@ def hellinger_fidelity(prob_dist_p, prob_dist_q, validate: bool = False, backend
     .. math::
         (1 - H^{2}(p, q))^{2} \\, ,
 
-    where :math:`H(p, q)` is the Hellinger distance
-    (:func:`qibo.quantum_info.utils.hellinger_distance`).
+    where :math:`H(p, q)` is the :func:`qibo.quantum_info.utils.hellinger_distance`.
 
     Args:
         prob_dist_p (ndarray or list): discrete probability distribution :math:`p`.
         prob_dist_q (ndarray or list): discrete probability distribution :math:`q`.
-        validate (bool, optional): if True, checks if :math:`p` and :math:`q` are proper
-            probability distributions. Default: False.
+        validate (bool, optional): if ``True``, checks if :math:`p` and :math:`q` are proper
+            probability distributions. Defaults to ``False``.
         backend (:class:`qibo.backends.abstract.Backend`, optional): backend to be
             used in the execution. If ``None``, it uses
             :class:`qibo.backends.GlobalBackend`. Defaults to ``None``.
@@ -274,11 +273,62 @@ def hellinger_fidelity(prob_dist_p, prob_dist_q, validate: bool = False, backend
         (float): Hellinger fidelity.
 
     """
+    backend = _check_backend(backend)
+
     distance = hellinger_distance(prob_dist_p, prob_dist_q, validate, backend=backend)
 
     return (1 - distance**2) ** 2
 
 
+def hellinger_shot_error(
+    prob_dist_p, prob_dist_q, nshots: int, validate: bool = False, backend=None
+):
+    """Calculates the Hellinger fidelity error between two discrete probability distributions estimated from finite statistics.
+
+    It is calculated propagating the probability error of each state of the system.
+    The complete formula is:
+
+    .. math::
+        \\frac{1 - H^{2}(p, q)}{\\sqrt{nshots}} \\, \\sum_{k} \\,
+            \\left(\\sqrt{p_{k} \\, (1 - q_{k})} + \\sqrt{q_{k} \\, (1 - p_{k})}\\right)
+
+    where :math:`H(p, q)` is the :func:`qibo.quantum_info.utils.hellinger_distance`,
+    and :math:`1 - H^{2}(p, q)` is the square root of the
+    :func:`qibo.quantum_info.utils.hellinger_fidelity`.
+
+    Args:
+        prob_dist_p (ndarray or list): discrete probability distribution :math:`p`.
+        prob_dist_q (ndarray or list): discrete probability distribution :math:`q`.
+        nshots (int): number of shots we used to run the circuit to obtain :math:`p` and :math:`q`.
+        validate (bool, optional): if ``True``, checks if :math:`p` and :math:`q` are proper
+            probability distributions. Defaults to ``False``.
+        backend (:class:`qibo.backends.abstract.Backend`, optional): backend to be
+            used in the execution. If ``None``, it uses
+            :class:`qibo.backends.GlobalBackend`. Defaults to ``None``.
+
+    Returns:
+        (float): Hellinger fidelity error.
+
+    """
+    backend = _check_backend(backend)
+
+    if isinstance(prob_dist_p, list):
+        prob_dist_p = backend.cast(prob_dist_p, dtype=np.float64)
+
+    if isinstance(prob_dist_q, list):
+        prob_dist_q = backend.cast(prob_dist_q, dtype=np.float64)
+
+    hellinger_error = hellinger_fidelity(
+        prob_dist_p, prob_dist_q, validate=validate, backend=backend
+    )
+    hellinger_error = np.sqrt(hellinger_error / nshots) * np.sum(
+        np.sqrt(prob_dist_q * (1 - prob_dist_p))
+        + np.sqrt(prob_dist_p * (1 - prob_dist_q))
+    )
+
+    return hellinger_error
+
+
 def haar_integral(
     nqubits: int,
     power_t: int,
diff --git a/tests/test_quantum_info_utils.py b/tests/test_quantum_info_utils.py
index 049f800275..0d845cf2f3 100644
--- a/tests/test_quantum_info_utils.py
+++ b/tests/test_quantum_info_utils.py
@@ -4,9 +4,10 @@
 import numpy as np
 import pytest
 
-from qibo import Circuit, matrices
+from qibo import Circuit, gates, matrices
 from qibo.config import PRECISION_TOL
 from qibo.quantum_info.metrics import fidelity
+from qibo.quantum_info.random_ensembles import random_clifford
 from qibo.quantum_info.utils import (
     haar_integral,
     hadamard_transform,
@@ -14,6 +15,7 @@
     hamming_weight,
     hellinger_distance,
     hellinger_fidelity,
+    hellinger_shot_error,
     pqc_integral,
 )
 
@@ -175,6 +177,34 @@ def test_hellinger(backend, validate, kind):
     assert fidelity == (1 - target**2) ** 2
 
 
+@pytest.mark.parametrize("kind", [None, list])
+@pytest.mark.parametrize("validate", [False, True])
+def test_hellinger_shot_error(backend, validate, kind):
+    nqubits, nshots = 5, 1000
+
+    circuit = random_clifford(nqubits, seed=1, backend=backend)
+    circuit.add(gates.M(qubit) for qubit in range(nqubits))
+
+    circuit_2 = random_clifford(nqubits, seed=2, backend=backend)
+    circuit_2.add(gates.M(qubit) for qubit in range(nqubits))
+
+    prob_dist_p = backend.execute_circuit(circuit, nshots=nshots).probabilities()
+    prob_dist_q = backend.execute_circuit(circuit_2, nshots=nshots).probabilities()
+
+    if kind is not None:
+        prob_dist_p = kind(prob_dist_p)
+        prob_dist_q = kind(prob_dist_q)
+
+    hellinger_error = hellinger_shot_error(
+        prob_dist_p, prob_dist_q, nshots, validate=validate, backend=backend
+    )
+    hellinger_fid = hellinger_fidelity(
+        prob_dist_p, prob_dist_q, validate=validate, backend=backend
+    )
+
+    assert 2 * hellinger_error < hellinger_fid
+
+
 def test_haar_integral_errors(backend):
     with pytest.raises(TypeError):
         nqubits, power_t, samples = 0.5, 2, 10