diff --git a/doc/source/api-reference/qibo.rst b/doc/source/api-reference/qibo.rst
index 7438c3df13..ff8308a5ad 100644
--- a/doc/source/api-reference/qibo.rst
+++ b/doc/source/api-reference/qibo.rst
@@ -1499,6 +1499,12 @@ Random Ensembles
 Functions that can generate random quantum objects.
 
 
+Haar-random :math:`U_{3}`
+"""""""""""""""""""""""""
+
+.. autofunction:: qibo.quantum_info.uniform_sampling_U3
+
+
 Random Gaussian matrix
 """"""""""""""""""""""
 
diff --git a/src/qibo/quantum_info/random_ensembles.py b/src/qibo/quantum_info/random_ensembles.py
index e13c981577..8c9d4f8379 100644
--- a/src/qibo/quantum_info/random_ensembles.py
+++ b/src/qibo/quantum_info/random_ensembles.py
@@ -4,6 +4,7 @@
 from typing import Optional, Union
 
 import numpy as np
+from scipy.stats import rv_continuous
 
 from qibo import Circuit, gates
 from qibo.backends import GlobalBackend, NumpyBackend
@@ -19,6 +20,66 @@
 )
 
 
+class _probability_distribution_sin(rv_continuous):
+    def _pdf(self, theta: float):
+        return 0.5 * np.sin(theta)
+
+    def _cdf(self, theta: float):
+        return np.sin(theta / 2) ** 2
+
+    def _ppf(self, theta: float):
+        return 2 * np.arcsin(np.sqrt(theta))
+
+
+def uniform_sampling_U3(ngates: int, seed=None, backend=None):
+    """Samples parameters for Haar-random :math:`U_{3}`s (:class:`qibo.gates.U3`).
+
+    Args:
+        ngates (int): Total number of :math:`U_{3}`s to be sampled.
+        seed (int or :class:`numpy.random.Generator`, optional): Either a generator of random
+            numbers or a fixed seed to initialize a generator. If ``None``, initializes
+            a generator with a random seed. Default: ``None``.
+        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:
+        (ndarray): array of shape (``ngates``, :math:`3`).
+    """
+    if not isinstance(ngates, int):
+        raise_error(
+            TypeError, f"ngates must be type int, but it is type {type(ngates)}."
+        )
+    elif ngates <= 0:
+        raise_error(ValueError, f"ngates must be non-negative, but it is {ngates}.")
+
+    if (
+        seed is not None
+        and not isinstance(seed, int)
+        and not isinstance(seed, np.random.Generator)
+    ):
+        raise_error(
+            TypeError, "seed must be either type int or numpy.random.Generator."
+        )
+
+    if backend is None:  # pragma: no cover
+        backend = GlobalBackend()
+
+    local_state = (
+        np.random.default_rng(seed) if seed is None or isinstance(seed, int) else seed
+    )
+
+    sampler = _probability_distribution_sin(a=0, b=np.pi, seed=local_state)
+    phases = local_state.random((ngates, 3))
+    phases[:, 0] = sampler.rvs(size=len(phases[:, 0]))
+    phases[:, 1] = phases[:, 1] * 2 * np.pi
+    phases[:, 2] = phases[:, 2] * 2 * np.pi
+
+    phases = backend.cast(phases, dtype=phases.dtype)
+
+    return phases
+
+
 def random_gaussian_matrix(
     dims: int,
     rank: Optional[int] = None,
diff --git a/tests/test_quantum_info_random.py b/tests/test_quantum_info_random.py
index 5f7331dad9..7704970cc6 100644
--- a/tests/test_quantum_info_random.py
+++ b/tests/test_quantum_info_random.py
@@ -5,10 +5,11 @@
 import numpy as np
 import pytest
 
-from qibo import matrices
+from qibo import Circuit, gates, matrices
 from qibo.config import PRECISION_TOL
 from qibo.quantum_info.metrics import purity
 from qibo.quantum_info.random_ensembles import (
+    _probability_distribution_sin,
     random_clifford,
     random_density_matrix,
     random_gaussian_matrix,
@@ -19,10 +20,54 @@
     random_statevector,
     random_stochastic_matrix,
     random_unitary,
+    uniform_sampling_U3,
 )
 
 
-@pytest.mark.parametrize("seed", [None, 10, np.random.Generator(np.random.MT19937(10))])
+@pytest.mark.parametrize("seed", [None, 10, np.random.default_rng(10)])
+def test_uniform_sampling_U3(backend, seed):
+    with pytest.raises(TypeError):
+        uniform_sampling_U3("1", seed=seed, backend=backend)
+    with pytest.raises(ValueError):
+        uniform_sampling_U3(0, seed=seed, backend=backend)
+    with pytest.raises(TypeError):
+        uniform_sampling_U3(2, seed="1", backend=backend)
+
+    X = backend.cast(matrices.X, dtype=matrices.X.dtype)
+    Y = backend.cast(matrices.Y, dtype=matrices.Y.dtype)
+    Z = backend.cast(matrices.Z, dtype=matrices.Z.dtype)
+
+    ngates = int(1e4)
+    phases = uniform_sampling_U3(ngates, seed=seed, backend=backend)
+
+    # expectation values in the 3 directions should be the same
+    expectation_values = []
+    for row in phases:
+        row = [float(phase) for phase in row]
+        circuit = Circuit(1)
+        circuit.add(gates.U3(0, *row))
+        state = backend.execute_circuit(circuit).state()
+
+        expectation_values.append(
+            [
+                np.conj(state) @ X @ state,
+                np.conj(state) @ Y @ state,
+                np.conj(state) @ Z @ state,
+            ]
+        )
+    expectation_values = backend.cast(expectation_values)
+    expectation_values = np.mean(expectation_values, axis=0)
+
+    backend.assert_allclose(expectation_values[0], expectation_values[1], atol=1e-1)
+    backend.assert_allclose(expectation_values[0], expectation_values[2], atol=1e-1)
+
+    # execution for coverage
+    sampler = _probability_distribution_sin(a=0, b=np.pi, seed=seed)
+    sampler.pdf(1)
+    sampler.cdf(1)
+
+
+@pytest.mark.parametrize("seed", [None, 10, np.random.default_rng(10)])
 def test_random_gaussian_matrix(backend, seed):
     with pytest.raises(TypeError):
         dims = np.array([2])