diff --git a/src/qibo/states.py b/src/qibo/states.py
index 4d3f2d575b..6cb0493f8c 100644
--- a/src/qibo/states.py
+++ b/src/qibo/states.py
@@ -362,9 +362,17 @@ def samples(self, binary=True, registers=False):
                     [gate.result.samples() for gate in self.measurements], axis=1
                 )
             else:
-                # generate new samples
-                probs = self.probabilities(qubits)
-                samples = self.backend.sample_shots(probs, self.nshots)
+                if self._frequencies is not None:
+                    # generate samples that respect the existing frequencies
+                    frequencies = self.frequencies(binary=False)
+                    samples = np.concatenate(
+                        [np.repeat(x, f) for x, f in frequencies.items()]
+                    )
+                    np.random.shuffle(samples)
+                else:
+                    # generate new samples
+                    probs = self.probabilities(qubits)
+                    samples = self.backend.sample_shots(probs, self.nshots)
                 samples = self.backend.samples_to_binary(samples, len(qubits))
                 if self.measurement_gate.has_bitflip_noise():
                     p0, p1 = self.measurement_gate.bitflip_map
diff --git a/tests/test_measurements_probabilistic.py b/tests/test_measurements_probabilistic.py
index 63b5710662..4beb085770 100644
--- a/tests/test_measurements_probabilistic.py
+++ b/tests/test_measurements_probabilistic.py
@@ -1,6 +1,4 @@
 """Test circuit measurements when outcome is probabilistic."""
-import sys
-
 import numpy as np
 import pytest
 
@@ -31,6 +29,23 @@ def test_probabilistic_measurement(backend, accelerators, use_samples):
     assert_result(backend, result, decimal_frequencies=decimal_frequencies)
 
 
+def test_sample_frequency_agreement(backend):
+    # set single-thread to fix the random values generated from the frequency custom op
+    backend.set_threads(1)
+    c = models.Circuit(2)
+    c.add(gates.H(0))
+    c.add(gates.H(1))
+    c.add(gates.M(0, 1))
+    result = backend.execute_circuit(c, nshots=1000)
+
+    backend.set_seed(1234)
+    target_frequencies = result.frequencies(binary=False)
+    samples = result.samples(binary=False)
+    outcomes, counts = np.unique(samples, return_counts=True)
+    frequencies = dict(zip(outcomes, counts))
+    assert frequencies == target_frequencies
+
+
 @pytest.mark.parametrize("use_samples", [True, False])
 def test_unbalanced_probabilistic_measurement(backend, use_samples):
     # set single-thread to fix the random values generated from the frequency custom op