Repository: IBM/differential-privacy-library (language: python)

Issue:
Seedable models
Currently, ` self._rng = secrets.SystemRandom()` is used to set the randomness for DPMechanisms. 

https://github.com/IBM/differential-privacy-library/blob/a889ba0f8d19c77e2b0369451ebc392969fac685/diffprivlib/mechanisms/base.py#L77

I understand that for real models you don't want to set the seed, but is adding a parameter for reproduction purposes an option?


Discussion / hints:
This is something we've been looking to implement, but there are a number of knock-on effects that we need to be wary of. Implementing a seeding on mechanisms would be relatively straight-forward, but special attention is required when looking to seed models or tools, particularly those that initiate more than one mechanism.

Using the same seed to instantiate multiple mechanisms may result in correlated noise (i.e., for a `sum()` over a particular axis), which will break differential privacy. Enabling reproducibility at the cost of differential privacy is undesirable, especially given the potential for unwitting misuse.

At present there is no easy way to communicate between different mechanisms initialised in a given function, so it may require a globally-accessible seed generator. If you have any ideas on this, your thoughts would be greatly appreciated.

[Stackexchange](https://stats.stackexchange.com/questions/233061/best-way-to-seed-n-independent-random-number-generators-from-1-value) has an interesting discussion on this
I'm not sure in python/numpy internals, but I'd think without explicitly setting a seed, a random number is already used as seed. So these effects, such as correlated noise, would also be an issue when not setting a seed, except it not being repeatable. 

Finding a way to truly independently seed each DP-mechanism sounds like an improvement over the current seeding mechanism, but might not be comparable to the current implementation?
Yes, it's also my understanding that RNGs are seeded at random when a seed isn't provided (such as, with the timestamp in milliseconds). In this case we may get correlation by chance. But, seeding many different RNGs with the same seed would result in correlation.

Additionally, it seems that the RNG given by the `secrets` module of Python does not support seeding (despite it being a parameter). It looks like it will be necessary to revert to a different RNG (numpy perhaps) in the instance of seeding, which presents another difficulty of a seeded model not being directly representative of a non-seeded model. Sadly, Numpy doesn't provide any cryptographically-secure RNGs (as desired for DP).

Fix patch:
diff --git a/diffprivlib/__init__.py b/diffprivlib/__init__.py
index fa7b7693..2fac59ee 100644
--- a/diffprivlib/__init__.py
+++ b/diffprivlib/__init__.py
@@ -29,4 +29,4 @@
 from diffprivlib import tools
 from diffprivlib.accountant import BudgetAccountant
 
-__version__ = '0.5.2'
+__version__ = '0.6.dev0'
diff --git a/diffprivlib/mechanisms/base.py b/diffprivlib/mechanisms/base.py
index de8671d0..2554cb58 100644
--- a/diffprivlib/mechanisms/base.py
+++ b/diffprivlib/mechanisms/base.py
@@ -22,7 +22,8 @@
 from copy import copy
 import inspect
 from numbers import Real
-import secrets
+
+from diffprivlib.utils import check_random_state
 
 
 class DPMachine(abc.ABC):
@@ -70,11 +71,16 @@ class DPMechanism(DPMachine, abc.ABC):
         Privacy parameter :math:`\delta` for the mechanism.  Must be in [0, 1].  Cannot be simultaneously zero with
         ``epsilon``.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the mechanism.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     """
-    def __init__(self, *, epsilon, delta):
+    def __init__(self, *, epsilon, delta, random_state=None):
         self.epsilon, self.delta = self._check_epsilon_delta(epsilon, delta)
+        self.random_state = random_state
 
-        self._rng = secrets.SystemRandom()
+        self._rng = check_random_state(random_state, True)
 
     def __repr__(self):
         attrs = inspect.getfullargspec(self.__class__).kwonlyargs
@@ -172,7 +178,7 @@ def _check_all(self, value):
         return True
 
 
-class TruncationAndFoldingMixin:
+class TruncationAndFoldingMixin:  # pylint: disable=too-few-public-methods
     """Mixin for truncating or folding the outputs of a mechanism.  Must be instantiated with a :class:`.DPMechanism`.
 
     Parameters
@@ -225,7 +231,7 @@ def _fold(self, value):
         return value
 
 
-def bernoulli_neg_exp(gamma, rng=None):
+def bernoulli_neg_exp(gamma, random_state=None):
     """Sample from Bernoulli(exp(-gamma)).
 
     Adapted from "The Discrete Gaussian for Differential Privacy", Canonne, Kamath, Steinke, 2020.
@@ -236,8 +242,9 @@ def bernoulli_neg_exp(gamma, rng=None):
     gamma : float
         Parameter to sample from Bernoulli(exp(-gamma)).  Must be non-negative.
 
-    rng : Random number generator, optional
-        Random number generator to use.  If not provided, uses SystemRandom from secrets by default.
+    random_state : int or RandomState, optional
+        Controls the randomness of the mechanism.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
 
     Returns
     -------
@@ -247,8 +254,7 @@ def bernoulli_neg_exp(gamma, rng=None):
     if gamma < 0:
         raise ValueError(f"Gamma must be non-negative, got {gamma}.")
 
-    if rng is None:
-        rng = secrets.SystemRandom()
+    rng = check_random_state(random_state, True)
 
     while gamma > 1:
         gamma -= 1
diff --git a/diffprivlib/mechanisms/binary.py b/diffprivlib/mechanisms/binary.py
index c6c43aeb..b574bcb4 100644
--- a/diffprivlib/mechanisms/binary.py
+++ b/diffprivlib/mechanisms/binary.py
@@ -44,14 +44,18 @@ class Binary(DPMechanism):
     value1 : str
         1st binary label.  Cannot be the same as ``value0``.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the mechanism.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     Notes
     -----
     * The binary attributes, known as `labels`, must be specified as strings.  If non-string labels are required (e.g.
       integer-valued labels), a :class:`.DPTransformer` can be used (e.g. :class:`.IntToString`).
 
     """
-    def __init__(self, *, epsilon, value0, value1):
-        super().__init__(epsilon=epsilon, delta=0.0)
+    def __init__(self, *, epsilon, value0, value1, random_state=None):
+        super().__init__(epsilon=epsilon, delta=0.0, random_state=random_state)
         self.value0, self.value1 = self._check_labels(value0, value1)
 
     @classmethod
diff --git a/diffprivlib/mechanisms/bingham.py b/diffprivlib/mechanisms/bingham.py
index 91f66bb3..81035022 100644
--- a/diffprivlib/mechanisms/bingham.py
+++ b/diffprivlib/mechanisms/bingham.py
@@ -18,6 +18,7 @@
 """
 The Bingham mechanism in differential privacy, for estimating the first eigenvector of a covariance matrix.
 """
+import secrets
 from numbers import Real
 
 import numpy as np
@@ -42,12 +43,17 @@ class Bingham(DPMechanism):
     sensitivity : float, default: 1
         The sensitivity of the mechanism.  Must be in [0, ∞).
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the mechanism.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     """
-    def __init__(self, *, epsilon, sensitivity=1.0):
-        super().__init__(epsilon=epsilon, delta=0)
+    def __init__(self, *, epsilon, sensitivity=1.0, random_state=None):
+        super().__init__(epsilon=epsilon, delta=0, random_state=random_state)
         self.sensitivity = self._check_sensitivity(sensitivity)
 
-        self._rng = np.random.default_rng()
+        if isinstance(self._rng, secrets.SystemRandom):
+            self._rng = np.random.default_rng()
 
     @classmethod
     def _check_epsilon_delta(cls, epsilon, delta):
diff --git a/diffprivlib/mechanisms/exponential.py b/diffprivlib/mechanisms/exponential.py
index 0ca3865c..139114c2 100644
--- a/diffprivlib/mechanisms/exponential.py
+++ b/diffprivlib/mechanisms/exponential.py
@@ -58,9 +58,14 @@ class Exponential(DPMechanism):
     measure : list, optional
         An optional list of measures for each candidate.  If omitted, a uniform measure is used.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the mechanism.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     """
-    def __init__(self, *, epsilon, sensitivity, utility, monotonic=False, candidates=None, measure=None):
-        super().__init__(epsilon=epsilon, delta=0.0)
+    def __init__(self, *, epsilon, sensitivity, utility, monotonic=False, candidates=None, measure=None,
+                 random_state=None):
+        super().__init__(epsilon=epsilon, delta=0.0, random_state=random_state)
         self.sensitivity = self._check_sensitivity(sensitivity)
         self.utility, self.candidates, self.measure = self._check_utility_candidates_measure(utility, candidates,
                                                                                              measure)
@@ -214,10 +219,14 @@ class PermuteAndFlip(Exponential):
     candidates : list, optional
         An optional list of candidate labels.  If omitted, the zero-indexed list [0, 1, ..., n] is used.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the mechanism.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     """
-    def __init__(self, *, epsilon, sensitivity, utility, monotonic=False, candidates=None):
+    def __init__(self, *, epsilon, sensitivity, utility, monotonic=False, candidates=None, random_state=None):
         super().__init__(epsilon=epsilon, sensitivity=sensitivity, utility=utility, monotonic=monotonic,
-                         candidates=candidates, measure=None)
+                         candidates=candidates, measure=None, random_state=random_state)
 
     @copy_docstring(DPMechanism.bias)
     def bias(self, value):
@@ -290,9 +299,13 @@ class ExponentialCategorical(DPMechanism):
         utility), where each `value` is a string and `utility` is a strictly positive float.  A `utility` must be
         specified for every pair of values given in the `utility_list`.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the mechanism.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     """
-    def __init__(self, *, epsilon, utility_list):
-        super().__init__(epsilon=epsilon, delta=0.0)
+    def __init__(self, *, epsilon, utility_list, random_state=None):
+        super().__init__(epsilon=epsilon, delta=0.0, random_state=random_state)
 
         self._balanced_tree = False
         self._utility_values, self._sensitivity, self._domain_values = self._build_utility(utility_list)
@@ -471,6 +484,10 @@ class ExponentialHierarchical(ExponentialCategorical):
         The hierarchy as specified as a nested list of string.  Each string must be a leaf node, and each leaf node
         must lie at the same depth in the hierarchy.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the mechanism.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     Examples
     --------
     Example hierarchies:
@@ -479,10 +496,10 @@ class ExponentialHierarchical(ExponentialCategorical):
     >>> nested_hierarchy = [["A"], ["B"], ["C"], ["D", "E"]]
 
     """
-    def __init__(self, *, epsilon, hierarchy):
+    def __init__(self, *, epsilon, hierarchy, random_state=None):
         self.hierarchy = hierarchy
         utility_list = self._build_utility_list(self._build_hierarchy(hierarchy))
-        super().__init__(epsilon=epsilon, utility_list=utility_list)
+        super().__init__(epsilon=epsilon, utility_list=utility_list, random_state=random_state)
         self._list_hierarchy = None
 
     def _build_hierarchy(self, nested_list, parent_node=None):
diff --git a/diffprivlib/mechanisms/gaussian.py b/diffprivlib/mechanisms/gaussian.py
index ea8008dc..9d094494 100644
--- a/diffprivlib/mechanisms/gaussian.py
+++ b/diffprivlib/mechanisms/gaussian.py
@@ -48,6 +48,10 @@ class Gaussian(DPMechanism):
     sensitivity : float
         The sensitivity of the mechanism.  Must be in [0, ∞).
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the mechanism.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     References
     ----------
     .. [DR14] Dwork, Cynthia, and Aaron Roth. "The algorithmic foundations of differential privacy." Found. Trends
@@ -57,8 +61,8 @@ class Gaussian(DPMechanism):
         arXiv:2107.10138 (2021).
 
     """
-    def __init__(self, *, epsilon, delta, sensitivity):
-        super().__init__(epsilon=epsilon, delta=delta)
+    def __init__(self, *, epsilon, delta, sensitivity, random_state=None):
+        super().__init__(epsilon=epsilon, delta=delta, random_state=random_state)
         self.sensitivity = self._check_sensitivity(sensitivity)
         self._scale = np.sqrt(2 * np.log(1.25 / self.delta)) * self.sensitivity / self.epsilon
 
@@ -105,7 +109,10 @@ def variance(self, value):
     def randomise(self, value):
         self._check_all(value)
 
-        standard_normal = (self._rng.normalvariate(0, 1) + self._rng.normalvariate(0, 1)) / np.sqrt(2)
+        try:
+            standard_normal = (self._rng.normalvariate(0, 1) + self._rng.normalvariate(0, 1)) / np.sqrt(2)
+        except AttributeError:  # random_state is a np.random.RandomState
+            standard_normal = (self._rng.standard_normal() + self._rng.standard_normal()) / np.sqrt(2)
 
         return value + standard_normal * self._scale
 
@@ -129,9 +136,13 @@ class GaussianAnalytic(Gaussian):
     sensitivity : float
         The sensitivity of the mechanism.  Must be in [0, ∞).
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the mechanism.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     """
-    def __init__(self, *, epsilon, delta, sensitivity):
-        super().__init__(epsilon=epsilon, delta=delta, sensitivity=sensitivity)
+    def __init__(self, *, epsilon, delta, sensitivity, random_state=None):
+        super().__init__(epsilon=epsilon, delta=delta, sensitivity=sensitivity, random_state=random_state)
         self._scale = self._find_scale()
 
     @classmethod
@@ -139,7 +150,7 @@ def _check_epsilon_delta(cls, epsilon, delta):
         if epsilon == 0 or delta == 0:
             raise ValueError("Neither Epsilon nor Delta can be zero")
 
-        return super(Gaussian, cls)._check_epsilon_delta(epsilon, delta)
+        return DPMechanism._check_epsilon_delta(epsilon, delta)  # pylint: disable=protected-access
 
     def _check_all(self, value):
         super()._check_all(value)
@@ -217,9 +228,13 @@ class GaussianDiscrete(DPMechanism):
     sensitivity : int, default: 1
         The sensitivity of the mechanism.  Must be in [0, ∞).
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the mechanism.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     """
-    def __init__(self, *, epsilon, delta, sensitivity=1):
-        super().__init__(epsilon=epsilon, delta=delta)
+    def __init__(self, *, epsilon, delta, sensitivity=1, random_state=None):
+        super().__init__(epsilon=epsilon, delta=delta, random_state=random_state)
         self.sensitivity = self._check_sensitivity(sensitivity)
         self._scale = self._find_scale()
 
@@ -272,7 +287,7 @@ def randomise(self, value):
             while bernoulli_neg_exp(tau, self._rng):
                 geom_x += 1
 
-            bern_b = np.random.binomial(1, 0.5)
+            bern_b = self._rng.random() < 0.5
             if bern_b and not geom_x:
                 continue
 
diff --git a/diffprivlib/mechanisms/geometric.py b/diffprivlib/mechanisms/geometric.py
index 2dbed02e..27e87c60 100644
--- a/diffprivlib/mechanisms/geometric.py
+++ b/diffprivlib/mechanisms/geometric.py
@@ -41,9 +41,13 @@ class Geometric(DPMechanism):
     sensitivity : float, default: 1
         The sensitivity of the mechanism.  Must be in [0, ∞).
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the mechanism.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     """
-    def __init__(self, *, epsilon, sensitivity=1):
-        super().__init__(epsilon=epsilon, delta=0.0)
+    def __init__(self, *, epsilon, sensitivity=1, random_state=None):
+        super().__init__(epsilon=epsilon, delta=0.0, random_state=random_state)
         self.sensitivity = self._check_sensitivity(sensitivity)
         self._scale = - self.epsilon / self.sensitivity if self.sensitivity > 0 else - float("inf")
 
@@ -128,9 +132,13 @@ class GeometricTruncated(Geometric, TruncationAndFoldingMixin):
     upper : int
         The upper bound of the mechanism.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the mechanism.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     """
-    def __init__(self, *, epsilon, sensitivity=1, lower, upper):
-        super().__init__(epsilon=epsilon, sensitivity=sensitivity)
+    def __init__(self, *, epsilon, sensitivity=1, lower, upper, random_state=None):
+        super().__init__(epsilon=epsilon, sensitivity=sensitivity, random_state=random_state)
         TruncationAndFoldingMixin.__init__(self, lower=lower, upper=upper)
 
     @classmethod
@@ -184,9 +192,13 @@ class GeometricFolded(Geometric, TruncationAndFoldingMixin):
     upper : int or float
         The upper bound of the mechanism.  Must be integer or half-integer -valued.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the mechanism.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     """
-    def __init__(self, *, epsilon, sensitivity=1, lower, upper):
-        super().__init__(epsilon=epsilon, sensitivity=sensitivity)
+    def __init__(self, *, epsilon, sensitivity=1, lower, upper, random_state=None):
+        super().__init__(epsilon=epsilon, sensitivity=sensitivity, random_state=random_state)
         TruncationAndFoldingMixin.__init__(self, lower=lower, upper=upper)
 
     @classmethod
diff --git a/diffprivlib/mechanisms/laplace.py b/diffprivlib/mechanisms/laplace.py
index a7b00461..9d6f9fa9 100644
--- a/diffprivlib/mechanisms/laplace.py
+++ b/diffprivlib/mechanisms/laplace.py
@@ -48,6 +48,10 @@ class Laplace(DPMechanism):
     sensitivity : float
         The sensitivity of the mechanism.  Must be in [0, ∞).
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the mechanism.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     References
     ----------
     .. [DMNS16] Dwork, Cynthia, Frank McSherry, Kobbi Nissim, and Adam Smith. "Calibrating noise to sensitivity in
@@ -60,8 +64,8 @@ class Laplace(DPMechanism):
         arXiv:2107.10138 (2021).
 
     """
-    def __init__(self, *, epsilon, delta=0.0, sensitivity):
-        super().__init__(epsilon=epsilon, delta=delta)
+    def __init__(self, *, epsilon, delta=0.0, sensitivity, random_state=None):
+        super().__init__(epsilon=epsilon, delta=delta, random_state=random_state)
         self.sensitivity = self._check_sensitivity(sensitivity)
         self._scale = None
 
@@ -168,9 +172,13 @@ class LaplaceTruncated(Laplace, TruncationAndFoldingMixin):
     upper : float
         The upper bound of the mechanism.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the mechanism.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     """
-    def __init__(self, *, epsilon, delta=0.0, sensitivity, lower, upper):
-        super().__init__(epsilon=epsilon, delta=delta, sensitivity=sensitivity)
+    def __init__(self, *, epsilon, delta=0.0, sensitivity, lower, upper, random_state=None):
+        super().__init__(epsilon=epsilon, delta=delta, sensitivity=sensitivity, random_state=random_state)
         TruncationAndFoldingMixin.__init__(self, lower=lower, upper=upper)
 
     @copy_docstring(Laplace.bias)
@@ -233,9 +241,13 @@ class LaplaceFolded(Laplace, TruncationAndFoldingMixin):
     upper : float
         The upper bound of the mechanism.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the mechanism.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     """
-    def __init__(self, *, epsilon, delta=0.0, sensitivity, lower, upper):
-        super().__init__(epsilon=epsilon, delta=delta, sensitivity=sensitivity)
+    def __init__(self, *, epsilon, delta=0.0, sensitivity, lower, upper, random_state=None):
+        super().__init__(epsilon=epsilon, delta=delta, sensitivity=sensitivity, random_state=random_state)
         TruncationAndFoldingMixin.__init__(self, lower=lower, upper=upper)
 
     @copy_docstring(Laplace.bias)
@@ -290,16 +302,16 @@ class LaplaceBoundedDomain(LaplaceTruncated):
     upper : float
         The upper bound of the mechanism.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the mechanism.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     References
     ----------
     .. [HABM20] Holohan, Naoise, Spiros Antonatos, Stefano Braghin, and Pól Mac Aonghusa. "The Bounded Laplace Mechanism
         in Differential Privacy." Journal of Privacy and Confidentiality 10, no. 1 (2020).
 
     """
-    def __init__(self, *, epsilon, delta=0.0, sensitivity, lower, upper):
-        super().__init__(epsilon=epsilon, delta=delta, sensitivity=sensitivity, lower=lower, upper=upper)
-        self._rng = np.random.default_rng()
-
     def _find_scale(self):
         eps = self.epsilon
         delta = self.delta
@@ -396,8 +408,12 @@ def randomise(self, value):
         samples = 1
 
         while True:
-            noisy = value + self._scale * self._laplace_sampler(self._rng.random(samples), self._rng.random(samples),
-                                                                self._rng.random(samples), self._rng.random(samples))
+            try:
+                unif = self._rng.random(4 * samples)
+            except TypeError:  # rng is secrets.SystemRandom
+                unif = [self._rng.random() for _ in range(4 * samples)]
+            noisy = value + self._scale * self._laplace_sampler(*np.array(unif).reshape(4, -1))
+
             if ((noisy >= self.lower) & (noisy <= self.upper)).any():
                 idx = np.argmax((noisy >= self.lower) & (noisy <= self.upper))
                 return noisy[idx]
@@ -424,16 +440,19 @@ class LaplaceBoundedNoise(Laplace):
     sensitivity : float
         The sensitivity of the mechanism.  Must be in [0, ∞).
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the mechanism.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     References
     ----------
     .. [GDGK18] Geng, Quan, Wei Ding, Ruiqi Guo, and Sanjiv Kumar. "Truncated Laplacian Mechanism for Approximate
         Differential Privacy." arXiv preprint arXiv:1810.00877v1 (2018).
 
     """
-    def __init__(self, *, epsilon, delta, sensitivity):
-        super().__init__(epsilon=epsilon, delta=delta, sensitivity=sensitivity)
+    def __init__(self, *, epsilon, delta, sensitivity, random_state=None):
+        super().__init__(epsilon=epsilon, delta=delta, sensitivity=sensitivity, random_state=random_state)
         self._noise_bound = None
-        self._rng = np.random.default_rng()
 
     @classmethod
     def _check_epsilon_delta(cls, epsilon, delta):
@@ -443,7 +462,7 @@ def _check_epsilon_delta(cls, epsilon, delta):
         if isinstance(delta, Real) and not 0 < delta < 0.5:
             raise ValueError("Delta must be strictly in the interval (0,0.5). For zero delta, use :class:`.Laplace`.")
 
-        return super(Laplace, cls)._check_epsilon_delta(epsilon, delta)
+        return super()._check_epsilon_delta(epsilon, delta)
 
     @copy_docstring(Laplace.bias)
     def bias(self, value):
@@ -468,8 +487,12 @@ def randomise(self, value):
         samples = 1
 
         while True:
-            noisy = self._scale * self._laplace_sampler(self._rng.random(samples), self._rng.random(samples),
-                                                        self._rng.random(samples), self._rng.random(samples))
+            try:
+                unif = self._rng.random(4 * samples)
+            except TypeError:  # rng is secrets.SystemRandom
+                unif = [self._rng.random() for _ in range(4 * samples)]
+            noisy = self._scale * self._laplace_sampler(*np.array(unif).reshape(4, -1))
+
             if ((noisy >= - self._noise_bound) & (noisy <= self._noise_bound)).any():
                 idx = np.argmax((noisy >= - self._noise_bound) & (noisy <= self._noise_bound))
                 return value + noisy[idx]
diff --git a/diffprivlib/mechanisms/snapping.py b/diffprivlib/mechanisms/snapping.py
index a86c0c7d..393f3e75 100644
--- a/diffprivlib/mechanisms/snapping.py
+++ b/diffprivlib/mechanisms/snapping.py
@@ -6,7 +6,7 @@
 
 import numpy as np
 try:
-    from crlibm import log_rn
+    from crlibm import log_rn  # pylint: disable=no-name-in-module
 except ModuleNotFoundError:
     log_rn = np.log
 
@@ -40,14 +40,19 @@ class Snapping(LaplaceTruncated):
     upper : float
         The upper bound of the mechanism.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the mechanism.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     References
     ----------
     .. [Mir12] Mironov, Ilya. "On significance of the least significant bits for differential privacy." Proceedings of
      the 2012 ACM conference on Computer and communications security (2012).
 
     """
-    def __init__(self, *, epsilon, sensitivity, lower, upper):
-        super().__init__(epsilon=epsilon, sensitivity=sensitivity, delta=0.0, lower=lower, upper=upper)
+    def __init__(self, *, epsilon, sensitivity, lower, upper, random_state=None):
+        super().__init__(epsilon=epsilon, sensitivity=sensitivity, delta=0.0, lower=lower, upper=upper,
+                         random_state=random_state)
         self._bound = self._scale_bound()
 
     @classmethod
@@ -186,14 +191,20 @@ def _uniform_sampler(self):
 
         """
         mantissa_size = np.finfo(float).nmant
-        mantissa = 1 << mantissa_size | self._rng.getrandbits(mantissa_size)
+        mantissa = 1 << mantissa_size | self._getrandbits(mantissa_size)
         exponent = -(mantissa_size + 1)
         x = 0
         while not x:
-            x = self._rng.getrandbits(32)
+            x = self._getrandbits(32)
             exponent += x.bit_length() - 32
         return np.ldexp(mantissa, exponent)
 
+    def _getrandbits(self, bits):
+        try:
+            return self._rng.getrandbits(bits)
+        except AttributeError:
+            return self._rng.randint(0, 2 ** bits)
+
     @staticmethod
     def _laplace_sampler(unif_bit, unif):
         r"""
@@ -232,8 +243,8 @@ def randomise(self, value):
         value_scaled_offset = self._scale_and_offset_value(value)
         value_clamped = self._truncate(value_scaled_offset)
 
-        scale = 1.0 / self.effective_epsilon()  # everything is scaled to sensitivity 1
+        scale = 1.0 / self.effective_epsilon()  # everything is already scaled to sensitivity 1
         lambda_ = self._get_nearest_power_of_2(scale)
-        laplace = scale * self._laplace_sampler(self._rng.getrandbits(1), self._uniform_sampler())
+        laplace = scale * self._laplace_sampler(self._getrandbits(1), self._uniform_sampler())
         value_rounded = self._round_to_nearest_power_of_2(value_clamped + laplace, lambda_)
         return self._reverse_scale_and_offset_value(self._truncate(value_rounded))
diff --git a/diffprivlib/mechanisms/staircase.py b/diffprivlib/mechanisms/staircase.py
index a18d541f..be294293 100644
--- a/diffprivlib/mechanisms/staircase.py
+++ b/diffprivlib/mechanisms/staircase.py
@@ -18,6 +18,7 @@
 """
 The staircase mechanism in differential privacy.
 """
+import secrets
 from numbers import Real
 
 import numpy as np
@@ -45,12 +46,17 @@ class Staircase(Laplace):
     gamma : float, default: 1 / (1 + exp(epsilon/2))
         Value of the tuning parameter gamma for the mechanism.  Must be in [0, 1].
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the mechanism.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     """
-    def __init__(self, *, epsilon, sensitivity, gamma=None):
-        super().__init__(epsilon=epsilon, delta=0, sensitivity=sensitivity)
+    def __init__(self, *, epsilon, sensitivity, gamma=None, random_state=None):
+        super().__init__(epsilon=epsilon, delta=0, sensitivity=sensitivity, random_state=random_state)
         self.gamma = self._check_gamma(gamma, epsilon=self.epsilon)
 
-        self._rng = np.random.default_rng()
+        if isinstance(self._rng, secrets.SystemRandom):
+            self._rng = np.random.default_rng()
 
     @classmethod
     def _check_gamma(cls, gamma, epsilon=None):
diff --git a/diffprivlib/mechanisms/uniform.py b/diffprivlib/mechanisms/uniform.py
index d393a1ca..f68cd1dc 100644
--- a/diffprivlib/mechanisms/uniform.py
+++ b/diffprivlib/mechanisms/uniform.py
@@ -40,9 +40,13 @@ class Uniform(DPMechanism):
     sensitivity : float
         The sensitivity of the mechanism.  Must be in [0, ∞).
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the mechanism.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     """
-    def __init__(self, *, delta, sensitivity):
-        super().__init__(epsilon=0.0, delta=delta)
+    def __init__(self, *, delta, sensitivity, random_state=None):
+        super().__init__(epsilon=0.0, delta=delta, random_state=random_state)
         self.sensitivity = self._check_sensitivity(sensitivity)
 
     @classmethod
diff --git a/diffprivlib/mechanisms/vector.py b/diffprivlib/mechanisms/vector.py
index b19db814..d59847aa 100644
--- a/diffprivlib/mechanisms/vector.py
+++ b/diffprivlib/mechanisms/vector.py
@@ -52,16 +52,19 @@ class Vector(DPMechanism):
     alpha : float, default: 0.01
         Regularisation parameter.  Must be in (0, ∞).
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the mechanism.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     """
-    def __init__(self, *, epsilon, function_sensitivity, data_sensitivity=1.0, dimension, alpha=0.01):
-        super().__init__(epsilon=epsilon, delta=0.0)
+    def __init__(self, *, epsilon, function_sensitivity, data_sensitivity=1.0, dimension, alpha=0.01,
+                 random_state=None):
+        super().__init__(epsilon=epsilon, delta=0.0, random_state=random_state)
         self.function_sensitivity, self.data_sensitivity = self._check_sensitivity(function_sensitivity,
                                                                                    data_sensitivity)
         self.dimension = self._check_dimension(dimension)
         self.alpha = self._check_alpha(alpha)
 
-        self._rng = np.random.default_rng()
-
     @classmethod
     def _check_epsilon_delta(cls, epsilon, delta):
         if not delta == 0:
@@ -147,10 +150,15 @@ def randomise(self, value):
 
         scale = self.data_sensitivity * 2 / epsilon_p
 
-        normed_noisy_vector = self._rng.standard_normal((self.dimension, 4)).sum(axis=1) / 2
-        norm = np.linalg.norm(normed_noisy_vector, 2)
-        noisy_norm = self._rng.gamma(self.dimension / 4, scale, 4).sum()
+        try:
+            normed_noisy_vector = self._rng.standard_normal((self.dimension, 4)).sum(axis=1) / 2
+            noisy_norm = self._rng.gamma(self.dimension / 4, scale, 4).sum()
+        except AttributeError:  # rng is secrets.SystemRandom
+            normed_noisy_vector = np.reshape([self._rng.normalvariate(0, 1) for _ in range(self.dimension * 4)],
+                                             (-1, 4)).sum(axis=1) / 2
+            noisy_norm = sum(self._rng.gammavariate(self.dimension / 4, scale) for _ in range(4)) if scale > 0 else 0.0
 
+        norm = np.linalg.norm(normed_noisy_vector, 2)
         normed_noisy_vector = normed_noisy_vector / norm * noisy_norm
 
         def output_func(*args):
diff --git a/diffprivlib/models/forest.py b/diffprivlib/models/forest.py
index 11e18ede..ccf7ff73 100644
--- a/diffprivlib/models/forest.py
+++ b/diffprivlib/models/forest.py
@@ -29,7 +29,7 @@
 from sklearn.tree import DecisionTreeClassifier as skDecisionTreeClassifier
 
 from diffprivlib.accountant import BudgetAccountant
-from diffprivlib.utils import PrivacyLeakWarning
+from diffprivlib.utils import PrivacyLeakWarning, check_random_state
 from diffprivlib.mechanisms import PermuteAndFlip
 from diffprivlib.validation import DiffprivlibMixin
 
@@ -66,6 +66,11 @@ class RandomForestClassifier(skRandomForestClassifier, DiffprivlibMixin):  # pyl
     verbose : int, default: 0
         Set to any positive number for verbosity.
 
+    random_state : int or RandomState, optional
+        Controls both the randomness of the shuffling of the samples used when building trees (if ``shuffle=True``) and
+        training of the differentially-private :class:`.DecisionTreeClassifier` to construct the forest.  To obtain a
+        deterministic behaviour during randomisation, ``random_state`` has to be fixed to an integer.
+
     accountant : BudgetAccountant, optional
         Accountant to keep track of privacy budget.
 
@@ -124,13 +129,13 @@ class RandomForestClassifier(skRandomForestClassifier, DiffprivlibMixin):  # pyl
     """
 
     def __init__(self, n_estimators=10, *, epsilon=1.0, bounds=None, classes=None, n_jobs=1, verbose=0, accountant=None,
-                 max_depth=5, warm_start=False, shuffle=False, **unused_args):
+                 random_state=None, max_depth=5, warm_start=False, shuffle=False, **unused_args):
         super().__init__(
             n_estimators=n_estimators,
             criterion=None,
             bootstrap=False,
             n_jobs=n_jobs,
-            random_state=None,
+            random_state=random_state,
             verbose=verbose,
             warm_start=warm_start)
         self.epsilon = epsilon
@@ -209,7 +214,7 @@ def fit(self, X, y, sample_weight=None):
         # Check parameters
         self._validate_estimator()
 
-        # random_state = check_random_state(self.random_state)
+        random_state = check_random_state(self.random_state)
 
         if not self.warm_start or not hasattr(self, "estimators_"):
             # Free allocated memory, if any
@@ -227,16 +232,16 @@ def fit(self, X, y, sample_weight=None):
         if self.warm_start and len(self.estimators_) > 0:
             # We draw from the random state to get the random state we
             # would have got if we hadn't used a warm_start.
-            np.random.randint(MAX_INT, size=len(self.estimators_))
+            random_state.randint(MAX_INT, size=len(self.estimators_))
 
         trees = [
-            self._make_estimator(append=False, random_state=None)
+            self._make_estimator(append=False, random_state=random_state)
             for _ in range(n_more_estimators)
         ]
 
         # Split samples between trees as evenly as possible (randomly if shuffle==True)
         n_samples = X.shape[0]
-        tree_idxs = np.random.permutation(n_samples) if self.shuffle else np.arange(n_samples)
+        tree_idxs = random_state.permutation(n_samples) if self.shuffle else np.arange(n_samples)
         tree_idxs = (tree_idxs // (n_samples / n_more_estimators)).astype(int)
 
         # Parallel loop: we prefer the threading backend as the Cython code
@@ -305,6 +310,12 @@ class DecisionTreeClassifier(skDecisionTreeClassifier, DiffprivlibMixin):
         Array of class labels. If not provided, the classes will be read from the data when ``.fit()`` is first called,
         resulting in a :class:`.PrivacyLeakWarning`.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the estimator.  At each split, the feature to split on is chosen randomly, as is the
+        threshold at which to split.  The classification label at each leaf is then randomised, subject to differential
+        privacy constraints. To obtain a deterministic behaviour during randomisation, ``random_state`` has to be fixed
+        to an integer.
+
     accountant : BudgetAccountant, optional
         Accountant to keep track of privacy budget.
 
@@ -321,7 +332,8 @@ class DecisionTreeClassifier(skDecisionTreeClassifier, DiffprivlibMixin):
 
     """
 
-    def __init__(self, max_depth=5, *, epsilon=1, bounds=None, classes=None, accountant=None, **unused_args):
+    def __init__(self, max_depth=5, *, epsilon=1, bounds=None, classes=None, random_state=None, accountant=None,
+                 **unused_args):
         # TODO: Remove try...except when sklearn v1.0 is min-requirement
         try:
             super().__init__(  # pylint: disable=unexpected-keyword-arg
@@ -332,7 +344,7 @@ def __init__(self, max_depth=5, *, epsilon=1, bounds=None, classes=None, account
                 min_samples_leaf=None,
                 min_weight_fraction_leaf=None,
                 max_features=None,
-                random_state=None,
+                random_state=random_state,
                 max_leaf_nodes=None,
                 min_impurity_decrease=None,
                 min_impurity_split=None
@@ -346,7 +358,7 @@ def __init__(self, max_depth=5, *, epsilon=1, bounds=None, classes=None, account
                 min_samples_leaf=None,
                 min_weight_fraction_leaf=None,
                 max_features=None,
-                random_state=None,
+                random_state=random_state,
                 max_leaf_nodes=None,
                 min_impurity_decrease=None
             )
@@ -379,6 +391,8 @@ def fit(self, X, y, sample_weight=None, check_input=True):
         self : DecisionTreeClassifier
             Fitted estimator.
         """
+        random_state = check_random_state(self.random_state)
+
         self.accountant.check(self.epsilon, 0)
 
         if sample_weight is not None:
@@ -406,7 +420,8 @@ def fit(self, X, y, sample_weight=None, check_input=True):
         self.n_features_in_ = X.shape[1]
 
         # Build and fit the _FittingTree
-        fitting_tree = _FittingTree(self.max_depth, self.n_features_in_, self.classes_, self.epsilon, self.bounds)
+        fitting_tree = _FittingTree(self.max_depth, self.n_features_in_, self.classes_, self.epsilon, self.bounds,
+                                    random_state)
         fitting_tree.build()
         fitting_tree.fit(X, y)
 
@@ -449,16 +464,20 @@ class _FittingTree(DiffprivlibMixin):
     epsilon : float
         Privacy parameter :math:`\epsilon`.
 
-    bounds : tuple, optional
+    bounds : tuple
         Bounds of the data, provided as a tuple of the form (min, max).  `min` and `max` can either be scalars, covering
         the min/max of the entire data.
 
+    random_state : RandomState
+        Controls the randomness of the building and training process: the feature to split at each node, the threshold
+        to split at and the randomisation of the label at each leaf.
+
     """
     _TREE_LEAF = -1
     _TREE_UNDEFINED = -2
     StackNode = namedtuple("StackNode", ["parent", "is_left", "depth", "bounds"])
 
-    def __init__(self, max_depth, n_features, classes, epsilon, bounds):
+    def __init__(self, max_depth, n_features, classes, epsilon, bounds, random_state):
         self.node_count = 0
         self.nodes = []
         self.max_depth = max_depth
@@ -466,6 +485,7 @@ def __init__(self, max_depth, n_features, classes, epsilon, bounds):
         self.classes = classes
         self.epsilon = epsilon
         self.bounds = bounds
+        self.random_state = random_state
 
     def __getstate__(self):
         """Get state of _FittingTree to feed into __setstate__ of sklearn.Tree"""
@@ -502,8 +522,8 @@ def build(self):
                 continue
 
             # We have a decision node, so pick feature and threshold
-            feature = np.random.randint(self.n_features)
-            threshold = np.random.uniform(bounds_lower[feature], bounds_upper[feature])
+            feature = self.random_state.randint(self.n_features)
+            threshold = self.random_state.uniform(bounds_lower[feature], bounds_upper[feature])
 
             left_bounds_upper = bounds_upper.copy()
             left_bounds_upper[feature] = threshold
@@ -532,6 +552,9 @@ def fit(self, X, y):
             Target vector relative to X.
 
         """
+        if not self.nodes:
+            raise ValueError("Fitting Tree must be built before calling fit().")
+
         leaves = self.apply(X)
         unique_leaves = np.unique(leaves)
         values = np.zeros(shape=(self.node_count, 1, len(self.classes)))
@@ -542,7 +565,8 @@ def fit(self, X, y):
             leaf_y = y[idxs]
 
             counts = [np.sum(leaf_y == cls) for cls in self.classes]
-            mech = PermuteAndFlip(epsilon=self.epsilon, sensitivity=1, monotonic=True, utility=counts)
+            mech = PermuteAndFlip(epsilon=self.epsilon, sensitivity=1, monotonic=True, utility=counts,
+                                  random_state=self.random_state)
             values[leaf, 0, mech.randomise()] = 1
 
         # Populate value of empty leaves
@@ -550,7 +574,7 @@ def fit(self, X, y):
             if values[node.node_id].sum() or node.left_child != self._TREE_LEAF:
                 continue
 
-            values[node.node_id, 0, np.random.randint(len(self.classes))] = 1
+            values[node.node_id, 0, self.random_state.randint(len(self.classes))] = 1
 
         self.values_ = values
 
diff --git a/diffprivlib/models/k_means.py b/diffprivlib/models/k_means.py
index 9a5eefd8..20f4283a 100644
--- a/diffprivlib/models/k_means.py
+++ b/diffprivlib/models/k_means.py
@@ -25,7 +25,7 @@
 
 from diffprivlib.accountant import BudgetAccountant
 from diffprivlib.mechanisms import LaplaceBoundedDomain, GeometricFolded
-from diffprivlib.utils import PrivacyLeakWarning
+from diffprivlib.utils import PrivacyLeakWarning, check_random_state
 from diffprivlib.validation import DiffprivlibMixin
 
 
@@ -48,6 +48,10 @@ class KMeans(sk_cluster.KMeans, DiffprivlibMixin):
         the min/max of the entire data, or vectors with one entry per feature.  If not provided, the bounds are computed
         on the data when ``.fit()`` is first called, resulting in a :class:`.PrivacyLeakWarning`.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the model.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     accountant : BudgetAccountant, optional
         Accountant to keep track of privacy budget.
 
@@ -74,8 +78,8 @@ class KMeans(sk_cluster.KMeans, DiffprivlibMixin):
 
     """
 
-    def __init__(self, n_clusters=8, *, epsilon=1.0, bounds=None, accountant=None, **unused_args):
-        super().__init__(n_clusters=n_clusters)
+    def __init__(self, n_clusters=8, *, epsilon=1.0, bounds=None, random_state=None, accountant=None, **unused_args):
+        super().__init__(n_clusters=n_clusters, random_state=random_state)
 
         self.epsilon = epsilon
         self.bounds = bounds
@@ -116,6 +120,8 @@ def fit(self, X, y=None, sample_weight=None):
 
         del y
 
+        random_state = check_random_state(self.random_state)
+
         X = self._validate_data(X, accept_sparse=False, dtype=[np.float64, np.float32])
         n_samples, n_dims = X.shape
 
@@ -133,14 +139,15 @@ def fit(self, X, y=None, sample_weight=None):
         self.bounds = self._check_bounds(self.bounds, n_dims, min_separation=1e-5)
         X = self._clip_to_bounds(X, self.bounds)
 
-        centers = self._init_centers(n_dims)
+        centers = self._init_centers(n_dims, random_state=random_state)
         labels = None
         distances = None
 
         # Run _update_centers first to ensure consistency of `labels` and `centers`, since convergence unlikely
         for _ in range(-1, iters):
             if labels is not None:
-                centers = self._update_centers(X, centers=centers, labels=labels, dims=n_dims, total_iters=iters)
+                centers = self._update_centers(X, centers=centers, labels=labels, dims=n_dims, total_iters=iters,
+                                               random_state=random_state)
 
             distances, labels = self._distances_labels(X, centers)
 
@@ -153,7 +160,7 @@ def fit(self, X, y=None, sample_weight=None):
 
         return self
 
-    def _init_centers(self, dims):
+    def _init_centers(self, dims, random_state):
         if self.bounds_processed is None:
             bounds_processed = np.zeros(shape=(dims, 2))
 
@@ -175,7 +182,7 @@ def _init_centers(self, dims):
                 if cluster >= self.n_clusters:
                     break
 
-                temp_center = np.random.random(dims) * (self.bounds_processed[:, 0] - 2 * cluster_proximity) + \
+                temp_center = random_state.random(dims) * (self.bounds_processed[:, 0] - 2 * cluster_proximity) + \
                     self.bounds_processed[:, 1] + cluster_proximity
 
                 if cluster == 0:
@@ -208,7 +215,7 @@ def _distances_labels(self, X, centers):
         labels = np.argmin(distances, axis=1)
         return distances, labels
 
-    def _update_centers(self, X, centers, labels, dims, total_iters):
+    def _update_centers(self, X, centers, labels, dims, total_iters, random_state):
         """Updates the centers of the KMeans algorithm for the current iteration, while satisfying differential
         privacy.
 
@@ -218,7 +225,8 @@ def _update_centers(self, X, centers, labels, dims, total_iters):
 
         """
         epsilon_0, epsilon_i = self._split_epsilon(dims, total_iters)
-        geometric_mech = GeometricFolded(epsilon=epsilon_0, sensitivity=1, lower=0.5, upper=float("inf"))
+        geometric_mech = GeometricFolded(epsilon=epsilon_0, sensitivity=1, lower=0.5, upper=float("inf"),
+                                         random_state=random_state)
 
         for cluster in range(self.n_clusters):
             if cluster not in labels:
@@ -234,7 +242,7 @@ def _update_centers(self, X, centers, labels, dims, total_iters):
                 laplace_mech = LaplaceBoundedDomain(epsilon=epsilon_i,
                                                     sensitivity=self.bounds[1][i] - self.bounds[0][i],
                                                     lower=noisy_count * self.bounds[0][i],
-                                                    upper=noisy_count * self.bounds[1][i])
+                                                    upper=noisy_count * self.bounds[1][i], random_state=random_state)
                 noisy_sum[i] = laplace_mech.randomise(cluster_sum[i])
 
             centers[cluster, :] = noisy_sum / noisy_count
diff --git a/diffprivlib/models/linear_regression.py b/diffprivlib/models/linear_regression.py
index bc54c3f7..00594feb 100644
--- a/diffprivlib/models/linear_regression.py
+++ b/diffprivlib/models/linear_regression.py
@@ -54,15 +54,17 @@
 from diffprivlib.accountant import BudgetAccountant
 from diffprivlib.mechanisms import Laplace, LaplaceFolded
 from diffprivlib.tools import mean
-from diffprivlib.utils import warn_unused_args, PrivacyLeakWarning
+from diffprivlib.utils import warn_unused_args, PrivacyLeakWarning, check_random_state
 from diffprivlib.validation import check_bounds, clip_to_bounds, DiffprivlibMixin
 
 
 # noinspection PyPep8Naming
 def _preprocess_data(X, y, fit_intercept, epsilon=1.0, bounds_X=None, bounds_y=None, copy=True, check_input=True,
-                     **unused_args):
+                     random_state=None, **unused_args):
     warn_unused_args(unused_args)
 
+    random_state = check_random_state(random_state)
+
     if check_input:
         X = check_array(X, copy=copy, accept_sparse=False, dtype=FLOAT_DTYPES)
     elif copy:
@@ -78,9 +80,11 @@ def _preprocess_data(X, y, fit_intercept, epsilon=1.0, bounds_X=None, bounds_y=N
         X = clip_to_bounds(X, bounds_X)
         y = clip_to_bounds(y, bounds_y)
 
-        X_offset = mean(X, axis=0, bounds=bounds_X, epsilon=epsilon, accountant=BudgetAccountant())
+        X_offset = mean(X, axis=0, bounds=bounds_X, epsilon=epsilon, random_state=random_state,
+                        accountant=BudgetAccountant())
         X -= X_offset
-        y_offset = mean(y, axis=0, bounds=bounds_y, epsilon=epsilon, accountant=BudgetAccountant())
+        y_offset = mean(y, axis=0, bounds=bounds_y, epsilon=epsilon, random_state=random_state,
+                        accountant=BudgetAccountant())
         y = y - y_offset
     else:
         X_offset = np.zeros(X.shape[1], dtype=X.dtype)
@@ -92,7 +96,7 @@ def _preprocess_data(X, y, fit_intercept, epsilon=1.0, bounds_X=None, bounds_y=N
     return X, y, X_offset, y_offset, X_scale
 
 
-def _construct_regression_obj(X, y, bounds_X, bounds_y, epsilon, alpha):
+def _construct_regression_obj(X, y, bounds_X, bounds_y, epsilon, alpha, random_state):
     if y.ndim == 1:
         y = y.reshape(-1, 1)
 
@@ -113,7 +117,8 @@ def get_max_sensitivity(y_lower, y_upper, x_lower, x_upper):
 
     for i in range(n_targets):
         sensitivity = np.abs([bounds_y[0][i], bounds_y[1][i]]).max() ** 2
-        mech = LaplaceFolded(epsilon=local_epsilon, sensitivity=sensitivity, lower=0, upper=float("inf"))
+        mech = LaplaceFolded(epsilon=local_epsilon, sensitivity=sensitivity, lower=0, upper=float("inf"),
+                             random_state=random_state)
         mono_coef_0[i] = mech.randomise(coefs[0][i])
 
     # Randomise 1st-degree monomial coefficients
@@ -122,7 +127,7 @@ def get_max_sensitivity(y_lower, y_upper, x_lower, x_upper):
     for i in range(n_targets):
         for j in range(n_features):
             sensitivity = get_max_sensitivity(bounds_y[0][i], bounds_y[1][i], bounds_X[0][j], bounds_X[1][j])
-            mech = Laplace(epsilon=local_epsilon, sensitivity=sensitivity)
+            mech = Laplace(epsilon=local_epsilon, sensitivity=sensitivity, random_state=random_state)
             mono_coef_1[j, i] = mech.randomise(coefs[1][j, i])
 
     # Randomise 2nd-degree monomial coefficients
@@ -130,12 +135,13 @@ def get_max_sensitivity(y_lower, y_upper, x_lower, x_upper):
 
     for i in range(n_features):
         sensitivity = np.max(np.abs([bounds_X[0][i], bounds_X[0][i]])) ** 2
-        mech = LaplaceFolded(epsilon=local_epsilon, sensitivity=sensitivity, lower=0, upper=float("inf"))
+        mech = LaplaceFolded(epsilon=local_epsilon, sensitivity=sensitivity, lower=0, upper=float("inf"),
+                             random_state=random_state)
         mono_coef_2[i, i] = mech.randomise(coefs[2][i, i])
 
         for j in range(i + 1, n_features):
             sensitivity = get_max_sensitivity(bounds_X[0][i], bounds_X[1][i], bounds_X[0][j], bounds_X[1][j])
-            mech = Laplace(epsilon=local_epsilon, sensitivity=sensitivity)
+            mech = Laplace(epsilon=local_epsilon, sensitivity=sensitivity, random_state=random_state)
             mono_coef_2[i, j] = mech.randomise(coefs[2][i, j])
             mono_coef_2[j, i] = mono_coef_2[i, j]  # Enforce symmetry
 
@@ -192,6 +198,10 @@ class LinearRegression(sk_lr.LinearRegression, DiffprivlibMixin):
     copy_X : bool, default: True
         If True, X will be copied; else, it may be overwritten.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the model.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     accountant : BudgetAccountant, optional
         Accountant to keep track of privacy budget.
 
@@ -211,13 +221,14 @@ class LinearRegression(sk_lr.LinearRegression, DiffprivlibMixin):
         regression analysis under differential privacy." arXiv preprint arXiv:1208.0219 (2012).
 
     """
-    def __init__(self, *, epsilon=1.0, bounds_X=None, bounds_y=None, fit_intercept=True, copy_X=True, accountant=None,
-                 **unused_args):
+    def __init__(self, *, epsilon=1.0, bounds_X=None, bounds_y=None, fit_intercept=True, copy_X=True, random_state=None,
+                 accountant=None, **unused_args):
         super().__init__(fit_intercept=fit_intercept, copy_X=copy_X, n_jobs=None)
 
         self.epsilon = epsilon
         self.bounds_X = bounds_X
         self.bounds_y = bounds_y
+        self.random_state = random_state
         self.accountant = BudgetAccountant.load_default(accountant)
 
         self._warn_unused_args(unused_args)
@@ -247,6 +258,8 @@ def fit(self, X, y, sample_weight=None):
         if sample_weight is not None:
             self._warn_unused_args("sample_weight")
 
+        random_state = check_random_state(self.random_state)
+
         X, y = self._validate_data(X, y, accept_sparse=False, y_numeric=True, multi_output=True)
 
         if self.bounds_X is None or self.bounds_y is None:
@@ -262,6 +275,7 @@ def fit(self, X, y, sample_weight=None):
             if self.bounds_y is None:
                 self.bounds_y = (np.min(y, axis=0), np.max(y, axis=0))
 
+        # pylint: disable=no-member
         self.bounds_X = self._check_bounds(self.bounds_X, X.shape[1])
         self.bounds_y = self._check_bounds(self.bounds_y, y.shape[1] if y.ndim > 1 else 1)
 
@@ -271,13 +285,14 @@ def fit(self, X, y, sample_weight=None):
 
         X, y, X_offset, y_offset, X_scale = self._preprocess_data(
             X, y, fit_intercept=self.fit_intercept, bounds_X=self.bounds_X, bounds_y=self.bounds_y,
-            epsilon=self.epsilon * epsilon_intercept_scale, copy=self.copy_X)
+            epsilon=self.epsilon * epsilon_intercept_scale, copy=self.copy_X, random_state=random_state)
 
         bounds_X = (self.bounds_X[0] - X_offset, self.bounds_X[1] - X_offset)
         bounds_y = (self.bounds_y[0] - y_offset, self.bounds_y[1] - y_offset)
 
         objs, obj_coefs = _construct_regression_obj(
-            X, y, bounds_X, bounds_y, epsilon=self.epsilon * (1 - epsilon_intercept_scale), alpha=0)
+            X, y, bounds_X, bounds_y, epsilon=self.epsilon * (1 - epsilon_intercept_scale), alpha=0,
+            random_state=random_state)
         coef = np.zeros((n_features, n_targets))
 
         for i, obj in enumerate(objs):
diff --git a/diffprivlib/models/logistic_regression.py b/diffprivlib/models/logistic_regression.py
index 20e7afbf..a775e50d 100644
--- a/diffprivlib/models/logistic_regression.py
+++ b/diffprivlib/models/logistic_regression.py
@@ -49,16 +49,23 @@
 import numpy as np
 from joblib import delayed, Parallel
 from scipy import optimize
-from scipy.special import expit
 from sklearn.exceptions import ConvergenceWarning
 from sklearn import linear_model
 from sklearn.utils import check_array, check_consistent_length
-from sklearn.utils.extmath import safe_sparse_dot, log_logistic
 from sklearn.utils.multiclass import check_classification_targets
 
+# todo: Remove when sklearn v1.1.0 is min requirement
+try:
+    from sklearn.linear_model._linear_loss import LinearModelLoss
+    from sklearn._loss import HalfBinomialLoss
+    SKL_LOSS_MODULE = True
+except (ModuleNotFoundError, ImportError):
+    from sklearn.linear_model._logistic import _logistic_loss_and_grad
+    SKL_LOSS_MODULE = False
+
 from diffprivlib.accountant import BudgetAccountant
 from diffprivlib.mechanisms import Vector
-from diffprivlib.utils import PrivacyLeakWarning, DiffprivlibCompatibilityWarning, warn_unused_args
+from diffprivlib.utils import PrivacyLeakWarning, warn_unused_args, check_random_state
 from diffprivlib.validation import DiffprivlibMixin
 
 
@@ -117,6 +124,10 @@ class LogisticRegression(linear_model.LogisticRegression, DiffprivlibMixin):
         Number of CPU cores used when parallelising over classes.  ``None`` means 1 unless in a context. ``-1`` means
         using all processors.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the model.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     accountant : BudgetAccountant, optional
         Accountant to keep track of privacy budget.
 
@@ -167,10 +178,10 @@ class LogisticRegression(linear_model.LogisticRegression, DiffprivlibMixin):
     """
 
     def __init__(self, *, epsilon=1.0, data_norm=None, tol=1e-4, C=1.0, fit_intercept=True, max_iter=100, verbose=0,
-                 warm_start=False, n_jobs=None, accountant=None, **unused_args):
+                 warm_start=False, n_jobs=None, random_state=None, accountant=None, **unused_args):
         super().__init__(penalty='l2', dual=False, tol=tol, C=C, fit_intercept=fit_intercept, intercept_scaling=1.0,
-                         class_weight=None, random_state=None, solver='lbfgs', max_iter=max_iter, multi_class='ovr',
-                         verbose=verbose, warm_start=warm_start, n_jobs=n_jobs)
+                         class_weight=None, random_state=random_state, solver='lbfgs', max_iter=max_iter,
+                         multi_class='ovr', verbose=verbose, warm_start=warm_start, n_jobs=n_jobs)
         self.epsilon = epsilon
         self.data_norm = data_norm
         self.classes_ = None
@@ -203,6 +214,8 @@ def fit(self, X, y, sample_weight=None):
         if sample_weight is not None:
             self._warn_unused_args("sample_weight")
 
+        random_state = check_random_state(self.random_state)
+
         if not isinstance(self.C, numbers.Real) or self.C < 0:
             raise ValueError(f"Penalty term must be positive; got (C={self.C})")
         if not isinstance(self.max_iter, numbers.Integral) or self.max_iter < 0:
@@ -210,9 +223,8 @@ def fit(self, X, y, sample_weight=None):
         if not isinstance(self.tol, numbers.Real) or self.tol < 0:
             raise ValueError(f"Tolerance for stopping criteria must be positive; got (tol={self.tol})")
 
-        solver = _check_solver(self.solver, self.penalty, self.dual)
         X, y = self._validate_data(X, y, accept_sparse='csr', dtype=float, order="C",
-                                   accept_large_sparse=solver != 'liblinear')
+                                   accept_large_sparse=True)
         check_classification_targets(y)
         self.classes_ = np.unique(y)
         _, n_features = X.shape
@@ -225,8 +237,6 @@ def fit(self, X, y, sample_weight=None):
 
         X = self._clip_to_norm(X, self.data_norm)
 
-        self.multi_class = _check_multi_class(self.multi_class, solver, len(self.classes_))
-
         n_classes = len(self.classes_)
         classes_ = self.classes_
         if n_classes < 2:
@@ -255,7 +265,7 @@ def fit(self, X, y, sample_weight=None):
         fold_coefs_ = Parallel(n_jobs=self.n_jobs, verbose=self.verbose, prefer='processes')(
             path_func(X, y, epsilon=self.epsilon / n_classes, data_norm=self.data_norm, pos_class=class_, Cs=[self.C],
                       fit_intercept=self.fit_intercept, max_iter=self.max_iter, tol=self.tol, verbose=self.verbose,
-                      coef=warm_start_coef_, check_input=False)
+                      coef=warm_start_coef_, random_state=random_state, check_input=False)
             for class_, warm_start_coef_ in zip(classes_, warm_start_coef))
 
         fold_coefs_, _, n_iter_ = zip(*fold_coefs_)
@@ -274,7 +284,7 @@ def fit(self, X, y, sample_weight=None):
 
 
 def _logistic_regression_path(X, y, epsilon, data_norm, pos_class=None, Cs=10, fit_intercept=True, max_iter=100,
-                              tol=1e-4, verbose=0, coef=None, check_input=True, **unused_args):
+                              tol=1e-4, verbose=0, coef=None, random_state=None, check_input=True, **unused_args):
     """Compute a Logistic Regression model with differential privacy for a list of regularization parameters.  Takes
     inspiration from ``_logistic_regression_path`` in scikit-learn, specified to the LBFGS solver and one-vs-rest
     multi class fitting.
@@ -318,6 +328,10 @@ def _logistic_regression_path(X, y, epsilon, data_norm, pos_class=None, Cs=10, f
     coef : array-like, shape (n_features,), optional
         Initialization value for coefficients of logistic regression.  Useless for liblinear solver.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the model.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     check_input : bool, default: True
         If False, the input arrays X and y will not be checked.
 
@@ -337,18 +351,18 @@ def _logistic_regression_path(X, y, epsilon, data_norm, pos_class=None, Cs=10, f
     """
     warn_unused_args(unused_args)
 
+    random_state = check_random_state(random_state)
+
     if isinstance(Cs, numbers.Integral):
         Cs = np.logspace(-4, 4, int(Cs))
 
-    solver = 'lbfgs'
-
     # Data norm increases if intercept is included
     if fit_intercept:
         data_norm = np.sqrt(data_norm ** 2 + 1)
 
     # Pre-processing.
     if check_input:
-        X = check_array(X, accept_sparse='csr', dtype=np.float64, accept_large_sparse=solver != 'liblinear')
+        X = check_array(X, accept_sparse='csr', dtype=np.float64, accept_large_sparse=True)
         y = check_array(y, ensure_2d=False, dtype=None)
         check_consistent_length(X, y)
     _, n_features = X.shape
@@ -367,8 +381,7 @@ def _logistic_regression_path(X, y, epsilon, data_norm, pos_class=None, Cs=10, f
     output_vec = np.zeros(n_features + int(fit_intercept), dtype=X.dtype)
     mask = (y == pos_class)
     y_bin = np.ones(y.shape, dtype=X.dtype)
-    y_bin[~mask] = -1.
-    # for compute_class_weight
+    y_bin[~mask] = 0.0 if SKL_LOSS_MODULE else -1.0
 
     if coef is not None:
         # it must work both giving the bias term and not
@@ -379,17 +392,23 @@ def _logistic_regression_path(X, y, epsilon, data_norm, pos_class=None, Cs=10, f
 
     target = y_bin
 
+    if SKL_LOSS_MODULE:
+        func = LinearModelLoss(base_loss=HalfBinomialLoss(), fit_intercept=fit_intercept).loss_gradient
+    else:
+        func = _logistic_loss_and_grad
+
     coefs = []
     n_iter = np.zeros(len(Cs), dtype=np.int32)
     for i, C in enumerate(Cs):
         vector_mech = Vector(epsilon=epsilon, dimension=n_features + int(fit_intercept), alpha=1. / C,
-                             function_sensitivity=0.25, data_sensitivity=data_norm)
-        noisy_logistic_loss = vector_mech.randomise(_logistic_loss_and_grad)
+                             function_sensitivity=0.25, data_sensitivity=data_norm, random_state=random_state)
+        noisy_logistic_loss = vector_mech.randomise(func)
+
+        args = (X, target, sample_weight, 1. / C) if SKL_LOSS_MODULE else (X, target, 1. / C, sample_weight)
 
         iprint = [-1, 50, 1, 100, 101][np.searchsorted(np.array([0, 1, 2, 3]), verbose)]
         output_vec, _, info = optimize.fmin_l_bfgs_b(noisy_logistic_loss, output_vec, fprime=None,
-                                                     args=(X, target, 1. / C, sample_weight), iprint=iprint, pgtol=tol,
-                                                     maxiter=max_iter)
+                                                     args=args, iprint=iprint, pgtol=tol, maxiter=max_iter)
         if info["warnflag"] == 1:
             warnings.warn("lbfgs failed to converge. Increase the number of iterations.", ConvergenceWarning)
 
@@ -398,115 +417,3 @@ def _logistic_regression_path(X, y, epsilon, data_norm, pos_class=None, Cs=10, f
         n_iter[i] = info['nit']
 
     return np.array(coefs), np.array(Cs), n_iter
-
-
-def _logistic_loss_and_grad(w, X, y, alpha, sample_weight=None):
-    """Computes the logistic loss and gradient.
-
-    Parameters
-    ----------
-    w : ndarray of shape (n_features,) or (n_features + 1,)
-        Coefficient vector.
-
-    X : {array-like, sparse matrix} of shape (n_samples, n_features)
-        Training data.
-
-    y : ndarray of shape (n_samples,)
-        Array of labels.
-
-    alpha : float
-        Regularization parameter. alpha is equal to 1 / C.
-
-    sample_weight : array-like of shape (n_samples,), default=None
-        Array of weights that are assigned to individual samples.  If not provided, then each sample is given unit
-        weight.
-
-    Returns
-    -------
-    out : float
-        Logistic loss.
-
-    grad : ndarray of shape (n_features,) or (n_features + 1,)
-        Logistic gradient.
-
-    """
-    n_samples, n_features = X.shape
-    grad = np.empty_like(w)
-
-    w, _, yz = _intercept_dot(w, X, y)
-
-    if sample_weight is None:
-        sample_weight = np.ones(n_samples)
-
-    # Logistic loss is the negative of the log of the logistic function.
-    out = -np.sum(sample_weight * log_logistic(yz)) + 0.5 * alpha * np.dot(w, w)
-
-    z = expit(yz)
-    z0 = sample_weight * (z - 1) * y
-
-    grad[:n_features] = safe_sparse_dot(X.T, z0) + alpha * w
-
-    # Case where we fit the intercept.
-    if grad.shape[0] > n_features:
-        grad[-1] = z0.sum()
-    return out, grad
-
-
-def _intercept_dot(w, X, y):
-    """Computes y * np.dot(X, w).
-    It takes into consideration if the intercept should be fit or not.
-
-    Parameters
-    ----------
-    w : ndarray of shape (n_features,) or (n_features + 1,)
-        Coefficient vector.
-
-    X : {array-like, sparse matrix} of shape (n_samples, n_features)
-        Training data.
-
-    y : ndarray of shape (n_samples,)
-        Array of labels.
-
-    Returns
-    -------
-    w : ndarray of shape (n_features,)
-        Coefficient vector without the intercept weight (w[-1]) if the intercept should be fit. Unchanged otherwise.
-
-    c : float
-        The intercept.
-
-    yz : float
-        y * np.dot(X, w).
-
-    """
-    c = 0.0
-    if w.size == X.shape[1] + 1:
-        c = w[-1]
-        w = w[:-1]
-
-    z = safe_sparse_dot(X, w) + c
-    yz = y * z
-    return w, c, yz
-
-
-def _check_solver(solver, penalty, dual):
-    if solver != 'lbfgs':
-        warnings.warn("For diffprivlib, solver must be 'lbfgs'.", DiffprivlibCompatibilityWarning)
-        solver = 'lbfgs'
-
-    if penalty != 'l2':
-        raise ValueError(f"Solver {solver} supports only l2 penalties, got {penalty} penalty.")
-    if dual:
-        raise ValueError(f"Solver {solver} supports only dual=False, got dual={dual}")
-    return solver
-
-
-def _check_multi_class(multi_class, solver, n_classes):
-    del solver, n_classes
-
-    if multi_class != 'ovr':
-        warnings.warn(f"For diffprivlib, multi_class must be 'ovr', got {multi_class}.",
-                      DiffprivlibCompatibilityWarning)
-        multi_class = 'ovr'
-
-    return multi_class
diff --git a/diffprivlib/models/naive_bayes.py b/diffprivlib/models/naive_bayes.py
index b6f0468e..290d4ace 100644
--- a/diffprivlib/models/naive_bayes.py
+++ b/diffprivlib/models/naive_bayes.py
@@ -26,7 +26,7 @@
 
 from diffprivlib.accountant import BudgetAccountant
 from diffprivlib.mechanisms import LaplaceBoundedDomain, GeometricTruncated, LaplaceTruncated
-from diffprivlib.utils import PrivacyLeakWarning, warn_unused_args
+from diffprivlib.utils import PrivacyLeakWarning, warn_unused_args, check_random_state
 from diffprivlib.validation import DiffprivlibMixin
 
 
@@ -52,6 +52,10 @@ class GaussianNB(sk_nb.GaussianNB, DiffprivlibMixin):
     var_smoothing : float, default: 1e-9
         Portion of the largest variance of all features that is added to variances for calculation stability.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the model.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     accountant : BudgetAccountant, optional
         Accountant to keep track of privacy budget.
 
@@ -80,11 +84,13 @@ class GaussianNB(sk_nb.GaussianNB, DiffprivlibMixin):
 
     """
 
-    def __init__(self, *, epsilon=1.0, bounds=None, priors=None, var_smoothing=1e-9, accountant=None):
+    def __init__(self, *, epsilon=1.0, bounds=None, priors=None, var_smoothing=1e-9, random_state=None,
+                 accountant=None):
         super().__init__(priors=priors, var_smoothing=var_smoothing)
 
         self.epsilon = epsilon
         self.bounds = bounds
+        self.random_state = random_state
         self.accountant = BudgetAccountant.load_default(accountant)
 
     def _partial_fit(self, X, y, classes=None, _refit=False, sample_weight=None):
@@ -93,6 +99,8 @@ def _partial_fit(self, X, y, classes=None, _refit=False, sample_weight=None):
         if sample_weight is not None:
             warn_unused_args("sample_weight")
 
+        random_state = check_random_state(self.random_state)
+
         X, y = self._validate_data(X, y)
 
         if self.bounds is None:
@@ -146,7 +154,7 @@ def _partial_fit(self, X, y, classes=None, _refit=False, sample_weight=None):
             raise ValueError(f"The target label(s) {unique_y[~unique_y_in_classes]} in y do not exist in the initial "
                              f"classes {classes}")
 
-        noisy_class_counts = self._noisy_class_counts(y)
+        noisy_class_counts = self._noisy_class_counts(y, random_state=random_state)
 
         for _i, y_i in enumerate(unique_y):
             i = classes.searchsorted(y_i)
@@ -154,8 +162,8 @@ def _partial_fit(self, X, y, classes=None, _refit=False, sample_weight=None):
 
             n_i = noisy_class_counts[_i]
 
-            new_theta, new_var = self._update_mean_variance(self.class_count_[i], self.theta_[i, :],
-                                                            self.var_[i, :], X_i, n_noisy=n_i)
+            new_theta, new_var = self._update_mean_variance(self.class_count_[i], self.theta_[i, :], self.var_[i, :],
+                                                            X_i, random_state=random_state, n_noisy=n_i)
 
             self.theta_[i, :] = new_theta
             self.var_[i, :] = new_var
@@ -172,7 +180,7 @@ def _partial_fit(self, X, y, classes=None, _refit=False, sample_weight=None):
 
         return self
 
-    def _update_mean_variance(self, n_past, mu, var, X, sample_weight=None, n_noisy=None):
+    def _update_mean_variance(self, n_past, mu, var, X, random_state, sample_weight=None, n_noisy=None):
         """Compute online update of Gaussian mean and variance.
 
         Given starting sample count, mean, and variance, a new set of points X return the updated mean and variance.
@@ -197,6 +205,10 @@ def _update_mean_variance(self, n_past, mu, var, X, sample_weight=None, n_noisy=
         var : array-like, shape (number of Gaussians,)
             Variances for Gaussians in original set.
 
+        random_state : RandomState
+            Controls the randomness of the model.  To obtain a deterministic behaviour during randomisation,
+            ``random_state`` has to be fixed to an integer.
+
         sample_weight : ignored
             Ignored in diffprivlib.
 
@@ -235,12 +247,12 @@ def _update_mean_variance(self, n_past, mu, var, X, sample_weight=None, n_noisy=
             local_diameter = upper - lower
 
             mech_mu = LaplaceTruncated(epsilon=local_epsilon, delta=0, sensitivity=local_diameter,
-                                       lower=lower * n_noisy, upper=upper * n_noisy)
+                                       lower=lower * n_noisy, upper=upper * n_noisy, random_state=random_state)
             _mu = mech_mu.randomise(temp_x.sum()) / n_noisy
 
             local_sq_sens = max(_mu - lower, upper - _mu) ** 2
             mech_var = LaplaceBoundedDomain(epsilon=local_epsilon, delta=0, sensitivity=local_sq_sens, lower=0,
-                                            upper=local_sq_sens * n_noisy)
+                                            upper=local_sq_sens * n_noisy, random_state=random_state)
             _var = mech_var.randomise(((temp_x - _mu) ** 2).sum()) / n_noisy
 
             new_mu[feature] = _mu
@@ -265,12 +277,13 @@ def _update_mean_variance(self, n_past, mu, var, X, sample_weight=None, n_noisy=
 
         return total_mu, total_var
 
-    def _noisy_class_counts(self, y):
+    def _noisy_class_counts(self, y, random_state):
         unique_y = np.unique(y)
         n_total = y.shape[0]
 
         # Use 1/3 of total epsilon budget for getting noisy class counts
-        mech = GeometricTruncated(epsilon=self.epsilon / 3, sensitivity=1, lower=1, upper=n_total)
+        mech = GeometricTruncated(epsilon=self.epsilon / 3, sensitivity=1, lower=1, upper=n_total,
+                                  random_state=random_state)
         noisy_counts = np.array([mech.randomise((y == y_i).sum()) for y_i in unique_y])
 
         argsort = np.argsort(noisy_counts)
diff --git a/diffprivlib/models/pca.py b/diffprivlib/models/pca.py
index 34a5acba..4b9f0a21 100644
--- a/diffprivlib/models/pca.py
+++ b/diffprivlib/models/pca.py
@@ -53,7 +53,7 @@
 from diffprivlib.accountant import BudgetAccountant
 from diffprivlib.models.utils import covariance_eig
 from diffprivlib.tools import mean
-from diffprivlib.utils import copy_docstring, PrivacyLeakWarning
+from diffprivlib.utils import copy_docstring, PrivacyLeakWarning, check_random_state
 from diffprivlib.validation import DiffprivlibMixin
 
 
@@ -122,9 +122,9 @@ class PCA(sk_pca.PCA, DiffprivlibMixin):
         components) but can sometime improve the predictive accuracy of the downstream estimators by making their
         data respect some hard-wired assumptions.
 
-    random_state : int or RandomState instance, optional
-        If int, random_state is the seed used by the random number generator; If RandomState instance, random_state
-        is the random number generator.
+    random_state : int or RandomState, optional
+        Controls the randomness of the model.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
 
     accountant : BudgetAccountant, optional
         Accountant to keep track of privacy budget.
@@ -199,6 +199,8 @@ def __init__(self, n_components=None, *, epsilon=1.0, data_norm=None, centered=F
     def _fit_full(self, X, n_components):
         self.accountant.check(self.epsilon, 0)
 
+        random_state = check_random_state(self.random_state)
+
         n_samples, n_features = X.shape
 
         if self.centered:
@@ -214,7 +216,8 @@ def _fit_full(self, X, n_components):
                 self.bounds = (np.min(X, axis=0), np.max(X, axis=0))
 
             self.bounds = self._check_bounds(self.bounds, n_features)
-            self.mean_ = mean(X, epsilon=self.epsilon / 2, bounds=self.bounds, axis=0, accountant=BudgetAccountant())
+            self.mean_ = mean(X, epsilon=self.epsilon / 2, bounds=self.bounds, axis=0, random_state=random_state,
+                              accountant=BudgetAccountant())
 
         X -= self.mean_
 
@@ -227,7 +230,7 @@ def _fit_full(self, X, n_components):
         X = self._clip_to_norm(X, self.data_norm)
 
         sigma_vec, u_mtx = covariance_eig(X, epsilon=self.epsilon if self.centered else self.epsilon / 2,
-                                          norm=self.data_norm,
+                                          norm=self.data_norm, random_state=random_state,
                                           dims=n_components if isinstance(n_components, Integral) else None)
         u_mtx, _ = svd_flip(u_mtx, np.zeros_like(u_mtx).T)
         sigma_vec = np.sqrt(sigma_vec)
@@ -242,7 +245,7 @@ def _fit_full(self, X, n_components):
 
         # Post-process the number of components required
         if n_components == 'mle':
-            n_components = sk_pca._infer_dimension(explained_variance_, n_samples)
+            n_components = sk_pca._infer_dimension(explained_variance_, n_samples)  # pylint: disable=protected-access
         elif 0 < n_components < 1.0:
             # number of components for which the cumulated explained
             # variance percentage is superior to the desired threshold
diff --git a/diffprivlib/models/standard_scaler.py b/diffprivlib/models/standard_scaler.py
index 2943ad10..3990323d 100644
--- a/diffprivlib/models/standard_scaler.py
+++ b/diffprivlib/models/standard_scaler.py
@@ -50,12 +50,12 @@
 from sklearn.preprocessing._data import _handle_zeros_in_scale
 
 from diffprivlib.accountant import BudgetAccountant
-from diffprivlib.utils import PrivacyLeakWarning
+from diffprivlib.utils import PrivacyLeakWarning, check_random_state
 from diffprivlib.tools import nanvar, nanmean
 from diffprivlib.validation import DiffprivlibMixin
 
 
-def _incremental_mean_and_var(X, epsilon, bounds, last_mean, last_variance, last_sample_count):
+def _incremental_mean_and_var(X, epsilon, bounds, last_mean, last_variance, last_sample_count, random_state=None):
     # Initialising new accountant, as budget is tracked in main class. Subject to review in line with GH issue #21
     temp_acc = BudgetAccountant()
 
@@ -64,7 +64,7 @@ def _incremental_mean_and_var(X, epsilon, bounds, last_mean, last_variance, last
     # updated = the aggregated stats
     last_sum = last_mean * last_sample_count
 
-    new_mean = nanmean(X, epsilon=epsilon, axis=0, bounds=bounds, accountant=temp_acc)
+    new_mean = nanmean(X, epsilon=epsilon, axis=0, bounds=bounds, random_state=random_state, accountant=temp_acc)
     new_sample_count = np.sum(~np.isnan(X), axis=0)
     new_sum = new_mean * new_sample_count
     updated_sample_count = last_sample_count + new_sample_count
@@ -74,7 +74,7 @@ def _incremental_mean_and_var(X, epsilon, bounds, last_mean, last_variance, last
     if last_variance is None:
         updated_variance = None
     else:
-        new_unnormalized_variance = nanvar(X, epsilon=epsilon, axis=0, bounds=bounds,
+        new_unnormalized_variance = nanvar(X, epsilon=epsilon, axis=0, bounds=bounds, random_state=random_state,
                                            accountant=temp_acc) * new_sample_count
         last_unnormalized_variance = last_variance * last_sample_count
 
@@ -133,6 +133,10 @@ class StandardScaler(sk_pp.StandardScaler, DiffprivlibMixin):
     with_std : boolean, True by default
         If True, scale the data to unit variance (or equivalently, unit standard deviation).
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the model.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     accountant : BudgetAccountant, optional
         Accountant to keep track of privacy budget.
 
@@ -167,10 +171,12 @@ class StandardScaler(sk_pp.StandardScaler, DiffprivlibMixin):
     NaNs are treated as missing values: disregarded in fit, and maintained in transform.
 
     """  # noqa
-    def __init__(self, *, epsilon=1.0, bounds=None, copy=True, with_mean=True, with_std=True, accountant=None):
+    def __init__(self, *, epsilon=1.0, bounds=None, copy=True, with_mean=True, with_std=True, random_state=None,
+                 accountant=None):
         super().__init__(copy=copy, with_mean=with_mean, with_std=with_std)
         self.epsilon = epsilon
         self.bounds = bounds
+        self.random_state = random_state
         self.accountant = BudgetAccountant.load_default(accountant)
 
     def partial_fit(self, X, y=None, sample_weight=None):
@@ -199,6 +205,8 @@ def partial_fit(self, X, y=None, sample_weight=None):
         if sample_weight is not None:
             self._warn_unused_args("sample_weight")
 
+        random_state = check_random_state(self.random_state)
+
         epsilon_0 = self.epsilon / 2 if self.with_std else self.epsilon
 
         X = self._validate_data(X, accept_sparse=False, copy=self.copy, estimator=self, dtype=float,
@@ -239,7 +247,7 @@ def partial_fit(self, X, y=None, sample_weight=None):
             self.n_samples_seen_ += X.shape[0] - np.isnan(X).sum(axis=0)
         else:
             self.mean_, self.var_, self.n_samples_seen_ = _incremental_mean_and_var(
-                X, epsilon_0, self.bounds, self.mean_, self.var_, self.n_samples_seen_
+                X, epsilon_0, self.bounds, self.mean_, self.var_, self.n_samples_seen_, random_state
             )
 
         # for backward-compatibility, reduce n_samples_seen_ to an integer
diff --git a/diffprivlib/models/utils.py b/diffprivlib/models/utils.py
index 8ef3a33a..6595a511 100644
--- a/diffprivlib/models/utils.py
+++ b/diffprivlib/models/utils.py
@@ -25,10 +25,10 @@
 from scipy.linalg import null_space
 
 from diffprivlib.mechanisms import LaplaceBoundedDomain, Bingham
-from diffprivlib.utils import PrivacyLeakWarning
+from diffprivlib.utils import PrivacyLeakWarning, check_random_state
 
 
-def covariance_eig(array, epsilon=1.0, norm=None, dims=None, eigvals_only=False):
+def covariance_eig(array, epsilon=1.0, norm=None, dims=None, eigvals_only=False, random_state=None):
     r"""
     Return the eigenvalues and eigenvectors of the covariance matrix of `array`, satisfying differential privacy.
 
@@ -56,6 +56,10 @@ def covariance_eig(array, epsilon=1.0, norm=None, dims=None, eigvals_only=False)
     eigvals_only : bool, default: False
         Only return the eigenvalue estimates.  If True, all the privacy budget is spent on estimating the eigenvalues.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the model.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     Returns
     -------
     w : (n_features) array
@@ -67,6 +71,8 @@ def covariance_eig(array, epsilon=1.0, norm=None, dims=None, eigvals_only=False)
 
     """
 
+    random_state = check_random_state(random_state)
+
     n_features = array.shape[1]
     dims = n_features if dims is None else min(dims, n_features)
     if not isinstance(dims, Integral):
@@ -87,7 +93,8 @@ def covariance_eig(array, epsilon=1.0, norm=None, dims=None, eigvals_only=False)
     eigvals = np.sort(np.linalg.eigvalsh(cov))[::-1]
     epsilon_0 = epsilon if eigvals_only else epsilon / (dims + (dims != n_features))
 
-    mech_eigvals = LaplaceBoundedDomain(epsilon=epsilon_0, lower=0, upper=float("inf"), sensitivity=2)
+    mech_eigvals = LaplaceBoundedDomain(epsilon=epsilon_0, lower=0, upper=float("inf"), sensitivity=2,
+                                        random_state=random_state)
     noisy_eigvals = np.array([mech_eigvals.randomise(eigval) for eigval in eigvals]) * (norm ** 2)
 
     if eigvals_only:
@@ -99,7 +106,7 @@ def covariance_eig(array, epsilon=1.0, norm=None, dims=None, eigvals_only=False)
     proj_i = np.eye(n_features)
 
     theta = np.zeros((0, n_features))
-    mech_cov = Bingham(epsilon=epsilon_i)
+    mech_cov = Bingham(epsilon=epsilon_i, random_state=random_state)
 
     for _ in range(dims):
         if cov_i.size > 1:
diff --git a/diffprivlib/tools/histograms.py b/diffprivlib/tools/histograms.py
index af33297a..36376c91 100644
--- a/diffprivlib/tools/histograms.py
+++ b/diffprivlib/tools/histograms.py
@@ -50,11 +50,12 @@
 
 from diffprivlib.accountant import BudgetAccountant
 from diffprivlib.mechanisms import GeometricTruncated
-from diffprivlib.utils import PrivacyLeakWarning, warn_unused_args
+from diffprivlib.utils import PrivacyLeakWarning, warn_unused_args, check_random_state
 
 
 # noinspection PyShadowingBuiltins
-def histogram(sample, epsilon=1.0, bins=10, range=None, weights=None, density=None, accountant=None, **unused_args):
+def histogram(sample, epsilon=1.0, bins=10, range=None, weights=None, density=None, random_state=None, accountant=None,
+              **unused_args):
     r"""
     Compute the differentially private histogram of a set of data.
 
@@ -95,6 +96,10 @@ def histogram(sample, epsilon=1.0, bins=10, range=None, weights=None, density=No
         Note that the sum of the histogram values will not be equal to 1 unless bins of unity width are chosen; it is
         not a probability *mass* function.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the algorithm.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     accountant : BudgetAccountant, optional
         Accountant to keep track of privacy budget.
 
@@ -123,6 +128,8 @@ def histogram(sample, epsilon=1.0, bins=10, range=None, weights=None, density=No
     """
     warn_unused_args(unused_args)
 
+    random_state = check_random_state(random_state)
+
     accountant = BudgetAccountant.load_default(accountant)
     accountant.check(epsilon, 0)
 
@@ -133,7 +140,7 @@ def histogram(sample, epsilon=1.0, bins=10, range=None, weights=None, density=No
 
     hist, bin_edges = np.histogram(sample, bins=bins, range=range, weights=weights, density=None)
 
-    dp_mech = GeometricTruncated(epsilon=epsilon, sensitivity=1, lower=0, upper=maxsize)
+    dp_mech = GeometricTruncated(epsilon=epsilon, sensitivity=1, lower=0, upper=maxsize, random_state=random_state)
 
     dp_hist = np.zeros_like(hist)
 
@@ -152,7 +159,8 @@ def histogram(sample, epsilon=1.0, bins=10, range=None, weights=None, density=No
 
 
 # noinspection PyShadowingBuiltins
-def histogramdd(sample, epsilon=1.0, bins=10, range=None, weights=None, density=None, accountant=None, **unused_args):
+def histogramdd(sample, epsilon=1.0, bins=10, range=None, weights=None, density=None, random_state=None,
+                accountant=None, **unused_args):
     r"""
     Compute the differentially private multidimensional histogram of some data.
 
@@ -200,6 +208,10 @@ def histogramdd(sample, epsilon=1.0, bins=10, range=None, weights=None, density=
         True.  If normed is False, the values of the returned histogram are equal to the sum of the weights belonging to
         the samples falling into each bin.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the algorithm.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     accountant : BudgetAccountant, optional
         Accountant to keep track of privacy budget.
 
@@ -218,6 +230,8 @@ def histogramdd(sample, epsilon=1.0, bins=10, range=None, weights=None, density=
     """
     warn_unused_args(unused_args)
 
+    random_state = check_random_state(random_state)
+
     accountant = BudgetAccountant.load_default(accountant)
     accountant.check(epsilon, 0)
 
@@ -232,7 +246,7 @@ def histogramdd(sample, epsilon=1.0, bins=10, range=None, weights=None, density=
 
     hist, bin_edges = np.histogramdd(sample, bins=bins, range=range, normed=None, weights=weights, density=None)
 
-    dp_mech = GeometricTruncated(epsilon=epsilon, sensitivity=1, lower=0, upper=maxsize)
+    dp_mech = GeometricTruncated(epsilon=epsilon, sensitivity=1, lower=0, upper=maxsize, random_state=random_state)
 
     dp_hist = np.zeros_like(hist)
     iterator = np.nditer(hist, flags=['multi_index'])
@@ -262,8 +276,8 @@ def histogramdd(sample, epsilon=1.0, bins=10, range=None, weights=None, density=
 
 
 # noinspection PyShadowingBuiltins
-def histogram2d(array_x, array_y, epsilon=1.0, bins=10, range=None, weights=None, density=None, accountant=None,
-                **unused_args):
+def histogram2d(array_x, array_y, epsilon=1.0, bins=10, range=None, weights=None, density=None, random_state=None,
+                accountant=None, **unused_args):
     r"""
     Compute the differentially private bi-dimensional histogram of two data samples.
 
@@ -301,6 +315,10 @@ def histogram2d(array_x, array_y, epsilon=1.0, bins=10, range=None, weights=None
         True.  If `normed` is False, the values of the returned histogram are equal to the sum of the weights belonging
         to the samples falling into each bin.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the algorithm.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     accountant : BudgetAccountant, optional
         Accountant to keep track of privacy budget.
 
@@ -343,5 +361,5 @@ def histogram2d(array_x, array_y, epsilon=1.0, bins=10, range=None, weights=None
         bins = [xedges, yedges]
 
     hist, edges = histogramdd([array_x, array_y], epsilon=epsilon, bins=bins, range=range, weights=weights,
-                              density=density, accountant=accountant)
+                              density=density, random_state=random_state, accountant=accountant)
     return hist, edges[0], edges[1]
diff --git a/diffprivlib/tools/quantiles.py b/diffprivlib/tools/quantiles.py
index 9a6e635f..0311e5f3 100644
--- a/diffprivlib/tools/quantiles.py
+++ b/diffprivlib/tools/quantiles.py
@@ -24,12 +24,13 @@
 
 from diffprivlib.accountant import BudgetAccountant
 from diffprivlib.mechanisms import Exponential
-from diffprivlib.utils import warn_unused_args, PrivacyLeakWarning
+from diffprivlib.utils import warn_unused_args, PrivacyLeakWarning, check_random_state
 from diffprivlib.validation import clip_to_bounds, check_bounds
 from diffprivlib.tools.utils import _wrap_axis
 
 
-def quantile(array, quant, epsilon=1.0, bounds=None, axis=None, keepdims=False, accountant=None, **unused_args):
+def quantile(array, quant, epsilon=1.0, bounds=None, axis=None, keepdims=False, random_state=None, accountant=None,
+             **unused_args):
     r"""
     Compute the differentially private quantile of the array.
 
@@ -70,6 +71,10 @@ def quantile(array, quant, epsilon=1.0, bounds=None, axis=None, keepdims=False,
         of `ndarray`, however any non-default value will be.  If the sub-class' method does not implement `keepdims` any
         exceptions will be raised.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the algorithm.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     accountant : BudgetAccountant, optional
         Accountant to keep track of privacy budget.
 
@@ -87,6 +92,8 @@ def quantile(array, quant, epsilon=1.0, bounds=None, axis=None, keepdims=False,
     """
     warn_unused_args(unused_args)
 
+    random_state = check_random_state(random_state)
+
     if bounds is None:
         warnings.warn("Bounds have not been specified and will be calculated on the data provided. This will "
                       "result in additional privacy leakage. To ensure differential privacy and no additional "
@@ -107,7 +114,7 @@ def quantile(array, quant, epsilon=1.0, bounds=None, axis=None, keepdims=False,
 
     if axis is not None or keepdims:
         return _wrap_axis(quantile, array, quant=quant, epsilon=epsilon, bounds=bounds, axis=axis, keepdims=keepdims,
-                          accountant=accountant)
+                          random_state=random_state, accountant=accountant)
 
     # Dealing with a scalar output from now on
     bounds = check_bounds(bounds, shape=0, min_separation=1e-5)
@@ -129,16 +136,17 @@ def quantile(array, quant, epsilon=1.0, bounds=None, axis=None, keepdims=False,
         return np.nan
 
     mech = Exponential(epsilon=epsilon, sensitivity=1, utility=list(-np.abs(np.arange(0, k + 1) - quant * k)),
-                       measure=list(interval_sizes))
+                       measure=list(interval_sizes), random_state=random_state)
     idx = mech.randomise()
-    output = mech._rng.random() * (array[idx+1] - array[idx]) + array[idx]
+    output = random_state.random() * (array[idx+1] - array[idx]) + array[idx]
 
     accountant.spend(epsilon, 0)
 
     return output
 
 
-def percentile(array, percent, epsilon=1.0, bounds=None, axis=None, keepdims=False, accountant=None, **unused_args):
+def percentile(array, percent, epsilon=1.0, bounds=None, axis=None, keepdims=False, random_state=None, accountant=None,
+               **unused_args):
     r"""
     Compute the differentially private percentile of the array.
 
@@ -175,6 +183,10 @@ def percentile(array, percent, epsilon=1.0, bounds=None, axis=None, keepdims=Fal
         of `ndarray`, however any non-default value will be.  If the sub-class' method does not implement `keepdims` any
         exceptions will be raised.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the algorithm.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     accountant : BudgetAccountant, optional
         Accountant to keep track of privacy budget.
 
@@ -197,10 +209,12 @@ def percentile(array, percent, epsilon=1.0, bounds=None, axis=None, keepdims=Fal
     if np.any(quant < 0) or np.any(quant > 1):
         raise ValueError("Percentiles must be between 0 and 100 inclusive")
 
-    return quantile(array, quant, epsilon=epsilon, bounds=bounds, axis=axis, keepdims=keepdims, accountant=accountant)
+    return quantile(array, quant, epsilon=epsilon, bounds=bounds, axis=axis, keepdims=keepdims,
+                    random_state=random_state, accountant=accountant)
 
 
-def median(array, epsilon=1.0, bounds=None, axis=None, keepdims=False, accountant=None, **unused_args):
+def median(array, epsilon=1.0, bounds=None, axis=None, keepdims=False, random_state=None, accountant=None,
+           **unused_args):
     r"""
     Compute the differentially private median of the array.
 
@@ -234,6 +248,10 @@ def median(array, epsilon=1.0, bounds=None, axis=None, keepdims=False, accountan
         of `ndarray`, however any non-default value will be.  If the sub-class' method does not implement `keepdims` any
         exceptions will be raised.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the algorithm.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     accountant : BudgetAccountant, optional
         Accountant to keep track of privacy budget.
 
@@ -251,4 +269,5 @@ def median(array, epsilon=1.0, bounds=None, axis=None, keepdims=False, accountan
     """
     warn_unused_args(unused_args)
 
-    return quantile(array, 0.5, epsilon=epsilon, bounds=bounds, axis=axis, keepdims=keepdims, accountant=accountant)
+    return quantile(array, 0.5, epsilon=epsilon, bounds=bounds, axis=axis, keepdims=keepdims, random_state=random_state,
+                    accountant=accountant)
diff --git a/diffprivlib/tools/utils.py b/diffprivlib/tools/utils.py
index 4fe34a1b..99188e0c 100644
--- a/diffprivlib/tools/utils.py
+++ b/diffprivlib/tools/utils.py
@@ -50,7 +50,7 @@
 
 from diffprivlib.accountant import BudgetAccountant
 from diffprivlib.mechanisms import LaplaceBoundedDomain, GeometricTruncated, LaplaceTruncated
-from diffprivlib.utils import PrivacyLeakWarning, warn_unused_args
+from diffprivlib.utils import PrivacyLeakWarning, warn_unused_args, check_random_state
 from diffprivlib.validation import check_bounds, clip_to_bounds
 
 _sum_ = sum
@@ -97,7 +97,7 @@ def _wrap_axis(func, array, *, axis, keepdims, epsilon, bounds, **kwargs):
     return func(array, bounds=bounds, epsilon=epsilon, **kwargs)
 
 
-def count_nonzero(array, epsilon=1.0, accountant=None, axis=None, keepdims=False):
+def count_nonzero(array, epsilon=1.0, axis=None, keepdims=False, random_state=None, accountant=None):
     r"""Counts the number of non-zero values in the array ``array`` with differential privacy.
 
     It is typical to use this function on the result of binary operations, such as ``count_nonzero(array >= 0)``.  If
@@ -117,9 +117,6 @@ def count_nonzero(array, epsilon=1.0, accountant=None, axis=None, keepdims=False
     epsilon : float, default: 1.0
         Privacy parameter :math:`\epsilon`.
 
-    accountant : BudgetAccountant, optional
-        Accountant to keep track of privacy budget.
-
     axis : int or tuple, optional
         Axis or tuple of axes along which to count non-zeros.  Default is None, meaning that non-zeros will be counted
         along a flattened version of ``array``.
@@ -128,6 +125,13 @@ def count_nonzero(array, epsilon=1.0, accountant=None, axis=None, keepdims=False
         If this is set to True, the axes that are counted are left in the result as dimensions with size one. With this
         option, the result will broadcast correctly against the input array.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the algorithm.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
+    accountant : BudgetAccountant, optional
+        Accountant to keep track of privacy budget.
+
     Returns
     -------
     count : int or array of int
@@ -142,11 +146,12 @@ def count_nonzero(array, epsilon=1.0, accountant=None, axis=None, keepdims=False
     else:
         array_bool = array.astype(np.bool_, copy=False)
 
-    return sum(array_bool, axis=axis, dtype=np.intp, bounds=(0, 1), epsilon=epsilon, accountant=accountant,
-               keepdims=keepdims)
+    return sum(array_bool, axis=axis, dtype=np.intp, bounds=(0, 1), epsilon=epsilon, keepdims=keepdims,
+               random_state=random_state, accountant=accountant)
 
 
-def mean(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False, accountant=None, **unused_args):
+def mean(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False, random_state=None, accountant=None,
+         **unused_args):
     r"""
     Compute the differentially private arithmetic mean along the specified axis.
 
@@ -180,6 +185,10 @@ def mean(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False,
         If this is set to True, the axes which are reduced are left in the result as dimensions with size one.  With
         this option, the result will broadcast correctly against the input array.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the algorithm.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     accountant : BudgetAccountant, optional
         Accountant to keep track of privacy budget.
 
@@ -196,10 +205,11 @@ def mean(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False,
     warn_unused_args(unused_args)
 
     return _mean(array, epsilon=epsilon, bounds=bounds, axis=axis, dtype=dtype, keepdims=keepdims,
-                 accountant=accountant, nan=False)
+                 random_state=random_state, accountant=accountant, nan=False)
 
 
-def nanmean(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False, accountant=None, **unused_args):
+def nanmean(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False, random_state=None, accountant=None,
+            **unused_args):
     r"""
     Compute the differentially private arithmetic mean along the specified axis, ignoring NaNs.
 
@@ -235,6 +245,10 @@ def nanmean(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=Fal
         If this is set to True, the axes which are reduced are left in the result as dimensions with size one.  With
         this option, the result will broadcast correctly against the input array.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the algorithm.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     accountant : BudgetAccountant, optional
         Accountant to keep track of privacy budget.
 
@@ -251,10 +265,13 @@ def nanmean(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=Fal
     warn_unused_args(unused_args)
 
     return _mean(array, epsilon=epsilon, bounds=bounds, axis=axis, dtype=dtype, keepdims=keepdims,
-                 accountant=accountant, nan=True)
+                 random_state=random_state, accountant=accountant, nan=True)
+
 
+def _mean(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False, random_state=None,
+          accountant=None, nan=False):
+    random_state = check_random_state(random_state)
 
-def _mean(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False, accountant=None, nan=False):
     if bounds is None:
         warnings.warn("Bounds have not been specified and will be calculated on the data provided. This will "
                       "result in additional privacy leakage. To ensure differential privacy and no additional "
@@ -263,7 +280,7 @@ def _mean(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False
 
     if axis is not None or keepdims:
         return _wrap_axis(_mean, array, epsilon=epsilon, bounds=bounds, axis=axis, dtype=dtype, keepdims=keepdims,
-                          accountant=accountant, nan=nan)
+                          random_state=random_state, accountant=accountant, nan=nan)
 
     lower, upper = check_bounds(bounds, shape=0, dtype=dtype)
 
@@ -276,7 +293,7 @@ def _mean(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False
     actual_mean = _func(array, axis=axis, dtype=dtype, keepdims=keepdims)
 
     mech = LaplaceTruncated(epsilon=epsilon, delta=0, sensitivity=(upper - lower) / array.size, lower=lower,
-                            upper=upper)
+                            upper=upper, random_state=random_state)
     output = mech.randomise(actual_mean)
 
     accountant.spend(epsilon, 0)
@@ -284,7 +301,8 @@ def _mean(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False
     return output
 
 
-def var(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False, accountant=None, **unused_args):
+def var(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False, random_state=None, accountant=None,
+        **unused_args):
     r"""
     Compute the differentially private variance along the specified axis.
 
@@ -320,6 +338,10 @@ def var(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False,
         If this is set to True, the axes which are reduced are left in the result as dimensions with size one.  With
         this option, the result will broadcast correctly against the input array.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the algorithm.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     accountant : BudgetAccountant, optional
         Accountant to keep track of privacy budget.
 
@@ -335,11 +357,12 @@ def var(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False,
     """
     warn_unused_args(unused_args)
 
-    return _var(array, epsilon=epsilon, bounds=bounds, axis=axis, dtype=dtype, keepdims=keepdims, accountant=accountant,
-                nan=False)
+    return _var(array, epsilon=epsilon, bounds=bounds, axis=axis, dtype=dtype, keepdims=keepdims,
+                random_state=random_state, accountant=accountant, nan=False)
 
 
-def nanvar(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False, accountant=None, **unused_args):
+def nanvar(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False, random_state=None, accountant=None,
+           **unused_args):
     r"""
     Compute the differentially private variance along the specified axis, ignoring NaNs.
 
@@ -377,6 +400,10 @@ def nanvar(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=Fals
         If this is set to True, the axes which are reduced are left in the result as dimensions with size one.  With
         this option, the result will broadcast correctly against the input array.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the algorithm.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     accountant : BudgetAccountant, optional
         Accountant to keep track of privacy budget.
 
@@ -393,11 +420,14 @@ def nanvar(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=Fals
     """
     warn_unused_args(unused_args)
 
-    return _var(array, epsilon=epsilon, bounds=bounds, axis=axis, dtype=dtype, keepdims=keepdims, accountant=accountant,
-                nan=True)
+    return _var(array, epsilon=epsilon, bounds=bounds, axis=axis, dtype=dtype, keepdims=keepdims,
+                random_state=random_state, accountant=accountant, nan=True)
 
 
-def _var(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False, accountant=None, nan=False):
+def _var(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False, random_state=None, accountant=None,
+         nan=False):
+    random_state = check_random_state(random_state)
+
     if bounds is None:
         warnings.warn("Bounds have not been specified and will be calculated on the data provided. This will "
                       "result in additional privacy leakage. To ensure differential privacy and no additional "
@@ -406,7 +436,7 @@ def _var(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False,
 
     if axis is not None or keepdims:
         return _wrap_axis(_var, array, epsilon=epsilon, bounds=bounds, axis=axis, dtype=dtype, keepdims=keepdims,
-                          accountant=accountant, nan=nan)
+                          random_state=random_state, accountant=accountant, nan=nan)
 
     lower, upper = check_bounds(bounds, shape=0, dtype=dtype)
 
@@ -421,7 +451,7 @@ def _var(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False,
 
     dp_mech = LaplaceBoundedDomain(epsilon=epsilon, delta=0,
                                    sensitivity=((upper - lower) / array.size) ** 2 * (array.size - 1), lower=0,
-                                   upper=((upper - lower) ** 2) / 4)
+                                   upper=((upper - lower) ** 2) / 4, random_state=random_state)
     output = dp_mech.randomise(actual_var)
 
     accountant.spend(epsilon, 0)
@@ -429,7 +459,8 @@ def _var(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False,
     return output
 
 
-def std(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False, accountant=None, **unused_args):
+def std(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False, random_state=None, accountant=None,
+        **unused_args):
     r"""
     Compute the standard deviation along the specified axis.
 
@@ -465,6 +496,10 @@ def std(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False,
         If this is set to True, the axes which are reduced are left in the result as dimensions with size one.  With
         this option, the result will broadcast correctly against the input array.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the algorithm.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     accountant : BudgetAccountant, optional
         Accountant to keep track of privacy budget.
 
@@ -481,10 +516,11 @@ def std(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False,
     warn_unused_args(unused_args)
 
     return _std(array, epsilon=epsilon, bounds=bounds, axis=axis, dtype=dtype, keepdims=keepdims,
-                accountant=accountant, nan=False)
+                random_state=random_state, accountant=accountant, nan=False)
 
 
-def nanstd(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False, accountant=None, **unused_args):
+def nanstd(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False, random_state=None, accountant=None,
+           **unused_args):
     r"""
     Compute the standard deviation along the specified axis, ignoring NaNs.
 
@@ -522,6 +558,10 @@ def nanstd(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=Fals
         If this is set to True, the axes which are reduced are left in the result as dimensions with size one.  With
         this option, the result will broadcast correctly against the input array.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the algorithm.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
     accountant : BudgetAccountant, optional
         Accountant to keep track of privacy budget.
 
@@ -537,13 +577,14 @@ def nanstd(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=Fals
     """
     warn_unused_args(unused_args)
 
-    return _std(array, epsilon=epsilon, bounds=bounds, axis=axis, dtype=dtype, keepdims=keepdims, accountant=accountant,
-                nan=True)
+    return _std(array, epsilon=epsilon, bounds=bounds, axis=axis, dtype=dtype, keepdims=keepdims,
+                random_state=random_state, accountant=accountant, nan=True)
 
 
-def _std(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False, accountant=None, nan=False):
-    ret = _var(array, epsilon=epsilon, bounds=bounds, axis=axis, dtype=dtype, keepdims=keepdims, accountant=accountant,
-               nan=nan)
+def _std(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False, random_state=None, accountant=None,
+         nan=False):
+    ret = _var(array, epsilon=epsilon, bounds=bounds, axis=axis, dtype=dtype, keepdims=keepdims,
+               random_state=random_state, accountant=accountant, nan=nan)
 
     if isinstance(ret, mu.ndarray):
         ret = um.sqrt(ret)
@@ -555,7 +596,8 @@ def _std(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False,
     return ret
 
 
-def sum(array, epsilon=1.0, bounds=None, accountant=None, axis=None, dtype=None, keepdims=False, **unused_args):
+def sum(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False, random_state=None, accountant=None,
+        **unused_args):
     r"""Sum of array elements over a given axis with differential privacy.
 
     Parameters
@@ -569,9 +611,6 @@ def sum(array, epsilon=1.0, bounds=None, accountant=None, axis=None, dtype=None,
     bounds : tuple, optional
         Bounds of the values of the array, of the form (min, max).
 
-    accountant : BudgetAccountant, optional
-        Accountant to keep track of privacy budget.
-
     axis : None or int or tuple of ints, optional
         Axis or axes along which a sum is performed.  The default, axis=None, will sum all of the elements of the input
         array.  If axis is negative it counts from the last to the first axis.
@@ -589,6 +628,13 @@ def sum(array, epsilon=1.0, bounds=None, accountant=None, axis=None, dtype=None,
         If this is set to True, the axes which are reduced are left in the result as dimensions with size one.  With
         this option, the result will broadcast correctly against the input array.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the algorithm.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
+    accountant : BudgetAccountant, optional
+        Accountant to keep track of privacy budget.
+
     Returns
     -------
     sum_along_axis : ndarray
@@ -604,11 +650,12 @@ def sum(array, epsilon=1.0, bounds=None, accountant=None, axis=None, dtype=None,
     """
     warn_unused_args(unused_args)
 
-    return _sum(array, epsilon=epsilon, bounds=bounds, accountant=accountant, axis=axis, dtype=dtype, keepdims=keepdims,
-                nan=False)
+    return _sum(array, epsilon=epsilon, bounds=bounds, axis=axis, dtype=dtype, keepdims=keepdims,
+                random_state=random_state, accountant=accountant, nan=False)
 
 
-def nansum(array, epsilon=1.0, bounds=None, accountant=None, axis=None, dtype=None, keepdims=False, **unused_args):
+def nansum(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False, random_state=None, accountant=None,
+           **unused_args):
     r"""Sum of array elements over a given axis with differential privacy, ignoring NaNs.
 
     Parameters
@@ -622,9 +669,6 @@ def nansum(array, epsilon=1.0, bounds=None, accountant=None, axis=None, dtype=No
     bounds : tuple, optional
         Bounds of the values of the array, of the form (min, max).
 
-    accountant : BudgetAccountant, optional
-        Accountant to keep track of privacy budget.
-
     axis : None or int or tuple of ints, optional
         Axis or axes along which a sum is performed.  The default, axis=None, will sum all of the elements of the input
         array.  If axis is negative it counts from the last to the first axis.
@@ -642,6 +686,13 @@ def nansum(array, epsilon=1.0, bounds=None, accountant=None, axis=None, dtype=No
         If this is set to True, the axes which are reduced are left in the result as dimensions with size one.  With
         this option, the result will broadcast correctly against the input array.
 
+    random_state : int or RandomState, optional
+        Controls the randomness of the algorithm.  To obtain a deterministic behaviour during randomisation,
+        ``random_state`` has to be fixed to an integer.
+
+    accountant : BudgetAccountant, optional
+        Accountant to keep track of privacy budget.
+
     Returns
     -------
     sum_along_axis : ndarray
@@ -657,11 +708,14 @@ def nansum(array, epsilon=1.0, bounds=None, accountant=None, axis=None, dtype=No
     """
     warn_unused_args(unused_args)
 
-    return _sum(array, epsilon=epsilon, bounds=bounds, accountant=accountant, axis=axis, dtype=dtype,
-                keepdims=keepdims, nan=True)
+    return _sum(array, epsilon=epsilon, bounds=bounds, axis=axis, dtype=dtype, keepdims=keepdims,
+                random_state=random_state, accountant=accountant, nan=True)
+
 
+def _sum(array, epsilon=1.0, bounds=None, axis=None, dtype=None, keepdims=False, random_state=None, accountant=None,
+         nan=False):
+    random_state = check_random_state(random_state)
 
-def _sum(array, epsilon=1.0, bounds=None, accountant=None, axis=None, dtype=None, keepdims=False, nan=False):
     if bounds is None:
         warnings.warn("Bounds have not been specified and will be calculated on the data provided. This will "
                       "result in additional privacy leakage. To ensure differential privacy and no additional "
@@ -669,8 +723,8 @@ def _sum(array, epsilon=1.0, bounds=None, accountant=None, axis=None, dtype=None
         bounds = (np.min(array), np.max(array))
 
     if axis is not None or keepdims:
-        return _wrap_axis(_sum, array, epsilon=epsilon, bounds=bounds, accountant=accountant, axis=axis, dtype=dtype,
-                          keepdims=keepdims, nan=nan)
+        return _wrap_axis(_sum, array, epsilon=epsilon, bounds=bounds, axis=axis, dtype=dtype, keepdims=keepdims,
+                          random_state=random_state, accountant=accountant, nan=nan)
 
     lower, upper = check_bounds(bounds, shape=0, dtype=dtype)
 
@@ -684,7 +738,8 @@ def _sum(array, epsilon=1.0, bounds=None, accountant=None, axis=None, dtype=None
     actual_sum = _func(array, axis=axis, dtype=dtype, keepdims=keepdims)
 
     mech = GeometricTruncated if dtype is not None and issubclass(dtype, Integral) else LaplaceTruncated
-    mech = mech(epsilon=epsilon, sensitivity=upper - lower, lower=lower * array.size, upper=upper * array.size)
+    mech = mech(epsilon=epsilon, sensitivity=upper - lower, lower=lower * array.size, upper=upper * array.size,
+                random_state=random_state)
     output = mech.randomise(actual_sum)
 
     accountant.spend(epsilon, 0)
diff --git a/diffprivlib/utils.py b/diffprivlib/utils.py
index 04bf09ad..467bc899 100644
--- a/diffprivlib/utils.py
+++ b/diffprivlib/utils.py
@@ -18,25 +18,11 @@
 """
 Basic functions and other utilities for the differential privacy library
 """
+import secrets
 import warnings
 
 import numpy as np
-
-
-def global_seed(seed):
-    """Sets the seed for all random number generators, to guarantee reproducibility in experiments.
-
-    Parameters
-    ----------
-    seed : int
-        The seed value for the random number generators.
-
-    Returns
-    -------
-    None
-
-    """
-    np.random.seed(seed)
+from sklearn.utils import check_random_state as skl_check_random_state
 
 
 def copy_docstring(source):
@@ -86,6 +72,36 @@ def warn_unused_args(args):
                       "warning.", DiffprivlibCompatibilityWarning)
 
 
+def check_random_state(seed, secure=False):
+    """Turn seed into a np.random.RandomState or secrets.SystemRandom instance.
+
+    If secure=True, and seed is None (or was generated from a previous None seed), then secrets is used.  Otherwise a
+    np.random.RandomState is used.
+
+    Parameters
+    ----------
+    seed : None, int or instance of RandomState
+        If seed is None and secure is False, return the RandomState singleton used by np.random.
+        If seed is None and secure is True, return a SystemRandom instance from secrets.
+        If seed is an int, return a new RandomState instance seeded with seed.
+        If seed is already a RandomState or SystemRandom instance, return it.
+        Otherwise raise ValueError.
+
+    secure : bool, default: False
+        Specifies if a secure random number generator from secrets can be used.
+    """
+    if secure:
+        if isinstance(seed, secrets.SystemRandom):
+            return seed
+
+        if seed is None or seed is np.random.mtrand._rand:  # pylint: disable=protected-access
+            return secrets.SystemRandom()
+    elif isinstance(seed, secrets.SystemRandom):
+        raise ValueError("secrets.SystemRandom instance cannot be passed when secure is False.")
+
+    return skl_check_random_state(seed)
+
+
 class Budget(tuple):
     """Custom tuple subclass for privacy budgets of the form (epsilon, delta).
 
@@ -106,9 +122,12 @@ class Budget(tuple):
 
     """
     def __new__(cls, epsilon, delta):
-        from diffprivlib.validation import check_epsilon_delta
+        if epsilon < 0:
+            raise ValueError("Epsilon must be non-negative")
+
+        if not 0 <= delta <= 1:
+            raise ValueError("Delta must be in [0, 1]")
 
-        check_epsilon_delta(epsilon, delta, allow_zero=True)
         return tuple.__new__(cls, (epsilon, delta))
 
     def __gt__(self, other):
diff --git a/diffprivlib/validation.py b/diffprivlib/validation.py
index 13532b56..53b52033 100644
--- a/diffprivlib/validation.py
+++ b/diffprivlib/validation.py
@@ -200,7 +200,7 @@ def clip_to_bounds(array, bounds):
     return clipped_array
 
 
-class DiffprivlibMixin:
+class DiffprivlibMixin:  # pylint: disable=too-few-public-methods
     """Mixin for Diffprivlib models."""
     _check_bounds = staticmethod(check_bounds)
     _clip_to_norm = staticmethod(clip_to_norm)
diff --git a/docs/modules/utils.rst b/docs/modules/utils.rst
index fa292c30..0745fef9 100644
--- a/docs/modules/utils.rst
+++ b/docs/modules/utils.rst
@@ -17,6 +17,6 @@ General classes
 General functions
 -----------------
 
+.. autofunction:: check_random_state
 .. autofunction:: copy_docstring
-.. autofunction:: global_seed
 .. autofunction:: warn_unused_args


Test patch:
diff --git a/tests/mechanisms/test_BernoulliNegExp.py b/tests/mechanisms/test_BernoulliNegExp.py
index 3f1ba802..9b409e99 100644
--- a/tests/mechanisms/test_BernoulliNegExp.py
+++ b/tests/mechanisms/test_BernoulliNegExp.py
@@ -30,7 +30,7 @@ def test_no_rng(self):
             self.assertIn(bernoulli_neg_exp(1), [0, 1])
 
     def test_bad_rng(self):
-        self.assertRaises(AttributeError, bernoulli_neg_exp, 1, "rng")
+        self.assertRaises(ValueError, bernoulli_neg_exp, 1, "rng")
 
     def test_output_domain(self):
         for i in range(1000):
@@ -43,3 +43,14 @@ def test_bernoulli_exp_prob(self):
                                0.5, delta=0.025)
         self.assertAlmostEqual(sum(bernoulli_neg_exp(-np.log(0.1), self.rng) for _ in range(runs)) / runs,
                                0.1, delta=0.025)
+
+    def test_random_state(self):
+        rng_42a, rng_42b = np.random.RandomState(42), np.random.RandomState(42)
+        self.assertEqual([bernoulli_neg_exp(1, rng_42a) for _ in range(100)],
+                         [bernoulli_neg_exp(1, rng_42b) for _ in range(100)])
+
+        rng_42a = np.random.RandomState(42)
+        self.assertEqual([bernoulli_neg_exp(1, rng_42a)] * 100, [bernoulli_neg_exp(1, 42) for _ in range(100)])
+
+        self.assertNotEqual([bernoulli_neg_exp(1, rng_42a) for _ in range(100)],
+                            [bernoulli_neg_exp(1, rng_42a) for _ in range(100)])
diff --git a/tests/mechanisms/test_Binary.py b/tests/mechanisms/test_Binary.py
index 852dc4d8..a47cda3d 100644
--- a/tests/mechanisms/test_Binary.py
+++ b/tests/mechanisms/test_Binary.py
@@ -2,14 +2,10 @@
 from unittest import TestCase
 
 from diffprivlib.mechanisms import Binary
-from diffprivlib.utils import global_seed
 
 
 class TestBinary(TestCase):
     def setup_method(self, method):
-        if method.__name__ .endswith("prob"):
-            global_seed(314159)
-
         self.mech = Binary
 
     def teardown_method(self, method):
@@ -59,7 +55,7 @@ def test_randomise_wrong_label(self):
     def test_distrib_prob(self):
         epsilon = np.log(2)
         runs = 20000
-        mech = self.mech(epsilon=epsilon, value0="0", value1="1")
+        mech = self.mech(epsilon=epsilon, value0="0", value1="1", random_state=0)
         count = [0, 0]
 
         for i in range(runs):
@@ -69,6 +65,16 @@ def test_distrib_prob(self):
         # print("%d / %d = %f" % (count[0], count[1], count[0] / count[1]))
         self.assertAlmostEqual(count[0] / count[1], np.exp(epsilon), delta=0.1)
 
+    def test_random_state(self):
+        mech1 = self.mech(epsilon=1, value0="0", value1="1", random_state=42)
+        mech2 = self.mech(epsilon=1, value0="0", value1="1", random_state=42)
+        self.assertEqual([mech1.randomise("1") for _ in range(100)], [mech2.randomise("1") for _ in range(100)])
+
+        self.assertNotEqual([mech1.randomise("1")] * 100, [mech1.randomise("1") for _ in range(100)])
+
+        mech2 = self.mech(epsilon=1, value0="0", value1="1", random_state=np.random.RandomState(42))
+        self.assertNotEqual([mech1.randomise("1") for _ in range(100)], [mech2.randomise("1") for _ in range(100)])
+
     def test_repr(self):
         repr_ = repr(self.mech(epsilon=1, value0="0", value1="1"))
         self.assertIn(".Binary(", repr_)
diff --git a/tests/mechanisms/test_Bingham.py b/tests/mechanisms/test_Bingham.py
index ed8831a4..f96418d3 100644
--- a/tests/mechanisms/test_Bingham.py
+++ b/tests/mechanisms/test_Bingham.py
@@ -2,14 +2,11 @@
 from unittest import TestCase
 
 from diffprivlib.mechanisms import Bingham
-from diffprivlib.utils import global_seed
+from diffprivlib.utils import check_random_state
 
 
 class TestBingham(TestCase):
     def setup_method(self, method):
-        if method.__name__ .endswith("prob"):
-            global_seed(314159)
-
         self.mech = Bingham
         a = np.random.random((5, 3))
         self.random_array = a.T.dot(a)
@@ -132,12 +129,13 @@ def test_3D_input(self):
             mech.randomise(np.ones((3, 3, 3)))
 
     def test_large_input(self):
-        X = np.random.randn(10000, 21)
+        rng = check_random_state(0)
+        X = rng.randn(5000, 21)
         X -= np.mean(X, axis=0)
         X /= np.linalg.norm(X, axis=1).max()
         XtX = X.T.dot(X)
 
-        mech = self.mech(epsilon=1)
+        mech = self.mech(epsilon=1, random_state=rng)
         self.assertIsNotNone(mech.randomise(XtX))
 
     def test_different_result(self):
@@ -152,6 +150,19 @@ def test_different_result(self):
             self.assertTrue(np.isclose(noisy_data.dot(noisy_data), 1.0))
             self.assertFalse(np.allclose(noisy_data, old_noisy_data))
 
+    def test_random_state(self):
+        mech1 = self.mech(epsilon=1, sensitivity=1, random_state=0)
+        mech2 = self.mech(epsilon=1, sensitivity=1, random_state=0)
+        self.assertTrue(np.allclose(mech1.randomise(self.random_array), mech2.randomise(self.random_array)))
+
+        mech1 = self.mech(epsilon=1, sensitivity=1, random_state=0)
+        self.assertFalse(np.allclose(mech1.randomise(self.random_array), mech2.randomise(self.random_array)))
+
+        rng = np.random.RandomState(0)
+        mech1 = self.mech(epsilon=1, sensitivity=1, random_state=rng)
+        mech2 = self.mech(epsilon=1, sensitivity=1, random_state=rng)
+        self.assertFalse(np.allclose(mech1.randomise(self.random_array), mech2.randomise(self.random_array)))
+
     def test_repr(self):
         repr_ = repr(self.mech(epsilon=1, sensitivity=1))
         self.assertIn(".Bingham(", repr_)
diff --git a/tests/mechanisms/test_DPMechanism.py b/tests/mechanisms/test_DPMechanism.py
index 30774082..f6db774f 100644
--- a/tests/mechanisms/test_DPMechanism.py
+++ b/tests/mechanisms/test_DPMechanism.py
@@ -1,5 +1,9 @@
+import secrets
 from unittest import TestCase
 
+import numpy as np
+from sklearn.utils import check_random_state
+
 from diffprivlib.mechanisms import DPMachine, DPMechanism
 
 
@@ -59,3 +63,19 @@ def variance(self, value):
         mse_mech = TestMSE(epsilon=1, delta=0)
 
         self.assertEqual(mse_mech.mse(1), 2)
+
+    def test_random_state(self):
+        rng = check_random_state(None)
+        rng2 = check_random_state(None)
+        self.assertIs(rng, np.random.mtrand._rand)
+        self.assertIs(rng, rng2)
+
+        mech = self.mech(epsilon=1, delta=0, random_state=None)
+        self.assertIsInstance(mech._rng, secrets.SystemRandom)
+
+        mech = self.mech(epsilon=1, delta=0, random_state=rng)
+        self.assertIsInstance(mech._rng, secrets.SystemRandom)
+
+        mech = self.mech(epsilon=1, delta=0, random_state=1)
+        self.assertNotIsInstance(mech._rng, secrets.SystemRandom)
+        self.assertIsInstance(mech._rng, np.random.RandomState)
diff --git a/tests/mechanisms/test_Exponential.py b/tests/mechanisms/test_Exponential.py
index 81320747..eb6128b3 100644
--- a/tests/mechanisms/test_Exponential.py
+++ b/tests/mechanisms/test_Exponential.py
@@ -2,14 +2,10 @@
 from unittest import TestCase
 
 from diffprivlib.mechanisms import Exponential
-from diffprivlib.utils import global_seed
 
 
 class TestExponential(TestCase):
     def setup_method(self, method):
-        if method.__name__ .endswith("prob"):
-            global_seed(314159)
-
         self.mech = Exponential
 
     def teardown_method(self, method):
@@ -125,7 +121,7 @@ def test_non_uniform_measure(self):
         measure = [1, 2, 1, 1]
         utility = [1, 1, 0, 0]
         runs = 20000
-        mech = self.mech(epsilon=2 * np.log(2), utility=utility, measure=measure, sensitivity=1)
+        mech = self.mech(epsilon=2 * np.log(2), utility=utility, measure=measure, sensitivity=1, random_state=0)
         count = [0] * 4
 
         for i in range(runs):
@@ -158,9 +154,10 @@ def test_inf_utility_measure(self):
 
     def test_distrib_prob(self):
         epsilon = np.log(2)
-        runs = 20000
-        mech1 = self.mech(epsilon=epsilon, utility=[2, 1, 0], sensitivity=1, monotonic=False)
-        mech2 = self.mech(epsilon=epsilon, utility=[2, 1, 1], sensitivity=1, monotonic=False)
+        runs = 2000
+        rng = np.random.RandomState(42)
+        mech1 = self.mech(epsilon=epsilon, utility=[2, 1, 0], sensitivity=1, monotonic=False, random_state=rng)
+        mech2 = self.mech(epsilon=epsilon, utility=[2, 1, 1], sensitivity=1, monotonic=False, random_state=rng)
         counts = np.zeros((2, 3))
 
         for i in range(runs):
@@ -173,9 +170,10 @@ def test_distrib_prob(self):
 
     def test_monotonic_distrib(self):
         epsilon = np.log(2)
-        runs = 20000
-        mech1 = self.mech(epsilon=epsilon, utility=[2, 1, 0], sensitivity=1, monotonic=True)
-        mech2 = self.mech(epsilon=epsilon, utility=[2, 1, 1], sensitivity=1, monotonic=True)
+        runs = 2000
+        rng = np.random.RandomState(42)
+        mech1 = self.mech(epsilon=epsilon, utility=[2, 1, 0], sensitivity=1, monotonic=True, random_state=rng)
+        mech2 = self.mech(epsilon=epsilon, utility=[2, 1, 1], sensitivity=1, monotonic=True, random_state=rng)
         counts = np.zeros((2, 3))
 
         for i in range(runs):
@@ -183,9 +181,18 @@ def test_monotonic_distrib(self):
             counts[1, mech2.randomise()] += 1
 
         for vec in counts.T:
-            # print(vec.max() / vec.min())
             self.assertLessEqual(vec.max() / vec.min(), np.exp(epsilon) + 0.1)
 
+    def test_random_state(self):
+        mech1 = self.mech(epsilon=1, utility=[2, 1, 0], sensitivity=1, random_state=42)
+        mech2 = self.mech(epsilon=1, utility=[2, 1, 0], sensitivity=1, random_state=42)
+        self.assertEqual([mech1.randomise() for _ in range(100)], [mech2.randomise() for _ in range(100)])
+
+        self.assertNotEqual([mech1.randomise()] * 100, [mech1.randomise() for _ in range(100)])
+
+        mech2 = self.mech(epsilon=1, utility=[2, 1, 0], sensitivity=1, random_state=np.random.RandomState(0))
+        self.assertNotEqual([mech1.randomise() for _ in range(100)], [mech2.randomise() for _ in range(100)])
+
     def test_repr(self):
         repr_ = repr(self.mech(epsilon=1, utility=[1], sensitivity=1))
         self.assertIn(".Exponential(", repr_)
diff --git a/tests/mechanisms/test_ExponentialCategorical.py b/tests/mechanisms/test_ExponentialCategorical.py
index 6e48db0c..7b5e5c2d 100644
--- a/tests/mechanisms/test_ExponentialCategorical.py
+++ b/tests/mechanisms/test_ExponentialCategorical.py
@@ -2,14 +2,10 @@
 from unittest import TestCase
 
 from diffprivlib.mechanisms import ExponentialCategorical
-from diffprivlib.utils import global_seed
 
 
 class TestExponential(TestCase):
     def setup_method(self, method):
-        if method.__name__ .endswith("prob"):
-            global_seed(314159)
-
         self.mech = ExponentialCategorical
 
     def teardown_method(self, method):
@@ -156,6 +152,22 @@ def test_distrib_prob(self):
         self.assertLessEqual(count[0] / runs, np.exp(epsilon) * count[2] / runs + 0.05)
         self.assertAlmostEqual(count[0] / count[1], count[1] / count[2], delta=0.15)
 
+    def test_random_state(self):
+        utility_list = [
+            ["A", "B", 1],
+            ["A", "C", 2],
+            ["B", "C", 2]
+        ]
+
+        mech1 = self.mech(epsilon=1, utility_list=utility_list, random_state=42)
+        mech2 = self.mech(epsilon=1, utility_list=utility_list, random_state=42)
+        self.assertEqual([mech1.randomise("A") for _ in range(100)], [mech2.randomise("A") for _ in range(100)])
+
+        self.assertNotEqual([mech1.randomise("A")] * 100, [mech1.randomise("A") for _ in range(100)])
+
+        mech2 = self.mech(epsilon=1, utility_list=utility_list, random_state=np.random.RandomState(0))
+        self.assertNotEqual([mech1.randomise("A") for _ in range(100)], [mech2.randomise("A") for _ in range(100)])
+
     def test_repr(self):
         repr_ = repr(self.mech(epsilon=1, utility_list=[]))
         self.assertIn(".ExponentialCategorical(", repr_)
diff --git a/tests/mechanisms/test_ExponentialHierarchical.py b/tests/mechanisms/test_ExponentialHierarchical.py
index 069a7eea..86b060a4 100644
--- a/tests/mechanisms/test_ExponentialHierarchical.py
+++ b/tests/mechanisms/test_ExponentialHierarchical.py
@@ -2,14 +2,10 @@
 from unittest import TestCase
 
 from diffprivlib.mechanisms import ExponentialHierarchical
-from diffprivlib.utils import global_seed
 
 
 class TestExponentialHierarchical(TestCase):
     def setup_method(self, method):
-        if method.__name__ .endswith("prob"):
-            global_seed(314159)
-
         self.mech = ExponentialHierarchical
 
     def teardown_method(self, method):
@@ -67,7 +63,7 @@ def test_distrib_prob(self):
         epsilon = np.log(2)
         runs = 20000
         balanced_tree = False
-        mech = self.mech(epsilon=epsilon, hierarchy=[["A", "B"], ["C"]])
+        mech = self.mech(epsilon=epsilon, hierarchy=[["A", "B"], ["C"]], random_state=0)
         count = [0, 0, 0]
 
         for i in range(runs):
@@ -90,7 +86,7 @@ def test_distrib_prob(self):
     def test_neighbours_prob(self):
         epsilon = np.log(2)
         runs = 20000
-        mech = self.mech(epsilon=epsilon, hierarchy=[["A", "B"], ["C"]])
+        mech = self.mech(epsilon=epsilon, hierarchy=[["A", "B"], ["C"]], random_state=0)
         count = [0, 0, 0]
 
         for i in range(runs):
@@ -114,7 +110,7 @@ def test_neighbours_prob(self):
     def test_neighbours_flat_hierarchy_prob(self):
         epsilon = np.log(2)
         runs = 20000
-        mech = self.mech(epsilon=epsilon, hierarchy=["A", "B", "C"])
+        mech = self.mech(epsilon=epsilon, hierarchy=["A", "B", "C"], random_state=0)
         count = [0, 0, 0]
 
         for i in range(runs):
@@ -135,6 +131,16 @@ def test_neighbours_flat_hierarchy_prob(self):
         self.assertLessEqual(count[0] / runs, np.exp(epsilon) * count[2] / runs + 0.05)
         self.assertLessEqual(count[1] / runs, np.exp(epsilon) * count[2] / runs + 0.05)
 
+    def test_random_state(self):
+        mech1 = self.mech(epsilon=1, hierarchy=["A", "B", "C"], random_state=42)
+        mech2 = self.mech(epsilon=1, hierarchy=["A", "B", "C"], random_state=42)
+        self.assertEqual([mech1.randomise("A") for _ in range(100)], [mech2.randomise("A") for _ in range(100)])
+
+        self.assertNotEqual([mech1.randomise("A")] * 100, [mech1.randomise("A") for _ in range(100)])
+
+        mech2 = self.mech(epsilon=1, hierarchy=["A", "B", "C"], random_state=np.random.RandomState(0))
+        self.assertNotEqual([mech1.randomise("A") for _ in range(100)], [mech2.randomise("A") for _ in range(100)])
+
     def test_repr(self):
         repr_ = repr(self.mech(epsilon=1, hierarchy=[]))
         self.assertIn(".ExponentialHierarchical(", repr_)
diff --git a/tests/mechanisms/test_Gaussian.py b/tests/mechanisms/test_Gaussian.py
index a7dc1ab0..0e2a4700 100644
--- a/tests/mechanisms/test_Gaussian.py
+++ b/tests/mechanisms/test_Gaussian.py
@@ -2,14 +2,10 @@
 from unittest import TestCase
 
 from diffprivlib.mechanisms import Gaussian
-from diffprivlib.utils import global_seed
 
 
 class TestGaussian(TestCase):
     def setup_method(self, method):
-        if method.__name__ .endswith("prob"):
-            global_seed(314159)
-
         self.mech = Gaussian
 
     def teardown_method(self, method):
@@ -54,7 +50,7 @@ def test_non_numeric(self):
             mech.randomise("Hello")
 
     def test_zero_median_prob(self):
-        mech = self.mech(epsilon=0.75, delta=0.1, sensitivity=1)
+        mech = self.mech(epsilon=0.75, delta=0.1, sensitivity=1, random_state=0)
         vals = []
 
         for i in range(20000):
@@ -66,7 +62,7 @@ def test_zero_median_prob(self):
     def test_neighbors_prob(self):
         epsilon = 1
         runs = 10000
-        mech = self.mech(epsilon=0.5, delta=0.1, sensitivity=1)
+        mech = self.mech(epsilon=0.5, delta=0.1, sensitivity=1, random_state=0)
         count = [0, 0]
 
         for i in range(runs):
@@ -81,6 +77,18 @@ def test_neighbors_prob(self):
         self.assertGreater(count[0], count[1])
         self.assertLessEqual(count[0] / runs, np.exp(epsilon) * count[1] / runs + 0.1)
 
+    def test_random_state(self):
+        mech1 = self.mech(epsilon=1, delta=1e-5, sensitivity=1, random_state=42)
+        mech2 = self.mech(epsilon=1, delta=1e-5, sensitivity=1, random_state=42)
+        self.assertEqual([mech1.randomise(0) for _ in range(100)], [mech2.randomise(0) for _ in range(100)])
+
+        self.assertNotEqual([mech1.randomise(0)] * 100, [mech1.randomise(0) for _ in range(100)])
+
+        rng = np.random.RandomState(0)
+        mech1 = self.mech(epsilon=1, delta=1e-5, sensitivity=1, random_state=rng)
+        mech2 = self.mech(epsilon=1, delta=1e-5, sensitivity=1, random_state=rng)
+        self.assertNotEqual([mech1.randomise(0) for _ in range(100)], [mech2.randomise(0) for _ in range(100)])
+
     def test_repr(self):
         repr_ = repr(self.mech(epsilon=0.5, delta=0.1, sensitivity=1))
         self.assertIn(".Gaussian(", repr_)
diff --git a/tests/mechanisms/test_GaussianAnalytic.py b/tests/mechanisms/test_GaussianAnalytic.py
index 0a79eeb1..80f548df 100644
--- a/tests/mechanisms/test_GaussianAnalytic.py
+++ b/tests/mechanisms/test_GaussianAnalytic.py
@@ -2,14 +2,10 @@
 from unittest import TestCase
 
 from diffprivlib.mechanisms import GaussianAnalytic
-from diffprivlib.utils import global_seed
 
 
 class TestGaussianAnalytic(TestCase):
     def setup_method(self, method):
-        if method.__name__ .endswith("prob"):
-            global_seed(314159)
-
         self.mech = GaussianAnalytic
 
     def teardown_method(self, method):
@@ -60,7 +56,7 @@ def test_simple(self):
         self.assertIsNotNone(mech.randomise(3))
 
     def test_zero_median_prob(self):
-        mech = self.mech(epsilon=0.5, delta=0.1, sensitivity=1)
+        mech = self.mech(epsilon=0.5, delta=0.1, sensitivity=1, random_state=0)
         vals = []
 
         for i in range(20000):
@@ -72,7 +68,7 @@ def test_zero_median_prob(self):
     def test_neighbors_prob(self):
         epsilon = 1
         runs = 10000
-        mech = self.mech(epsilon=epsilon, delta=0.1, sensitivity=1)
+        mech = self.mech(epsilon=epsilon, delta=0.1, sensitivity=1, random_state=0)
         count = [0, 0]
 
         for i in range(runs):
@@ -87,6 +83,18 @@ def test_neighbors_prob(self):
         self.assertGreater(count[0], count[1])
         self.assertLessEqual(count[0] / runs, np.exp(epsilon) * count[1] / runs + 0.1)
 
+    def test_random_state(self):
+        mech1 = self.mech(epsilon=1, delta=1e-5, sensitivity=1, random_state=42)
+        mech2 = self.mech(epsilon=1, delta=1e-5, sensitivity=1, random_state=42)
+        self.assertEqual([mech1.randomise(0) for _ in range(100)], [mech2.randomise(0) for _ in range(100)])
+
+        self.assertNotEqual([mech1.randomise(0)] * 100, [mech1.randomise(0) for _ in range(100)])
+
+        rng = np.random.RandomState(0)
+        mech1 = self.mech(epsilon=1, delta=1e-5, sensitivity=1, random_state=rng)
+        mech2 = self.mech(epsilon=1, delta=1e-5, sensitivity=1, random_state=rng)
+        self.assertNotEqual([mech1.randomise(0) for _ in range(100)], [mech2.randomise(0) for _ in range(100)])
+
     def test_repr(self):
         repr_ = repr(self.mech(epsilon=0.5, delta=0.1, sensitivity=1))
         self.assertIn(".GaussianAnalytic(", repr_)
diff --git a/tests/mechanisms/test_GaussianDiscrete.py b/tests/mechanisms/test_GaussianDiscrete.py
index f5b9532e..179b1f7f 100644
--- a/tests/mechanisms/test_GaussianDiscrete.py
+++ b/tests/mechanisms/test_GaussianDiscrete.py
@@ -2,14 +2,10 @@
 from unittest import TestCase
 
 from diffprivlib.mechanisms import GaussianDiscrete
-from diffprivlib.utils import global_seed
 
 
 class TestGaussianDiscrete(TestCase):
     def setup_method(self, method):
-        if method.__name__ .endswith("prob"):
-            global_seed(314159)
-
         self.mech = GaussianDiscrete
 
     def teardown_method(self, method):
@@ -75,7 +71,7 @@ def test_simple(self):
         self.assertIsNotNone(mech.randomise(3))
 
     def test_zero_median_prob(self):
-        mech = self.mech(epsilon=0.5, delta=0.1)
+        mech = self.mech(epsilon=0.5, delta=0.1, random_state=0)
         vals = []
 
         for i in range(10000):
@@ -85,7 +81,7 @@ def test_zero_median_prob(self):
         self.assertAlmostEqual(np.abs(median), 1, delta=0.1)
 
     def test_zero_median_sens_prob(self):
-        mech = self.mech(epsilon=0.5, delta=0.1, sensitivity=4)
+        mech = self.mech(epsilon=0.5, delta=0.1, sensitivity=4, random_state=0)
         vals = []
 
         for i in range(10000):
@@ -97,7 +93,7 @@ def test_zero_median_sens_prob(self):
     def test_neighbors_prob(self):
         epsilon = 1
         runs = 1000
-        mech = self.mech(epsilon=epsilon, delta=0.1)
+        mech = self.mech(epsilon=epsilon, delta=0.1, random_state=0)
         count = [0, 0]
 
         for i in range(runs):
@@ -123,6 +119,18 @@ def test_sanity_scale(self):
 
         self.assertAlmostEqual(2, scale3, places=5)
 
+    def test_random_state(self):
+        mech1 = self.mech(epsilon=1, delta=1e-5, sensitivity=1, random_state=42)
+        mech2 = self.mech(epsilon=1, delta=1e-5, sensitivity=1, random_state=42)
+        self.assertEqual([mech1.randomise(0) for _ in range(10)], [mech2.randomise(0) for _ in range(10)])
+
+        self.assertNotEqual([mech1.randomise(0)] * 100, [mech1.randomise(0) for _ in range(100)])
+
+        rng = np.random.RandomState(0)
+        mech1 = self.mech(epsilon=1, delta=1e-5, sensitivity=1, random_state=rng)
+        mech2 = self.mech(epsilon=1, delta=1e-5, sensitivity=1, random_state=rng)
+        self.assertNotEqual([mech1.randomise(0) for _ in range(100)], [mech2.randomise(0) for _ in range(100)])
+
     def test_bias(self):
         self.assertEqual(0, self.mech(epsilon=1, delta=0.5).bias(0))
 
diff --git a/tests/mechanisms/test_Geometric.py b/tests/mechanisms/test_Geometric.py
index fa2387b0..6aacfdad 100644
--- a/tests/mechanisms/test_Geometric.py
+++ b/tests/mechanisms/test_Geometric.py
@@ -2,14 +2,10 @@
 from unittest import TestCase
 
 from diffprivlib.mechanisms import Geometric
-from diffprivlib.utils import global_seed
 
 
 class TestGeometric(TestCase):
     def setup_method(self, method):
-        if method.__name__ .endswith("prob"):
-            global_seed(314159)
-
         self.mech = Geometric
 
     def teardown_method(self, method):
@@ -87,7 +83,7 @@ def test_zero_median(self):
     def test_neighbors_prob(self):
         epsilon = 1
         runs = 10000
-        mech = self.mech(epsilon=epsilon, sensitivity=1)
+        mech = self.mech(epsilon=epsilon, sensitivity=1, random_state=0)
         count = [0, 0]
 
         for i in range(runs):
@@ -102,6 +98,18 @@ def test_neighbors_prob(self):
         self.assertGreater(count[0], count[1])
         self.assertLessEqual(count[0] / runs, np.exp(epsilon) * count[1] / runs + 0.1)
 
+    def test_random_state(self):
+        mech1 = self.mech(epsilon=1, sensitivity=1, random_state=42)
+        mech2 = self.mech(epsilon=1, sensitivity=1, random_state=42)
+        self.assertEqual([mech1.randomise(0) for _ in range(100)], [mech2.randomise(0) for _ in range(100)])
+
+        self.assertNotEqual([mech1.randomise(0)] * 100, [mech1.randomise(0) for _ in range(100)])
+
+        rng = np.random.RandomState(0)
+        mech1 = self.mech(epsilon=1, sensitivity=1, random_state=rng)
+        mech2 = self.mech(epsilon=1, sensitivity=1, random_state=rng)
+        self.assertNotEqual([mech1.randomise(0) for _ in range(100)], [mech2.randomise(0) for _ in range(100)])
+
     def test_repr(self):
         repr_ = repr(self.mech(epsilon=1, sensitivity=1))
         self.assertIn(".Geometric(", repr_)
diff --git a/tests/mechanisms/test_GeometricFolded.py b/tests/mechanisms/test_GeometricFolded.py
index f2bb5e42..15a11579 100644
--- a/tests/mechanisms/test_GeometricFolded.py
+++ b/tests/mechanisms/test_GeometricFolded.py
@@ -2,14 +2,10 @@
 from unittest import TestCase
 
 from diffprivlib.mechanisms import GeometricFolded
-from diffprivlib.utils import global_seed
 
 
 class TestGeometricFolded(TestCase):
     def setup_method(self, method):
-        if method.__name__ .endswith("prob"):
-            global_seed(314159)
-
         self.mech = GeometricFolded
 
     def teardown_method(self, method):
@@ -85,7 +81,7 @@ def test_non_integer(self):
             mech.randomise(1.0)
 
     def test_zero_median_prob(self):
-        mech = self.mech(epsilon=1, sensitivity=1, lower=0, upper=4)
+        mech = self.mech(epsilon=1, sensitivity=1, lower=0, upper=4, random_state=0)
         vals = []
 
         for i in range(1000):
@@ -97,7 +93,7 @@ def test_zero_median_prob(self):
     def test_neighbors_prob(self):
         epsilon = np.log(2)
         runs = 1000
-        mech = self.mech(epsilon=epsilon, sensitivity=1, lower=0, upper=4)
+        mech = self.mech(epsilon=epsilon, sensitivity=1, lower=0, upper=4, random_state=0)
         count = [0, 0]
 
         for i in range(runs):
@@ -112,6 +108,18 @@ def test_neighbors_prob(self):
         self.assertGreater(count[0], count[1])
         self.assertLessEqual(count[0], np.exp(epsilon) * count[1] + 0.15 * runs)
 
+    def test_random_state(self):
+        mech1 = self.mech(epsilon=1, sensitivity=1, lower=0, upper=4, random_state=42)
+        mech2 = self.mech(epsilon=1, sensitivity=1, lower=0, upper=4, random_state=42)
+        self.assertEqual([mech1.randomise(0) for _ in range(100)], [mech2.randomise(0) for _ in range(100)])
+
+        self.assertNotEqual([mech1.randomise(0)] * 100, [mech1.randomise(0) for _ in range(100)])
+
+        rng = np.random.RandomState(0)
+        mech1 = self.mech(epsilon=1, sensitivity=1, lower=0, upper=4, random_state=rng)
+        mech2 = self.mech(epsilon=1, sensitivity=1, lower=0, upper=4, random_state=rng)
+        self.assertNotEqual([mech1.randomise(0) for _ in range(100)], [mech2.randomise(0) for _ in range(100)])
+
     def test_repr(self):
         repr_ = repr(self.mech(epsilon=1, sensitivity=1, lower=0, upper=4))
         self.assertIn(".GeometricFolded(", repr_)
diff --git a/tests/mechanisms/test_GeometricTruncated.py b/tests/mechanisms/test_GeometricTruncated.py
index 75aca57c..c702a99c 100644
--- a/tests/mechanisms/test_GeometricTruncated.py
+++ b/tests/mechanisms/test_GeometricTruncated.py
@@ -2,14 +2,10 @@
 from unittest import TestCase
 
 from diffprivlib.mechanisms import GeometricTruncated
-from diffprivlib.utils import global_seed
 
 
 class TestGeometricTruncated(TestCase):
     def setup_method(self, method):
-        if method.__name__ .endswith("prob"):
-            global_seed(314159)
-
         self.mech = GeometricTruncated
 
     def teardown_method(self, method):
@@ -83,7 +79,7 @@ def test_non_integer(self):
             mech.randomise(1.0)
 
     def test_zero_median_prob(self):
-        mech = self.mech(epsilon=1, sensitivity=1, lower=0, upper=4)
+        mech = self.mech(epsilon=1, sensitivity=1, lower=0, upper=4, random_state=0)
         vals = []
 
         for i in range(1000):
@@ -94,8 +90,8 @@ def test_zero_median_prob(self):
 
     def test_neighbors_prob(self):
         epsilon = 1
-        runs = 10000
-        mech = self.mech(epsilon=1, sensitivity=1, lower=0, upper=4)
+        runs = 5000
+        mech = self.mech(epsilon=1, sensitivity=1, lower=0, upper=4, random_state=0)
         count = [0, 0]
 
         for i in range(runs):
@@ -110,6 +106,18 @@ def test_neighbors_prob(self):
         self.assertGreater(count[0], count[1])
         self.assertLessEqual(count[0] / runs, np.exp(epsilon) * count[1] / runs + 0.1)
 
+    def test_random_state(self):
+        mech1 = self.mech(epsilon=1, sensitivity=1, lower=0, upper=4, random_state=42)
+        mech2 = self.mech(epsilon=1, sensitivity=1, lower=0, upper=4, random_state=42)
+        self.assertEqual([mech1.randomise(0) for _ in range(100)], [mech2.randomise(0) for _ in range(100)])
+
+        self.assertNotEqual([mech1.randomise(0)] * 100, [mech1.randomise(0) for _ in range(100)])
+
+        rng = np.random.RandomState(0)
+        mech1 = self.mech(epsilon=1, sensitivity=1, lower=0, upper=4, random_state=rng)
+        mech2 = self.mech(epsilon=1, sensitivity=1, lower=0, upper=4, random_state=rng)
+        self.assertNotEqual([mech1.randomise(0) for _ in range(100)], [mech2.randomise(0) for _ in range(100)])
+
     def test_repr(self):
         repr_ = repr(self.mech(epsilon=1, sensitivity=1, lower=0, upper=4))
         self.assertIn(".GeometricTruncated(", repr_)
diff --git a/tests/mechanisms/test_Laplace.py b/tests/mechanisms/test_Laplace.py
index 1c577169..0c60a4cb 100644
--- a/tests/mechanisms/test_Laplace.py
+++ b/tests/mechanisms/test_Laplace.py
@@ -2,14 +2,10 @@
 import numpy as np
 
 from diffprivlib.mechanisms import Laplace
-from diffprivlib.utils import global_seed
 
 
 class TestLaplace(TestCase):
     def setup_method(self, method):
-        if method.__name__ .endswith("prob"):
-            global_seed(314159)
-
         self.mech = Laplace
 
     def teardown_method(self, method):
@@ -73,7 +69,7 @@ def test_non_numeric(self):
             mech.randomise("Hello")
 
     def test_zero_median_prob(self):
-        mech = self.mech(epsilon=1, delta=0, sensitivity=1)
+        mech = self.mech(epsilon=1, delta=0, sensitivity=1, random_state=0)
         vals = []
 
         for i in range(10000):
@@ -85,7 +81,7 @@ def test_zero_median_prob(self):
     def test_neighbours_prob(self):
         epsilon = 1
         runs = 10000
-        mech = self.mech(epsilon=epsilon, delta=0, sensitivity=1)
+        mech = self.mech(epsilon=epsilon, delta=0, sensitivity=1, random_state=0)
         count = [0, 0]
 
         for i in range(runs):
@@ -100,6 +96,18 @@ def test_neighbours_prob(self):
         self.assertGreater(count[0], count[1])
         self.assertLessEqual(count[0] / runs, np.exp(epsilon) * count[1] / runs + 0.1)
 
+    def test_random_state(self):
+        mech1 = self.mech(epsilon=1, sensitivity=1, random_state=42)
+        mech2 = self.mech(epsilon=1, sensitivity=1, random_state=42)
+        self.assertEqual([mech1.randomise(0) for _ in range(100)], [mech2.randomise(0) for _ in range(100)])
+
+        self.assertNotEqual([mech1.randomise(0)] * 100, [mech1.randomise(0) for _ in range(100)])
+
+        rng = np.random.RandomState(0)
+        mech1 = self.mech(epsilon=1, sensitivity=1, random_state=rng)
+        mech2 = self.mech(epsilon=1, sensitivity=1, random_state=rng)
+        self.assertNotEqual([mech1.randomise(0) for _ in range(100)], [mech2.randomise(0) for _ in range(100)])
+
     def test_bias(self):
         self.assertEqual(0.0, self.mech(epsilon=1, delta=0, sensitivity=1).bias(0))
 
diff --git a/tests/mechanisms/test_LaplaceBoundedDomain.py b/tests/mechanisms/test_LaplaceBoundedDomain.py
index 69ad595a..5fda546e 100644
--- a/tests/mechanisms/test_LaplaceBoundedDomain.py
+++ b/tests/mechanisms/test_LaplaceBoundedDomain.py
@@ -6,14 +6,10 @@
 import pytest
 
 from diffprivlib.mechanisms import LaplaceBoundedDomain
-from diffprivlib.utils import global_seed
 
 
 class TestLaplaceBoundedDomain(TestCase):
     def setup_method(self, method):
-        if method.__name__ .endswith("prob"):
-            global_seed(314159)
-
         self.mech = LaplaceBoundedDomain
 
     def teardown_method(self, method):
@@ -56,7 +52,7 @@ def test_non_numeric(self):
             mech.randomise("Hello")
 
     def test_zero_median_prob(self):
-        mech = self.mech(epsilon=1, delta=0, sensitivity=1, lower=0, upper=1)
+        mech = self.mech(epsilon=1, delta=0, sensitivity=1, lower=0, upper=1, random_state=0)
         vals = []
 
         for i in range(10000):
@@ -67,21 +63,14 @@ def test_zero_median_prob(self):
 
     def test_neighbors_prob(self):
         epsilon = 1
-        runs = 10000
-        mech = self.mech(epsilon=1, delta=0, sensitivity=1, lower=0, upper=1)
-        count = [0, 0]
-
-        for i in range(runs):
-            val0 = mech.randomise(0)
-            if val0 <= 0.5:
-                count[0] += 1
+        runs = 1000
+        mech = self.mech(epsilon=1, delta=0, sensitivity=1, lower=0, upper=1, random_state=0)
 
-            val1 = mech.randomise(1)
-            if val1 <= 0.5:
-                count[1] += 1
+        count0 = (np.array([mech.randomise(0) for _ in range(runs)]) <= 0.5).sum()
+        count1 = (np.array([mech.randomise(1) for _ in range(runs)]) <= 0.5).sum()
 
-        self.assertGreater(count[0], count[1])
-        self.assertLessEqual(count[0] / runs, np.exp(epsilon) * count[1] / runs + 0.1)
+        self.assertGreater(count0, count1)
+        self.assertLessEqual(count0 / runs, np.exp(epsilon) * count1 / runs + 0.1)
 
     def test_within_bounds(self):
         mech = self.mech(epsilon=1, delta=0, sensitivity=1, lower=0, upper=1)
@@ -103,6 +92,18 @@ def test_semi_inf_domain_inf_epsilon(self):
 
         self.assertFalse(w, "Warning thrown for LaplaceBoundedDomain")
 
+    def test_random_state(self):
+        mech1 = self.mech(epsilon=1, sensitivity=1, lower=0, upper=4, random_state=42)
+        mech2 = self.mech(epsilon=1, sensitivity=1, lower=0, upper=4, random_state=42)
+        self.assertEqual([mech1.randomise(0) for _ in range(100)], [mech2.randomise(0) for _ in range(100)])
+
+        self.assertNotEqual([mech1.randomise(0)] * 100, [mech1.randomise(0) for _ in range(100)])
+
+        rng = np.random.RandomState(0)
+        mech1 = self.mech(epsilon=1, sensitivity=1, lower=0, upper=4, random_state=rng)
+        mech2 = self.mech(epsilon=1, sensitivity=1, lower=0, upper=4, random_state=rng)
+        self.assertNotEqual([mech1.randomise(0) for _ in range(100)], [mech2.randomise(0) for _ in range(100)])
+
     def test_repr(self):
         repr_ = repr(self.mech(epsilon=1, delta=0, sensitivity=1, lower=0, upper=1))
         self.assertIn(".LaplaceBoundedDomain(", repr_)
diff --git a/tests/mechanisms/test_LaplaceBoundedNoise.py b/tests/mechanisms/test_LaplaceBoundedNoise.py
index 57db7222..cc70b3ac 100644
--- a/tests/mechanisms/test_LaplaceBoundedNoise.py
+++ b/tests/mechanisms/test_LaplaceBoundedNoise.py
@@ -2,14 +2,10 @@
 from unittest import TestCase
 
 from diffprivlib.mechanisms import LaplaceBoundedNoise
-from diffprivlib.utils import global_seed
 
 
 class TestLaplaceBoundedNoise(TestCase):
     def setup_method(self, method):
-        if method.__name__ .endswith("prob"):
-            global_seed(314159)
-
         self.mech = LaplaceBoundedNoise
 
     def teardown_method(self, method):
@@ -65,7 +61,7 @@ def test_non_nan(self):
         self.assertTrue(np.isnan(mech.randomise(np.nan)))
 
     def test_zero_median_prob(self):
-        mech = self.mech(epsilon=1, delta=0.1, sensitivity=1)
+        mech = self.mech(epsilon=1, delta=0.1, sensitivity=1, random_state=0)
         vals = []
 
         for i in range(10000):
@@ -75,22 +71,15 @@ def test_zero_median_prob(self):
         self.assertAlmostEqual(np.abs(median), 0.0, delta=0.1)
 
     def test_neighbors_prob(self):
-        runs = 10000
+        runs = 1000
         delta = 0.1
-        mech = self.mech(epsilon=1, delta=delta, sensitivity=1)
-        count = [0, 0]
-
-        for i in range(runs):
-            val0 = mech.randomise(0)
-            if val0 <= 1 - mech._noise_bound:
-                count[0] += 1
+        mech = self.mech(epsilon=1, delta=delta, sensitivity=1, random_state=2)
 
-            val1 = mech.randomise(1)
-            if val1 >= mech._noise_bound:
-                count[1] += 1
+        count0 = (np.array([mech.randomise(0) for _ in range(runs)]) <= 1 - mech._noise_bound).sum()
+        count1 = (np.array([mech.randomise(1) for _ in range(runs)]) >= mech._noise_bound).sum()
 
-        self.assertAlmostEqual(count[0] / runs, delta, delta=delta/10)
-        self.assertAlmostEqual(count[1] / runs, delta, delta=delta/10)
+        self.assertAlmostEqual(count0 / runs, delta, delta=delta / 10)
+        self.assertAlmostEqual(count1 / runs, delta, delta=delta / 10)
 
     def test_within_bounds(self):
         mech = self.mech(epsilon=1, delta=0.1, sensitivity=1)
@@ -104,6 +93,18 @@ def test_within_bounds(self):
         self.assertTrue(np.all(vals >= -mech._noise_bound))
         self.assertTrue(np.all(vals <= mech._noise_bound))
 
+    def test_random_state(self):
+        mech1 = self.mech(epsilon=1, delta=1e-3, sensitivity=1, random_state=42)
+        mech2 = self.mech(epsilon=1, delta=1e-3, sensitivity=1, random_state=42)
+        self.assertEqual([mech1.randomise(0) for _ in range(100)], [mech2.randomise(0) for _ in range(100)])
+
+        self.assertNotEqual([mech1.randomise(0)] * 100, [mech1.randomise(0) for _ in range(100)])
+
+        rng = np.random.RandomState(0)
+        mech1 = self.mech(epsilon=1, delta=1e-3, sensitivity=1, random_state=rng)
+        mech2 = self.mech(epsilon=1, delta=1e-3, sensitivity=1, random_state=rng)
+        self.assertNotEqual([mech1.randomise(0) for _ in range(100)], [mech2.randomise(0) for _ in range(100)])
+
     def test_repr(self):
         repr_ = repr(self.mech(epsilon=1, delta=0.1, sensitivity=1))
         self.assertIn(".LaplaceBoundedNoise(", repr_)
diff --git a/tests/mechanisms/test_LaplaceFolded.py b/tests/mechanisms/test_LaplaceFolded.py
index 4ccc992b..e18da3ae 100644
--- a/tests/mechanisms/test_LaplaceFolded.py
+++ b/tests/mechanisms/test_LaplaceFolded.py
@@ -2,14 +2,10 @@
 from unittest import TestCase
 
 from diffprivlib.mechanisms import LaplaceFolded
-from diffprivlib.utils import global_seed
 
 
 class TestLaplaceFolded(TestCase):
     def setup_method(self, method):
-        if method.__name__ .endswith("prob"):
-            global_seed(314159)
-
         self.mech = LaplaceFolded
 
     def teardown_method(self, method):
@@ -52,7 +48,7 @@ def test_non_numeric(self):
             mech.randomise("Hello")
 
     def test_zero_median_prob(self):
-        mech = self.mech(epsilon=1, delta=0, sensitivity=1, lower=0, upper=1)
+        mech = self.mech(epsilon=1, delta=0, sensitivity=1, lower=0, upper=1, random_state=0)
         vals = []
 
         for i in range(10000):
@@ -64,7 +60,7 @@ def test_zero_median_prob(self):
     def test_neighbors_prob(self):
         epsilon = 1
         runs = 10000
-        mech = self.mech(epsilon=1, delta=0, sensitivity=1, lower=0, upper=1)
+        mech = self.mech(epsilon=1, delta=0, sensitivity=1, lower=0, upper=1, random_state=0)
         count = [0, 0]
 
         for i in range(runs):
@@ -91,6 +87,18 @@ def test_within_bounds(self):
         self.assertTrue(np.all(vals >= 0))
         self.assertTrue(np.all(vals <= 1))
 
+    def test_random_state(self):
+        mech1 = self.mech(epsilon=1, sensitivity=1, lower=0, upper=4, random_state=42)
+        mech2 = self.mech(epsilon=1, sensitivity=1, lower=0, upper=4, random_state=42)
+        self.assertEqual([mech1.randomise(0) for _ in range(100)], [mech2.randomise(0) for _ in range(100)])
+
+        self.assertNotEqual([mech1.randomise(0)] * 100, [mech1.randomise(0) for _ in range(100)])
+
+        rng = np.random.RandomState(0)
+        mech1 = self.mech(epsilon=1, sensitivity=1, lower=0, upper=4, random_state=rng)
+        mech2 = self.mech(epsilon=1, sensitivity=1, lower=0, upper=4, random_state=rng)
+        self.assertNotEqual([mech1.randomise(0) for _ in range(100)], [mech2.randomise(0) for _ in range(100)])
+
     def test_repr(self):
         repr_ = repr(self.mech(epsilon=1, delta=0, sensitivity=1, lower=0, upper=1))
         self.assertIn(".LaplaceFolded(", repr_)
diff --git a/tests/mechanisms/test_LaplaceTruncated.py b/tests/mechanisms/test_LaplaceTruncated.py
index 7cc76e44..babb527f 100644
--- a/tests/mechanisms/test_LaplaceTruncated.py
+++ b/tests/mechanisms/test_LaplaceTruncated.py
@@ -2,14 +2,10 @@
 from unittest import TestCase
 
 from diffprivlib.mechanisms import LaplaceTruncated
-from diffprivlib.utils import global_seed
 
 
 class TestLaplaceTruncated(TestCase):
     def setup_method(self, method):
-        if method.__name__ .endswith("prob"):
-            global_seed(314159)
-
         self.mech = LaplaceTruncated
 
     def teardown_method(self, method):
@@ -52,7 +48,7 @@ def test_non_numeric(self):
             mech.randomise("Hello")
 
     def test_zero_median_prob(self):
-        mech = self.mech(epsilon=1, delta=0, sensitivity=1, lower=0, upper=1)
+        mech = self.mech(epsilon=1, delta=0, sensitivity=1, lower=0, upper=1, random_state=0)
         vals = []
 
         for i in range(10000):
@@ -64,7 +60,7 @@ def test_zero_median_prob(self):
     def test_neighbors_prob(self):
         epsilon = 1
         runs = 10000
-        mech = self.mech(epsilon=1, delta=0, sensitivity=1, lower=0, upper=1)
+        mech = self.mech(epsilon=1, delta=0, sensitivity=1, lower=0, upper=1, random_state=0)
         count = [0, 0]
 
         for i in range(runs):
@@ -91,6 +87,18 @@ def test_within_bounds(self):
         self.assertTrue(np.all(vals >= 0))
         self.assertTrue(np.all(vals <= 1))
 
+    def test_random_state(self):
+        mech1 = self.mech(epsilon=1, sensitivity=1, lower=0, upper=4, random_state=42)
+        mech2 = self.mech(epsilon=1, sensitivity=1, lower=0, upper=4, random_state=42)
+        self.assertEqual([mech1.randomise(0) for _ in range(100)], [mech2.randomise(0) for _ in range(100)])
+
+        self.assertNotEqual([mech1.randomise(0)] * 100, [mech1.randomise(0) for _ in range(100)])
+
+        rng = np.random.RandomState(0)
+        mech1 = self.mech(epsilon=1, sensitivity=1, lower=0, upper=4, random_state=rng)
+        mech2 = self.mech(epsilon=1, sensitivity=1, lower=0, upper=4, random_state=rng)
+        self.assertNotEqual([mech1.randomise(0) for _ in range(100)], [mech2.randomise(0) for _ in range(100)])
+
     def test_repr(self):
         repr_ = repr(self.mech(epsilon=1, delta=0, sensitivity=1, lower=0, upper=1))
         self.assertIn(".LaplaceTruncated(", repr_)
diff --git a/tests/mechanisms/test_PermuteAndFlip.py b/tests/mechanisms/test_PermuteAndFlip.py
index 44013c62..937d2caa 100644
--- a/tests/mechanisms/test_PermuteAndFlip.py
+++ b/tests/mechanisms/test_PermuteAndFlip.py
@@ -2,6 +2,7 @@
 from unittest import TestCase
 
 from diffprivlib.mechanisms import PermuteAndFlip
+from diffprivlib.utils import check_random_state
 
 
 class TestPermuteAndFlip(TestCase):
@@ -162,8 +163,9 @@ def test_inf_utility_measure(self):
     def test_distrib_prob(self):
         epsilon = np.log(2)
         runs = 20000
-        mech1 = self.mech(epsilon=epsilon, utility=[2, 1, 0], sensitivity=1, monotonic=False)
-        mech2 = self.mech(epsilon=epsilon, utility=[2, 1, 1], sensitivity=1, monotonic=False)
+        rng = np.random.RandomState(42)
+        mech1 = self.mech(epsilon=epsilon, utility=[2, 1, 0], sensitivity=1, monotonic=False, random_state=rng)
+        mech2 = self.mech(epsilon=epsilon, utility=[2, 1, 1], sensitivity=1, monotonic=False, random_state=rng)
         counts = np.zeros((2, 3))
 
         for i in range(runs):
@@ -176,9 +178,10 @@ def test_distrib_prob(self):
 
     def test_monotonic_distrib(self):
         epsilon = np.log(2)
-        runs = 40000
-        mech1 = self.mech(epsilon=epsilon, utility=[2, 1, 0], sensitivity=1, monotonic=True)
-        mech2 = self.mech(epsilon=epsilon, utility=[2, 1, 1], sensitivity=1, monotonic=True)
+        runs = 1000
+        rng = check_random_state(0)
+        mech1 = self.mech(epsilon=epsilon, utility=[2, 1, 0], sensitivity=1, monotonic=True, random_state=rng)
+        mech2 = self.mech(epsilon=epsilon, utility=[2, 1, 1], sensitivity=1, monotonic=True, random_state=rng)
         counts = np.zeros((2, 3))
 
         for i in range(runs):
@@ -186,9 +189,18 @@ def test_monotonic_distrib(self):
             counts[1, mech2.randomise()] += 1
 
         for vec in counts.T:
-            # print(vec.max() / vec.min())
             self.assertLessEqual(vec.max() / vec.min(), np.exp(epsilon) + 0.1)
 
+    def test_random_state(self):
+        mech1 = self.mech(epsilon=1, utility=[2, 1, 0], sensitivity=1, random_state=42)
+        mech2 = self.mech(epsilon=1, utility=[2, 1, 0], sensitivity=1, random_state=42)
+        self.assertEqual([mech1.randomise() for _ in range(100)], [mech2.randomise() for _ in range(100)])
+
+        self.assertNotEqual([mech1.randomise()] * 100, [mech1.randomise() for _ in range(100)])
+
+        mech2 = self.mech(epsilon=1, utility=[2, 1, 0], sensitivity=1, random_state=np.random.RandomState(0))
+        self.assertNotEqual([mech1.randomise() for _ in range(100)], [mech2.randomise() for _ in range(100)])
+
     def test_repr(self):
         repr_ = repr(self.mech(epsilon=1, utility=[1], sensitivity=1))
         self.assertIn(".PermuteAndFlip(", repr_)
diff --git a/tests/mechanisms/test_Snapping.py b/tests/mechanisms/test_Snapping.py
index d289dc25..87d6e382 100644
--- a/tests/mechanisms/test_Snapping.py
+++ b/tests/mechanisms/test_Snapping.py
@@ -5,14 +5,10 @@
 from unittest import TestCase
 
 from diffprivlib.mechanisms import Snapping
-from diffprivlib.utils import global_seed
 
 
 class TestSnapping(TestCase):
     def setup_method(self, method):
-        if method.__name__ .endswith("prob"):
-            global_seed(314159)
-
         self.mech = Snapping
 
     def teardown_method(self, method):
@@ -60,7 +56,7 @@ def test_epsilon(self):
 
     def test_neg_effective_epsilon(self):
         with self.assertRaises(ValueError):
-            self.mech(epsilon=np.nextafter((2*np.finfo(float).epsneg), 0), sensitivity=1, lower=0, upper=1000)
+            self.mech(epsilon=np.nextafter((2 * np.finfo(float).epsneg), 0), sensitivity=1, lower=0, upper=1000)
 
     def test_zero_effective_epsilon(self):
         with self.assertRaises(ValueError):
@@ -70,20 +66,6 @@ def test_immediately_above_zero_effective_epsilon(self):
         mech = self.mech(epsilon=np.nextafter(2 * np.finfo(float).epsneg, math.inf), sensitivity=1, lower=0, upper=1000)
         self.assertIsNotNone(mech.randomise(1))
 
-    def test_repr(self):
-        repr_ = repr(self.mech(epsilon=1.0, sensitivity=1, lower=0, upper=1000))
-        self.assertIn(".Snapping(", repr_)
-
-    def test_bias(self):
-        mech = self.mech(epsilon=1.0, sensitivity=1, lower=0, upper=1000)
-        with self.assertRaises(NotImplementedError):
-            mech.bias(1)
-
-    def test_variance(self):
-        mech = self.mech(epsilon=1.0, sensitivity=1, lower=0, upper=1000)
-        with self.assertRaises(NotImplementedError):
-            mech.variance(1)
-
     def test_scale_bound_symmetric_sensitivity_1(self):
         mech = self.mech(epsilon=1.0, sensitivity=1, lower=-10, upper=10)
         self.assertAlmostEqual(mech._scale_bound(), 10)
@@ -164,28 +146,14 @@ def test_get_nearest_power_of_2(self):
         self.assertEqual(Snapping._get_nearest_power_of_2(np.nextafter(sys.float_info.min, -math.inf)),
                          sys.float_info.min)
 
-    def test_round_to_nearest_power_of_2_exact(self):
+    def test_round_to_nearest_power_of_2(self):
         mech = self.mech(epsilon=1.0, sensitivity=1, lower=0, upper=1000)
-        self.assertAlmostEqual(mech._round_to_nearest_power_of_2(2.0, 2.0), 2)
 
-    def test_round_to_nearest_power_of_2_below_exact(self):
-        mech = self.mech(epsilon=1.0, sensitivity=1, lower=0, upper=1000)
+        self.assertAlmostEqual(mech._round_to_nearest_power_of_2(2.0, 2.0), 2)
         self.assertAlmostEqual(mech._round_to_nearest_power_of_2(np.nextafter(2.0, -math.inf), 2.0), 2)
-
-    def test_round_to_nearest_power_of_2_above_exact(self):
-        mech = self.mech(epsilon=1.0, sensitivity=1, lower=0, upper=1000)
         self.assertAlmostEqual(mech._round_to_nearest_power_of_2(np.nextafter(2.0, math.inf), 2.0), 2.0)
-
-    def test_round_to_nearest_power_of_2_middle(self):
-        mech = self.mech(epsilon=1.0, sensitivity=1, lower=0, upper=1000)
         self.assertAlmostEqual(mech._round_to_nearest_power_of_2(3.0, 2.0), 4)
-
-    def test_round_to_nearest_power_of_2_below_middle(self):
-        mech = self.mech(epsilon=1.0, sensitivity=1, lower=0, upper=1000)
         self.assertAlmostEqual(mech._round_to_nearest_power_of_2(np.nextafter(3.0, -math.inf), 2.0), 2)
-
-    def test_round_to_nearest_power_of_2_above_middle(self):
-        mech = self.mech(epsilon=1.0, sensitivity=1, lower=0, upper=1000)
         self.assertAlmostEqual(mech._round_to_nearest_power_of_2(np.nextafter(3.0, math.inf), 2.0), 4.0)
 
     def test_non_numeric(self):
@@ -222,18 +190,37 @@ def test_within_bounds(self):
 
     def test_neighbours_prob(self):
         epsilon = 1
-        runs = 10000
-        mech = self.mech(epsilon=epsilon, sensitivity=1, lower=0, upper=1000)
-        count = [0, 0]
+        runs = 5000
+        mech = self.mech(epsilon=epsilon, sensitivity=1, lower=0, upper=1000, random_state=0)
+
+        count0 = (np.array([mech.randomise(0) for _ in range(runs)]) <= 0).sum()
+        count1 = (np.array([mech.randomise(1) for _ in range(runs)]) <= 0).sum()
 
-        for i in range(runs):
-            val0 = mech.randomise(0)
-            if val0 <= 0:
-                count[0] += 1
+        self.assertGreater(count0, count1)
+        self.assertLessEqual(count0 / runs, np.exp(epsilon) * count1 / runs + 0.1)
 
-            val1 = mech.randomise(1)
-            if val1 <= 0:
-                count[1] += 1
+    def test_random_state(self):
+        mech1 = self.mech(epsilon=1, sensitivity=1, lower=0, upper=1000, random_state=42)
+        mech2 = self.mech(epsilon=1, sensitivity=1, lower=0, upper=1000, random_state=42)
+        self.assertEqual([mech1.randomise(0) for _ in range(100)], [mech2.randomise(0) for _ in range(100)])
 
-        self.assertGreater(count[0], count[1])
-        self.assertLessEqual(count[0] / runs, np.exp(epsilon) * count[1] / runs + 0.1)
+        self.assertNotEqual([mech1.randomise(0)] * 100, [mech1.randomise(0) for _ in range(100)])
+
+        rng = np.random.RandomState(0)
+        mech1 = self.mech(epsilon=1, sensitivity=1, lower=0, upper=1000, random_state=rng)
+        mech2 = self.mech(epsilon=1, sensitivity=1, lower=0, upper=1000, random_state=rng)
+        self.assertNotEqual([mech1.randomise(0) for _ in range(100)], [mech2.randomise(0) for _ in range(100)])
+
+    def test_repr(self):
+        repr_ = repr(self.mech(epsilon=1.0, sensitivity=1, lower=0, upper=1000))
+        self.assertIn(".Snapping(", repr_)
+
+    def test_bias(self):
+        mech = self.mech(epsilon=1.0, sensitivity=1, lower=0, upper=1000)
+        with self.assertRaises(NotImplementedError):
+            mech.bias(1)
+
+    def test_variance(self):
+        mech = self.mech(epsilon=1.0, sensitivity=1, lower=0, upper=1000)
+        with self.assertRaises(NotImplementedError):
+            mech.variance(1)
diff --git a/tests/mechanisms/test_Staircase.py b/tests/mechanisms/test_Staircase.py
index 06d53ff5..ccb7f37b 100644
--- a/tests/mechanisms/test_Staircase.py
+++ b/tests/mechanisms/test_Staircase.py
@@ -2,14 +2,10 @@
 from unittest import TestCase
 
 from diffprivlib.mechanisms import Staircase
-from diffprivlib.utils import global_seed
 
 
 class TestStaircase(TestCase):
     def setup_method(self, method):
-        if method.__name__ .endswith("prob"):
-            global_seed(314159)
-
         self.mech = Staircase
 
     def teardown_method(self, method):
@@ -65,7 +61,7 @@ def test_non_numeric(self):
             mech.randomise("Hello")
 
     def test_zero_median_prob(self):
-        mech = self.mech(epsilon=1, sensitivity=1, gamma=0.5)
+        mech = self.mech(epsilon=1, sensitivity=1, gamma=0.5, random_state=0)
         vals = []
 
         for i in range(10000):
@@ -77,7 +73,7 @@ def test_zero_median_prob(self):
     def test_neighbors_prob(self):
         epsilon = 1
         runs = 10000
-        mech = self.mech(epsilon=epsilon, sensitivity=1, gamma=0.5)
+        mech = self.mech(epsilon=epsilon, sensitivity=1, gamma=0.5, random_state=0)
         count = [0, 0]
 
         for i in range(runs):
@@ -93,6 +89,18 @@ def test_neighbors_prob(self):
         self.assertGreater(count[0], count[1])
         self.assertLessEqual(count[0] / runs, np.exp(epsilon) * count[1] / runs + 0.1)
 
+    def test_random_state(self):
+        mech1 = self.mech(epsilon=1, sensitivity=1, gamma=0.5, random_state=42)
+        mech2 = self.mech(epsilon=1, sensitivity=1, gamma=0.5, random_state=42)
+        self.assertEqual([mech1.randomise(0) for _ in range(100)], [mech2.randomise(0) for _ in range(100)])
+
+        self.assertNotEqual([mech1.randomise(0)] * 100, [mech1.randomise(0) for _ in range(100)])
+
+        rng = np.random.RandomState(0)
+        mech1 = self.mech(epsilon=1, sensitivity=1, gamma=0.5, random_state=rng)
+        mech2 = self.mech(epsilon=1, sensitivity=1, gamma=0.5, random_state=rng)
+        self.assertNotEqual([mech1.randomise(0) for _ in range(100)], [mech2.randomise(0) for _ in range(100)])
+
     def test_repr(self):
         repr_ = repr(self.mech(epsilon=1, sensitivity=1, gamma=0.5))
         self.assertIn(".Staircase(", repr_)
diff --git a/tests/mechanisms/test_Uniform.py b/tests/mechanisms/test_Uniform.py
index bf8f36c4..53f11039 100644
--- a/tests/mechanisms/test_Uniform.py
+++ b/tests/mechanisms/test_Uniform.py
@@ -2,14 +2,10 @@
 from unittest import TestCase
 
 from diffprivlib.mechanisms import Uniform
-from diffprivlib.utils import global_seed
 
 
 class TestUniform(TestCase):
     def setup_method(self, method):
-        if method.__name__ .endswith("prob"):
-            global_seed(314159)
-
         self.mech = Uniform
 
     def teardown_method(self, method):
@@ -60,7 +56,7 @@ def test_non_numeric(self):
             mech.randomise("Hello")
 
     def test_zero_median_prob(self):
-        mech = self.mech(delta=0.2, sensitivity=1)
+        mech = self.mech(delta=0.2, sensitivity=1, random_state=0)
         vals = []
 
         for i in range(10000):
@@ -69,6 +65,18 @@ def test_zero_median_prob(self):
         median = float(np.median(vals))
         self.assertAlmostEqual(np.abs(median), 0.0, delta=0.1)
 
+    def test_random_state(self):
+        mech1 = self.mech(delta=0.2, sensitivity=1, random_state=42)
+        mech2 = self.mech(delta=0.2, sensitivity=1, random_state=42)
+        self.assertEqual([mech1.randomise(0) for _ in range(100)], [mech2.randomise(0) for _ in range(100)])
+
+        self.assertNotEqual([mech1.randomise(0)] * 100, [mech1.randomise(0) for _ in range(100)])
+
+        rng = np.random.RandomState(0)
+        mech1 = self.mech(delta=0.2, sensitivity=1, random_state=rng)
+        mech2 = self.mech(delta=0.2, sensitivity=1, random_state=rng)
+        self.assertNotEqual([mech1.randomise(0) for _ in range(100)], [mech2.randomise(0) for _ in range(100)])
+
     def test_repr(self):
         repr_ = repr(self.mech(delta=0.1, sensitivity=1))
         self.assertIn(".Uniform(", repr_)
diff --git a/tests/mechanisms/test_Vector.py b/tests/mechanisms/test_Vector.py
index 54486de2..6ebea67d 100644
--- a/tests/mechanisms/test_Vector.py
+++ b/tests/mechanisms/test_Vector.py
@@ -2,14 +2,10 @@
 from unittest import TestCase
 
 from diffprivlib.mechanisms import Vector
-from diffprivlib.utils import global_seed
 
 
 class TestVector(TestCase):
     def setup_method(self, method):
-        if method.__name__ .endswith("prob"):
-            global_seed(314159)
-
         self.mech = Vector
 
     def teardown_method(self, method):
@@ -109,6 +105,22 @@ def test_different_result(self):
             self.assertNotAlmostEqual(noisy_func(np.ones(3)), old_noisy_func(np.ones(3)))
             # print(noisy_func(np.ones(3)))
 
+    def test_random_state(self):
+        mech1 = self.mech(epsilon=1, dimension=3, function_sensitivity=1, random_state=0)
+        mech2 = self.mech(epsilon=1, dimension=3, function_sensitivity=1, random_state=0)
+        func1, func2 = mech1.randomise(self.func), mech2.randomise(self.func)
+        self.assertTrue(np.allclose(func1(np.ones(3)), func1(np.ones(3))))
+
+        mech1 = self.mech(epsilon=1, dimension=3, function_sensitivity=1, random_state=0)
+        func1, func2 = mech1.randomise(self.func), mech2.randomise(self.func)
+        self.assertFalse(np.allclose(func1(np.ones(3)), func2(np.ones(3))))
+
+        rng = np.random.RandomState(0)
+        mech1 = self.mech(epsilon=1, dimension=3, function_sensitivity=1, random_state=rng)
+        mech2 = self.mech(epsilon=1, dimension=3, function_sensitivity=1, random_state=rng)
+        func1, func2 = mech1.randomise(self.func), mech2.randomise(self.func)
+        self.assertFalse(np.allclose(func1(np.ones(3)), func2(np.ones(3))))
+
     def test_repr(self):
         repr_ = repr(self.mech(epsilon=1, dimension=4, function_sensitivity=1))
         self.assertIn(".Vector(", repr_)
diff --git a/tests/mechanisms/transforms/test_RoundedInteger.py b/tests/mechanisms/transforms/test_RoundedInteger.py
index dce62ee9..b5e688e8 100644
--- a/tests/mechanisms/transforms/test_RoundedInteger.py
+++ b/tests/mechanisms/transforms/test_RoundedInteger.py
@@ -29,8 +29,8 @@ def test_randomise(self):
 
     def test_distrib(self):
         epsilon = np.log(2)
-        runs = 10000
-        mech = RoundedInteger(Laplace(epsilon=epsilon, sensitivity=1))
+        runs = 5000
+        mech = RoundedInteger(Laplace(epsilon=epsilon, sensitivity=1, random_state=0))
         count = [0, 0]
 
         for _ in range(runs):
diff --git a/tests/mechanisms/transforms/test_StringToInt.py b/tests/mechanisms/transforms/test_StringToInt.py
index 541e21a8..6821ac6c 100644
--- a/tests/mechanisms/transforms/test_StringToInt.py
+++ b/tests/mechanisms/transforms/test_StringToInt.py
@@ -3,9 +3,6 @@
 
 from diffprivlib.mechanisms import Geometric
 from diffprivlib.mechanisms.transforms import StringToInt
-from diffprivlib.utils import global_seed
-
-global_seed(3141592653)
 
 
 class TestStringToInt(TestCase):
diff --git a/tests/models/test_DecisionTreeClassifier.py b/tests/models/test_DecisionTreeClassifier.py
index 5faaadf6..f43a4ad9 100644
--- a/tests/models/test_DecisionTreeClassifier.py
+++ b/tests/models/test_DecisionTreeClassifier.py
@@ -8,7 +8,7 @@
 from sklearn.exceptions import NotFittedError
 
 from diffprivlib.models.forest import DecisionTreeClassifier
-from diffprivlib.utils import PrivacyLeakWarning, global_seed, DiffprivlibCompatibilityWarning, BudgetError
+from diffprivlib.utils import PrivacyLeakWarning, DiffprivlibCompatibilityWarning, BudgetError
 
 
 X = [[-2, -1], [-1, -1], [-1, -2], [1, 1], [1, 2], [2, 1]]
@@ -19,9 +19,6 @@
 
 
 class TestDecisionTreeClassifier(TestCase):
-    def setUp(self):
-        global_seed(2718281828)
-
     def test_not_none(self):
         self.assertIsNotNone(DecisionTreeClassifier)
 
@@ -45,10 +42,10 @@ def test_bad_data(self):
         with self.assertRaises(ValueError):
             DecisionTreeClassifier(bounds=([0], [1])).fit([[1], [2]], [[1, 2], [2, 4]])
 
-    def test_simple(self):
+    def test_simple_prob(self):
         X = np.array([[12, 3, 14], [12, 3, 4], [12, 3, 4], [2, 13, 4], [2, 13, 14], [2, 3, 14], [3, 5, 15]] * 3)
         y = np.array([1, 1, 1, 0, 0, 0, 1] * 3)
-        model = DecisionTreeClassifier(epsilon=5, max_depth=3)
+        model = DecisionTreeClassifier(epsilon=5, max_depth=3, random_state=0)
         with self.assertRaises(NotFittedError):
             check_is_fitted(model)
         # when `bounds` is not provided, we should get a privacy leakage warning
@@ -60,23 +57,40 @@ def test_simple(self):
     def test_with_bounds(self):
         X = np.array([[12, 3, 14], [12, 3, 4], [12, 3, 4], [2, 13, 4], [2, 13, 14], [2, 3, 14], [3, 5, 15]] * 3)
         y = np.array([1, 1, 1, 0, 0, 0, 1] * 3)
-        model = DecisionTreeClassifier(epsilon=5, bounds=([2, 3, 4], [12, 13, 14]), classes=[0, 1], max_depth=3)
+        model = DecisionTreeClassifier(epsilon=5, bounds=([2, 3, 4], [12, 13, 14]), classes=[0, 1], max_depth=3,
+                                       random_state=0)
         with self.assertRaises(NotFittedError):
             check_is_fitted(model)
         model.fit(X, y)
         check_is_fitted(model)
         self.assertTrue(model.predict(np.array([[12, 3, 14]])))
 
-    def test_with_non_binary_labels(self):
+    def test_with_non_binary_labels_prob(self):
         X = np.array([[12, 3, 14], [12, 3, 4], [12, 3, 4], [2, 13, 4], [2, 13, 14], [2, 3, 14], [3, 5, 15]] * 3)
         y = np.array([3, 3, 3, 3, 5, 5, 3] * 3)
-        model = DecisionTreeClassifier(epsilon=5, bounds=([2, 3, 4], [12, 13, 14]), classes=[3, 5], max_depth=3)
+        model = DecisionTreeClassifier(epsilon=5, bounds=([2, 3, 4], [12, 13, 14]), classes=[3, 5], max_depth=3,
+                                       random_state=0)
         with self.assertRaises(NotFittedError):
             check_is_fitted(model)
         model.fit(X, y)
         check_is_fitted(model)
         self.assertEqual(model.predict(np.array([[12, 3, 14]])), 3)
 
+    def test_random_state(self):
+        model0 = DecisionTreeClassifier(epsilon=0.01, bounds=(-2, 2), classes=[-1, 1], max_depth=3, random_state=0)
+        model1 = DecisionTreeClassifier(epsilon=0.01, bounds=(-2, 2), classes=[-1, 1], max_depth=3, random_state=2)
+        model0.fit(X, y)
+        model1.fit(X, y)
+        self.assertNotEqual(model0.predict(np.array([[-2, -1]])), model1.predict(np.array([[-2, -1]])))
+
+        model1 = DecisionTreeClassifier(epsilon=0.01, bounds=(-2, 2), classes=[-1, 1], max_depth=3, random_state=0)
+        model1.fit(X, y)
+        self.assertTrue(np.isclose(model0.predict_proba(np.array([[-2, -1]])),
+                                   model1.predict_proba(np.array([[-2, -1]]))).all())
+
+        model1 = DecisionTreeClassifier(epsilon=0.01, bounds=(-2, 2), classes=[-1, 1], max_depth=3, random_state=None)
+        model1.fit(X, y)
+
     def test_sklearn_methods(self):
         depth = 3
         clf = DecisionTreeClassifier(epsilon=5, bounds=(-2, 2), classes=[-1, 1], max_depth=depth)
diff --git a/tests/models/test_FittingTree.py b/tests/models/test_FittingTree.py
index 77fcb45e..49ab2fdc 100644
--- a/tests/models/test_FittingTree.py
+++ b/tests/models/test_FittingTree.py
@@ -4,6 +4,7 @@
 from sklearn.tree._tree import Tree
 
 from diffprivlib.models.forest import _FittingTree
+from diffprivlib.utils import check_random_state
 
 
 class TestFittingTree(TestCase):
@@ -12,10 +13,16 @@ def test_not_none(self):
 
     def test_no_data(self):
         classes = np.array([0, 1, 2])
-        tree = _FittingTree(5, 3, classes, 1, (0, 1))
+        tree = _FittingTree(5, 3, classes, 1, (0, 1), check_random_state(None))
         tree.build()
         tree.values_ = np.zeros(shape=(tree.node_count, 1, len(tree.classes)))
         d = tree.__getstate__()
 
         sktree = Tree(3, np.array([len(classes)]), 1)
         sktree.__setstate__(d)
+
+    def test_fit_before_build(self):
+        tree = _FittingTree(5, 3, [0, 1], 1, (0, 1), check_random_state(None))
+
+        with self.assertRaises(ValueError):
+            tree.fit([[1, 1, 1]], [0])
diff --git a/tests/models/test_GaussianNB.py b/tests/models/test_GaussianNB.py
index 0e49a67c..fd15db97 100644
--- a/tests/models/test_GaussianNB.py
+++ b/tests/models/test_GaussianNB.py
@@ -5,51 +5,44 @@
 from sklearn.model_selection import train_test_split
 
 from diffprivlib.models.naive_bayes import GaussianNB
-from diffprivlib.utils import global_seed, PrivacyLeakWarning, DiffprivlibCompatibilityWarning, BudgetError
+from diffprivlib.utils import PrivacyLeakWarning, DiffprivlibCompatibilityWarning, BudgetError, check_random_state
+
+# Data is just 6 separable points in the plane
+X = np.array([[-2, -1], [-1, -1], [-1, -2], [1, 1], [1, 2], [2, 1]])
+y = np.array([1, 1, 1, 2, 2, 2])
 
 
 class TestGaussianNB(TestCase):
     def test_not_none(self):
-        clf = GaussianNB(epsilon=1, bounds=(0, 1))
+        clf = GaussianNB(epsilon=1, bounds=(-2, 2))
         self.assertIsNotNone(clf)
 
     def test_zero_epsilon(self):
-        X = np.random.random((10, 2))
-        y = np.random.randint(2, size=10)
-        clf = GaussianNB(epsilon=0, bounds=(0, 1))
+        clf = GaussianNB(epsilon=0, bounds=(-2, 2))
 
         with self.assertRaises(ValueError):
             clf.fit(X, y)
 
     def test_neg_epsilon(self):
-        X = np.random.random((10, 2))
-        y = np.random.randint(2, size=10)
-        clf = GaussianNB(epsilon=-1, bounds=(0, 1))
+        clf = GaussianNB(epsilon=-1, bounds=(-2, 2))
 
         with self.assertRaises(ValueError):
             clf.fit(X, y)
 
     def test_sample_weight_warning(self):
-        X = np.random.random((10, 2))
-        y = np.random.randint(2, size=10)
-        clf = GaussianNB(epsilon=1, bounds=([0, 0], [1, 1]))
+        clf = GaussianNB(epsilon=1, bounds=(-2, 2))
         w = abs(np.random.randn(10))
 
         with self.assertWarns(DiffprivlibCompatibilityWarning):
             clf.fit(X, y, sample_weight=w)
 
     def test_mis_ordered_bounds(self):
-        X = np.random.random((10, 2))
-        y = np.random.randint(2, size=10)
-
         clf = GaussianNB(epsilon=1, bounds=([0, 1], [1, 0]))
 
         with self.assertRaises(ValueError):
             clf.fit(X, y)
 
     def test_no_bounds(self):
-        X = np.random.random((10, 2))
-        y = np.random.randint(2, size=10)
         clf = GaussianNB()
 
         with self.assertWarns(PrivacyLeakWarning):
@@ -58,8 +51,8 @@ def test_no_bounds(self):
         self.assertIsNotNone(clf)
 
     def test_missing_bounds(self):
-        X = np.random.random((10, 3))
-        y = np.random.randint(2, size=10)
+        rng = check_random_state(0)
+        X = rng.random((10, 6))
 
         clf = GaussianNB(epsilon=1, bounds=([0, 0], [1, 1]))
 
@@ -67,16 +60,17 @@ def test_missing_bounds(self):
             clf.fit(X, y)
 
     def test_noisy_count(self):
-        y = np.random.randint(20, size=10000)
+        rng = check_random_state(0)
+        y = rng.randint(20, size=10000)
         actual_counts = np.array([(y == y_i).sum() for y_i in np.unique(y)])
 
         clf = GaussianNB(epsilon=3)
-        noisy_counts = clf._noisy_class_counts(y)
+        noisy_counts = clf._noisy_class_counts(y, random_state=None)
         self.assertEqual(y.shape[0], noisy_counts.sum())
         self.assertFalse(np.all(noisy_counts == actual_counts))
 
         clf = GaussianNB(epsilon=float("inf"))
-        noisy_counts = clf._noisy_class_counts(y)
+        noisy_counts = clf._noisy_class_counts(y, random_state=None)
         self.assertEqual(y.shape[0], noisy_counts.sum())
         self.assertTrue(np.all(noisy_counts == actual_counts))
 
@@ -85,14 +79,15 @@ def test_different_results(self):
         from sklearn.naive_bayes import GaussianNB as sk_nb
         from sklearn import datasets
 
-        global_seed(12345)
         dataset = datasets.load_iris()
+        rng = check_random_state(0)
 
-        x_train, x_test, y_train, y_test = train_test_split(dataset.data, dataset.target, test_size=.2)
+        x_train, x_test, y_train, y_test = train_test_split(dataset.data, dataset.target, test_size=.2,
+                                                            random_state=rng)
 
         bounds = ([4.3, 2.0, 1.0, 0.1], [7.9, 4.4, 6.9, 2.5])
 
-        clf_dp = GaussianNB(epsilon=1.0, bounds=bounds)
+        clf_dp = GaussianNB(epsilon=1.0, bounds=bounds, random_state=rng)
         clf_non_private = sk_nb()
 
         for clf in [clf_dp, clf_non_private]:
@@ -112,15 +107,16 @@ def test_different_results(self):
 
     @pytest.mark.filterwarnings('ignore: numpy.ufunc size changed')
     def test_with_iris(self):
-        global_seed(12345)
         from sklearn import datasets
         dataset = datasets.load_iris()
+        rng = check_random_state(0)
 
-        x_train, x_test, y_train, y_test = train_test_split(dataset.data, dataset.target, test_size=.2)
+        x_train, x_test, y_train, y_test = train_test_split(dataset.data, dataset.target, test_size=.2,
+                                                            random_state=rng)
 
         bounds = ([4.3, 2.0, 1.0, 0.1], [7.9, 4.4, 6.9, 2.5])
 
-        clf = GaussianNB(epsilon=5.0, bounds=bounds)
+        clf = GaussianNB(epsilon=5.0, bounds=bounds, random_state=rng)
         clf.fit(x_train, y_train)
 
         accuracy = clf.score(x_test, y_test)
@@ -135,81 +131,75 @@ def test_accountant(self):
         from diffprivlib.accountant import BudgetAccountant
         acc = BudgetAccountant()
 
-        x_train = np.random.random((10, 2))
-        y_train = np.random.randint(2, size=10)
-
-        clf = GaussianNB(epsilon=1.0, bounds=(0, 1), accountant=acc)
-        clf.fit(x_train, y_train)
+        clf = GaussianNB(epsilon=1.0, bounds=(-2, 2), accountant=acc)
+        clf.fit(X, y)
         self.assertEqual((1, 0), acc.total())
 
         with BudgetAccountant(1.5, 0) as acc2:
-            clf = GaussianNB(epsilon=1.0, bounds=(0, 1))
-            clf.fit(x_train, y_train)
+            clf = GaussianNB(epsilon=1.0, bounds=(-2, 2))
+            clf.fit(X, y)
             self.assertEqual((1, 0), acc2.total())
 
             with self.assertRaises(BudgetError):
-                clf.fit(x_train, y_train)
+                clf.fit(X, y)
 
     def test_priors(self):
-        X = np.random.random((10, 2))
-        y = np.random.randint(2, size=10)
-
-        clf = GaussianNB(epsilon=1, bounds=([0, 0], [1, 1]), priors=(0.75, 0.25))
+        clf = GaussianNB(epsilon=1, bounds=(-2, 2), priors=(0.75, 0.25))
         self.assertIsNotNone(clf.fit(X, y))
 
-        clf = GaussianNB(epsilon=1, bounds=([0, 0], [1, 1]), priors=(1,))
+        clf = GaussianNB(epsilon=1, bounds=(-2, 2), priors=(1,))
         with self.assertRaises(ValueError):
             clf.fit(X, y)
 
-        clf = GaussianNB(epsilon=1, bounds=([0, 0], [1, 1]), priors=(0.5, 0.7))
+        clf = GaussianNB(epsilon=1, bounds=(-2, 2), priors=(0.5, 0.7))
         with self.assertRaises(ValueError):
             clf.fit(X, y)
 
-        clf = GaussianNB(epsilon=1, bounds=([0, 0], [1, 1]), priors=(-0.5, 1.5))
+        clf = GaussianNB(epsilon=1, bounds=(-2, 2), priors=(-0.5, 1.5))
         with self.assertRaises(ValueError):
             clf.fit(X, y)
 
     def test_bad_refit_shape(self):
-        X = np.random.random((10, 2))
-        y = np.random.randint(2, size=10)
-
-        clf = GaussianNB(epsilon=1, bounds=([0, 0], [1, 1]))
+        clf = GaussianNB(epsilon=1, bounds=(-2, 2))
         clf.fit(X, y)
 
-        X2 = np.random.random((10, 3))
+        X2 = np.random.random((6, 3))
         clf.bounds = ([0, 0, 0], [1, 1, 1])
 
         with self.assertRaises(ValueError):
             clf.partial_fit(X2, y)
 
     def test_bad_refit_classes(self):
-        X = np.random.random((10, 2))
-        y = np.random.randint(2, size=10)
-
-        clf = GaussianNB(epsilon=1, bounds=([0, 0], [1, 1]))
+        clf = GaussianNB(epsilon=1, bounds=(-2, 2))
         clf.fit(X, y)
 
-        X2 = np.random.random((10, 2))
-        y2 = np.random.randint(3, size=10)
+        y2 = np.array([1, 1, 1, 2, 2, 3])
 
         with self.assertRaises(ValueError):
-            clf.partial_fit(X2, y2)
+            clf.partial_fit(X, y2)
 
     def test_update_mean_variance(self):
-        X = np.random.random((10, 2))
-        y = np.random.randint(2, size=10)
-
-        clf = GaussianNB(epsilon=1, bounds=([0, 0], [1, 1]))
-        self.assertIsNotNone(clf._update_mean_variance(0, 0, 0, X, n_noisy=5))
-        self.assertIsNotNone(clf._update_mean_variance(0, 0, 0, X, n_noisy=0))
-        self.assertWarns(PrivacyLeakWarning, clf._update_mean_variance, 0, 0, 0, X)
-        self.assertWarns(DiffprivlibCompatibilityWarning, clf._update_mean_variance, 0, 0, 0, X, n_noisy=1,
+        clf = GaussianNB(epsilon=1, bounds=([-2, -2], [2, 2]))
+        self.assertIsNotNone(clf._update_mean_variance(0, 0, 0, X, None, n_noisy=5))
+        self.assertIsNotNone(clf._update_mean_variance(0, 0, 0, X, None, n_noisy=0))
+        self.assertWarns(PrivacyLeakWarning, clf._update_mean_variance, 0, 0, 0, X, None)
+        self.assertWarns(DiffprivlibCompatibilityWarning, clf._update_mean_variance, 0, 0, 0, X, None, n_noisy=1,
                          sample_weight=1)
 
     def test_sigma(self):
-        X = np.random.random((10, 2))
-        y = np.random.randint(2, size=10)
-
-        clf = GaussianNB(epsilon=1, bounds=([0, 0], [1, 1]))
+        clf = GaussianNB(epsilon=1, bounds=(-2, 2))
         clf.fit(X, y)
         self.assertIsInstance(clf.sigma_, np.ndarray)
+
+    def test_random_state(self):
+        clf0 = GaussianNB(epsilon=1, bounds=(-2, 2), random_state=0)
+        clf1 = GaussianNB(epsilon=1, bounds=(-2, 2), random_state=1)
+        clf0.fit(X, y)
+        clf1.fit(X, y)
+        self.assertFalse(np.all(clf0.theta_ == clf1.theta_))
+        self.assertFalse(np.all(clf0.var_ == clf1.var_))
+
+        clf1 = GaussianNB(epsilon=1, bounds=(-2, 2), random_state=0)
+        clf1.fit(X, y)
+        self.assertTrue(np.all(clf0.theta_ == clf1.theta_))
+        self.assertTrue(np.all(clf0.var_ == clf1.var_))
diff --git a/tests/models/test_KMeans.py b/tests/models/test_KMeans.py
index a3de3393..4dd4788d 100644
--- a/tests/models/test_KMeans.py
+++ b/tests/models/test_KMeans.py
@@ -2,7 +2,7 @@
 from unittest import TestCase
 
 from diffprivlib.models.k_means import KMeans
-from diffprivlib.utils import global_seed, PrivacyLeakWarning, DiffprivlibCompatibilityWarning, BudgetError
+from diffprivlib.utils import PrivacyLeakWarning, DiffprivlibCompatibilityWarning, BudgetError
 
 
 class TestKMeans(TestCase):
@@ -10,12 +10,11 @@ def test_not_none(self):
         self.assertIsNotNone(KMeans)
 
     def test_simple(self):
-        global_seed(3141592653)
-        clf = KMeans(3, epsilon=10, bounds=(0, 1))
+        clf = KMeans(3, epsilon=5, bounds=(0, 1), random_state=0)
 
-        X = np.zeros(6000) + 0.1
-        X[:4000] = 0.5
-        X[:2000] = 0.9
+        X = np.zeros(600) + 0.1
+        X[:400] = 0.5
+        X[:200] = 0.9
         X = X.reshape(-1, 1)
 
         clf.fit(X)
@@ -46,7 +45,7 @@ def test_no_bounds(self):
             clf.fit(X)
 
     def test_predict(self):
-        clf = KMeans(3, epsilon=1, bounds=(0, 1))
+        clf = KMeans(3, epsilon=1, bounds=(0, 1), random_state=0)
 
         X = np.array([0.1, 0.1, 0.1, 0.1, 0.5, 0.5, 0.5, 0.5, 0.9, 0.9, 0.9]).reshape(-1, 1)
         clf.fit(X)
@@ -73,22 +72,22 @@ def test_too_many_clusters(self):
         self.assertRaises(ValueError, clf.fit, X)
 
     def test_inf_epsilon(self):
-        global_seed(3141592653)
-        clf = KMeans(3, epsilon=float("inf"), bounds=(0, 1))
+        clf = KMeans(3, epsilon=float("inf"), bounds=(0, 1), random_state=0)
 
-        X = np.array([0.1, 0.1, 0.1, 0.1, 0.5, 0.5, 0.5, 0.5, 0.9, 0.9, 0.9]).reshape(-1, 1)
+        X = np.array([0.1, 0.1, 0.1, 0.1, 0.5, 0.5, 0.5, 0.5, 0.9, 0.9, 0.9] * 3).reshape(-1, 1)
         clf.fit(X)
-        centers = clf.cluster_centers_
+        centers = clf.cluster_centers_.flatten()
 
-        self.assertIn(0.1, centers)
-        self.assertIn(0.5, centers)
-        self.assertIn(0.9, centers)
+        self.assertTrue(np.isclose(0.1, centers).any())
+        self.assertTrue(np.isclose(0.5, centers).any())
+        self.assertTrue(np.isclose(0.9, centers).any())
 
     def test_many_features(self):
-        X = np.random.random(size=(500, 3))
+        rng = np.random.RandomState(0)
+        X = rng.random(size=(500, 3))
         bounds = (0, 1)
 
-        clf = KMeans(4, bounds=bounds)
+        clf = KMeans(4, bounds=bounds, random_state=rng)
 
         clf.fit(X)
         centers = clf.cluster_centers_
@@ -98,6 +97,21 @@ def test_many_features(self):
 
         self.assertTrue(np.all(clf.transform(X).argmin(axis=1) == clf.predict(X)))
 
+    def test_random_state(self):
+        rng = np.random.RandomState(0)
+        X = rng.random(size=(500, 3))
+        bounds = (0, 1)
+
+        clf0 = KMeans(4, bounds=bounds, random_state=0)
+        clf1 = KMeans(4, bounds=bounds, random_state=1)
+        clf0.fit(X)
+        clf1.fit(X)
+        self.assertFalse(np.any(clf0.cluster_centers_ == clf1.cluster_centers_))
+
+        clf1 = KMeans(4, bounds=bounds, random_state=0)
+        clf1.fit(X)
+        self.assertTrue(np.all(clf0.cluster_centers_ == clf1.cluster_centers_))
+
     def test_accountant(self):
         from diffprivlib.accountant import BudgetAccountant
 
diff --git a/tests/models/test_LinearRegression.py b/tests/models/test_LinearRegression.py
index d57a85fc..15f286d8 100644
--- a/tests/models/test_LinearRegression.py
+++ b/tests/models/test_LinearRegression.py
@@ -4,14 +4,10 @@
 import pytest
 
 from diffprivlib.models.linear_regression import LinearRegression
-from diffprivlib.models.logistic_regression import _check_solver, _check_multi_class
-from diffprivlib.utils import global_seed, PrivacyLeakWarning, DiffprivlibCompatibilityWarning, BudgetError
+from diffprivlib.utils import PrivacyLeakWarning, DiffprivlibCompatibilityWarning, BudgetError
 
 
 class TestLinearRegression(TestCase):
-    def setup_method(self, method):
-        global_seed(3141592653)
-
     def test_not_none(self):
         self.assertIsNotNone(LinearRegression)
 
@@ -68,15 +64,17 @@ def test_different_results(self):
         from sklearn import linear_model
         from sklearn.model_selection import train_test_split
 
+        rng = np.random.RandomState(0)
         dataset = datasets.load_iris()
-        X_train, X_test, y_train, y_test = train_test_split(dataset.data, dataset.target, test_size=0.2)
+        X_train, X_test, y_train, y_test = train_test_split(dataset.data, dataset.target, test_size=0.2,
+                                                            random_state=rng)
 
-        clf = LinearRegression(bounds_X=([4.3, 2.0, 1.1, 0.1], [7.9, 4.4, 6.9, 2.5]), bounds_y=(0, 2))
+        clf = LinearRegression(bounds_X=([4.3, 2.0, 1.1, 0.1], [7.9, 4.4, 6.9, 2.5]), bounds_y=(0, 2), random_state=rng)
         clf.fit(X_train, y_train)
 
         predict1 = clf.predict(X_test)
 
-        clf = LinearRegression(bounds_X=([4.3, 2.0, 1.1, 0.1], [7.9, 4.4, 6.9, 2.5]), bounds_y=(0, 2))
+        clf = LinearRegression(bounds_X=([4.3, 2.0, 1.1, 0.1], [7.9, 4.4, 6.9, 2.5]), bounds_y=(0, 2), random_state=rng)
         clf.fit(X_train, y_train)
 
         predict2 = clf.predict(X_test)
@@ -117,7 +115,7 @@ def test_simple(self):
         y = X.copy()
         X = X[:, np.newaxis]
 
-        clf = LinearRegression(epsilon=2, fit_intercept=False, bounds_X=(-1, 1), bounds_y=(-1, 1))
+        clf = LinearRegression(epsilon=2, fit_intercept=False, bounds_X=(-1, 1), bounds_y=(-1, 1), random_state=0)
         clf.fit(X, y)
         print(clf.predict(np.array([0.5]).reshape(-1, 1)))
 
@@ -163,23 +161,20 @@ def test_multiple_targets(self):
         self.assertEqual(clf_dp2.coef_.shape, clf_sk.coef_.shape)
         self.assertEqual(clf_dp2.predict(x0).shape, clf_sk.predict(x0).shape)
 
-    def test_check_solver(self):
-        with self.assertWarns(DiffprivlibCompatibilityWarning):
-            _check_solver("wrong_solver", "l2", False)
-
-        with self.assertRaises(ValueError):
-            _check_solver("lbfgs", "l1", False)
-
-        with self.assertRaises(ValueError):
-            _check_solver("lbfgs", "l2", True)
-
-        self.assertIsNotNone(_check_solver("lbfgs", "l2", False))
+    def test_random_state(self):
+        X = np.linspace(-1, 1, 1000)
+        y = X.copy()
+        X = X[:, np.newaxis]
 
-    def test_check_multi_class(self):
-        with self.assertWarns(DiffprivlibCompatibilityWarning):
-            _check_multi_class("multinomial", "lbfgs", 2)
+        clf0 = LinearRegression(epsilon=2, fit_intercept=False, bounds_X=(-1, 1), bounds_y=(-1, 1), random_state=0)
+        clf1 = LinearRegression(epsilon=2, fit_intercept=False, bounds_X=(-1, 1), bounds_y=(-1, 1), random_state=1)
+        clf0.fit(X, y)
+        clf1.fit(X, y)
+        self.assertFalse(np.any(clf0.coef_ == clf1.coef_))
 
-        self.assertIsNotNone(_check_multi_class("ovr", "lbfgs", 2))
+        clf1 = LinearRegression(epsilon=2, fit_intercept=False, bounds_X=(-1, 1), bounds_y=(-1, 1), random_state=0)
+        clf1.fit(X, y)
+        self.assertTrue(np.all(clf0.coef_ == clf1.coef_))
 
     def test_accountant(self):
         from diffprivlib.accountant import BudgetAccountant
diff --git a/tests/models/test_LogisticRegression.py b/tests/models/test_LogisticRegression.py
index d32d9d4f..06d4ba15 100644
--- a/tests/models/test_LogisticRegression.py
+++ b/tests/models/test_LogisticRegression.py
@@ -2,13 +2,10 @@
 from unittest import TestCase
 
 from diffprivlib.models.logistic_regression import LogisticRegression
-from diffprivlib.utils import global_seed, PrivacyLeakWarning, DiffprivlibCompatibilityWarning, BudgetError
+from diffprivlib.utils import PrivacyLeakWarning, DiffprivlibCompatibilityWarning, BudgetError
 
 
 class TestLogisticRegression(TestCase):
-    def setup_method(self, method):
-        global_seed(3141592653)
-
     def test_not_none(self):
         self.assertIsNotNone(LogisticRegression)
 
@@ -131,20 +128,23 @@ def test_different_results(self):
         from sklearn import linear_model
         from sklearn.model_selection import train_test_split
 
+        rng = np.random.RandomState(0)
+
         dataset = datasets.load_iris()
-        X_train, X_test, y_train, y_test = train_test_split(dataset.data, dataset.target, test_size=0.2)
+        X_train, X_test, y_train, y_test = train_test_split(dataset.data, dataset.target, test_size=0.2,
+                                                            random_state=rng)
 
-        clf = LogisticRegression(data_norm=12)
+        clf = LogisticRegression(data_norm=12, random_state=rng)
         clf.fit(X_train, y_train)
 
         predict1 = clf.predict(X_test)
 
-        clf = LogisticRegression(data_norm=12)
+        clf = LogisticRegression(data_norm=12, random_state=rng)
         clf.fit(X_train, y_train)
 
         predict2 = clf.predict(X_test)
 
-        clf = linear_model.LogisticRegression(solver="lbfgs", multi_class="ovr")
+        clf = linear_model.LogisticRegression(solver="lbfgs", multi_class="ovr", random_state=rng)
         clf.fit(X_train, y_train)
 
         predict3 = clf.predict(X_test)
@@ -174,19 +174,40 @@ def test_same_results(self):
     def test_simple(self):
         X = np.array(
             [0.50, 0.75, 1.00, 1.25, 1.50, 1.75, 1.75, 2.00, 2.25, 2.50, 2.75, 3.00, 3.25, 3.50, 4.00, 4.25, 4.50, 4.75,
-             5.00, 5.50] * 5)
-        y = np.array([0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1] * 5)
+             5.00, 5.50] * 3)
+        y = np.array([0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1] * 3)
         X = X[:, np.newaxis]
         X -= 3.0
         X /= 2.5
 
-        clf = LogisticRegression(epsilon=2, data_norm=1.0)
+        clf = LogisticRegression(epsilon=2, data_norm=1.0, random_state=0)
         clf.fit(X, y)
 
         self.assertIsNotNone(clf)
         self.assertFalse(clf.predict(np.array([(0.5 - 3) / 2.5]).reshape(-1, 1)))
         self.assertTrue(clf.predict(np.array([(5.5 - 3) / 2.5]).reshape(-1, 1)))
 
+    def test_random_state(self):
+        X = np.array(
+            [0.50, 0.75, 1.00, 1.25, 1.50, 1.75, 1.75, 2.00, 2.25, 2.50, 2.75, 3.00, 3.25, 3.50, 4.00, 4.25, 4.50, 4.75,
+             5.00, 5.50] * 3)
+        y = np.array([0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1] * 3)
+        X = X[:, np.newaxis]
+        X -= 3.0
+        X /= 2.5
+
+        clf0 = LogisticRegression(epsilon=2, data_norm=1.0, random_state=0)
+        clf1 = LogisticRegression(epsilon=2, data_norm=1.0, random_state=1)
+        clf0.fit(X, y)
+        clf1.fit(X, y)
+        self.assertFalse(np.any(clf0.coef_ == clf1.coef_))
+        self.assertFalse(np.any(clf0.intercept_ == clf1.intercept_))
+
+        clf1 = LogisticRegression(epsilon=2, data_norm=1.0, random_state=0)
+        clf1.fit(X, y)
+        self.assertTrue(np.all(clf0.coef_ == clf1.coef_))
+        self.assertTrue(np.all(clf0.intercept_ == clf1.intercept_))
+
     def test_accountant(self):
         from diffprivlib.accountant import BudgetAccountant
         acc = BudgetAccountant()
diff --git a/tests/models/test_PCA.py b/tests/models/test_PCA.py
index c65da3c9..70cb8e2c 100644
--- a/tests/models/test_PCA.py
+++ b/tests/models/test_PCA.py
@@ -9,7 +9,7 @@
     import sklearn.decomposition.pca as sk_pca
 
 from diffprivlib.models.pca import PCA
-from diffprivlib.utils import PrivacyLeakWarning, DiffprivlibCompatibilityWarning, BudgetError
+from diffprivlib.utils import PrivacyLeakWarning, DiffprivlibCompatibilityWarning, BudgetError, check_random_state
 
 
 class TestPCA(TestCase):
@@ -73,15 +73,16 @@ def test_solver_warning(self):
             PCA(svd_solver='full')
 
     def test_inf_epsilon(self):
-        X = np.random.randn(250, 21)
+        rng = np.random.RandomState(0)
+        X = rng.randn(250, 21)
         X -= np.mean(X, axis=0)
         X /= np.linalg.norm(X, axis=1).max()
 
         for components in range(1, 10):
-            clf = PCA(components, epsilon=float("inf"), centered=True, data_norm=1)
+            clf = PCA(components, epsilon=float("inf"), centered=True, data_norm=1, random_state=rng)
             clf.fit(X)
 
-            sk_clf = sk_pca.PCA(components, svd_solver='full')
+            sk_clf = sk_pca.PCA(components, svd_solver='full', random_state=rng)
             sk_clf.fit(X)
 
             self.assertAlmostEqual(clf.score(X), sk_clf.score(X))
@@ -90,15 +91,16 @@ def test_inf_epsilon(self):
             self.assertTrue(np.allclose(np.abs((clf.components_ / sk_clf.components_).sum(axis=1)), clf.n_features_))
 
     def test_big_epsilon(self):
-        X = np.random.randn(5000, 10)
+        rng = check_random_state(2)
+        X = rng.randn(2000, 10)
         X -= np.mean(X, axis=0)
         X /= np.linalg.norm(X, axis=1).max()
 
         for components in range(1, 3):
-            clf = PCA(components * 3, epsilon=10, centered=True, data_norm=1)
+            clf = PCA(components * 3, epsilon=10, centered=True, data_norm=1, random_state=rng)
             clf.fit(X)
 
-            sk_clf = sk_pca.PCA(components * 3, svd_solver='full')
+            sk_clf = sk_pca.PCA(components * 3, svd_solver='full', random_state=rng)
             sk_clf.fit(X)
 
             self.assertAlmostEqual(clf.score(X) / sk_clf.score(X), 1, places=2)
@@ -124,24 +126,27 @@ def test_different_results(self):
         from sklearn import datasets
         from sklearn.model_selection import train_test_split
 
+        rng = np.random.RandomState(0)
+
         dataset = datasets.load_iris()
-        X_train, X_test, y_train, y_test = train_test_split(dataset.data, dataset.target, test_size=0.2)
+        X_train, X_test, y_train, y_test = train_test_split(dataset.data, dataset.target, test_size=0.2,
+                                                            random_state=rng)
 
         mean = np.mean(X_train, axis=0)
         X_train -= mean
         X_test -= mean
 
-        clf = PCA(data_norm=12, centered=True)
+        clf = PCA(data_norm=12, centered=True, random_state=rng)
         clf.fit(X_train, y_train)
 
         transform1 = clf.transform(X_test)
 
-        clf = PCA(data_norm=12, centered=True)
+        clf = PCA(data_norm=12, centered=True, random_state=rng)
         clf.fit(X_train, y_train)
 
         transform2 = clf.transform(X_test)
 
-        clf = sk_pca.PCA(svd_solver='full')
+        clf = sk_pca.PCA(svd_solver='full', random_state=rng)
         clf.fit(X_train, y_train)
 
         transform3 = clf.transform(X_test)
@@ -149,6 +154,24 @@ def test_different_results(self):
         self.assertFalse(np.all(transform1 == transform2))
         self.assertFalse(np.all(transform3 == transform1) and np.all(transform3 == transform2))
 
+    def test_random_state(self):
+        X = np.array(
+            [0.50, 0.75, 1.00, 1.25, 1.50, 1.75, 1.75, 2.00, 2.25, 2.50, 2.75, 3.00, 3.25, 3.50, 4.00, 4.25, 4.50, 4.75,
+             5.00, 5.50])
+        X = X[:, np.newaxis]
+
+        clf0 = PCA(data_norm=1.0, centered=True, random_state=0)
+        clf1 = PCA(data_norm=1.0, centered=True, random_state=1)
+        clf0.fit(X)
+        clf1.fit(X)
+        self.assertFalse(np.any(clf1.singular_values_ == clf0.singular_values_))
+        self.assertFalse(np.any(clf1.explained_variance_ == clf0.explained_variance_))
+
+        clf1 = PCA(data_norm=1.0, centered=True, random_state=0)
+        clf1.fit(X)
+        self.assertTrue(np.all(clf1.singular_values_ == clf0.singular_values_))
+        self.assertTrue(np.all(clf1.explained_variance_ == clf0.explained_variance_))
+
     def test_accountant(self):
         from diffprivlib.accountant import BudgetAccountant
         acc = BudgetAccountant()
diff --git a/tests/models/test_RandomForestClassifier.py b/tests/models/test_RandomForestClassifier.py
index faaad36e..a8331d4e 100644
--- a/tests/models/test_RandomForestClassifier.py
+++ b/tests/models/test_RandomForestClassifier.py
@@ -6,7 +6,7 @@
 from sklearn.exceptions import NotFittedError
 
 from diffprivlib.models.forest import RandomForestClassifier, DecisionTreeClassifier
-from diffprivlib.utils import PrivacyLeakWarning, global_seed, BudgetError, DiffprivlibCompatibilityWarning
+from diffprivlib.utils import PrivacyLeakWarning, BudgetError, DiffprivlibCompatibilityWarning
 
 X = [[-2, -1], [-1, -1], [-1, -2], [1, 1], [1, 2], [2, 1]] * 3
 y = [-1, -1, -1, 1, 1, 1] * 3
@@ -15,9 +15,6 @@
 
 
 class TestRandomForestClassifier(TestCase):
-    def setUp(self):
-        global_seed(2718281828)
-
     def test_not_none(self):
         self.assertIsNotNone(RandomForestClassifier)
 
@@ -25,7 +22,7 @@ def test_bad_params(self):
         X = [[1]]
         y = [0]
 
-        with self.assertRaises((ValueError, TypeError)):  # Should be TypeError
+        with self.assertRaises((ValueError, TypeError)):  # Should be TypeError?
             RandomForestClassifier(n_estimators="10", bounds=(0, 1), classes=[0, 1]).fit(X, y)
 
         with self.assertWarns(DiffprivlibCompatibilityWarning):
@@ -55,22 +52,24 @@ def test_multi_output(self):
     def test_simple(self):
         X = np.array([[12, 3, 14], [12, 3, 4], [12, 3, 4], [2, 13, 4], [2, 13, 14], [2, 3, 14], [3, 5, 15]] * 3)
         y = np.array([1, 1, 1, 0, 0, 0, 1] * 3)
-        model = RandomForestClassifier(epsilon=5, n_estimators=5, max_depth=3)
+        model = RandomForestClassifier(epsilon=5, n_estimators=5, max_depth=3, random_state=1)
 
         with self.assertRaises(NotFittedError):
             check_is_fitted(model)
-        # when `feature_domains` is not provided, we should get a privacy leakage warning
 
+        # when `bounds` is not provided, we should get a privacy leakage warning
         with self.assertWarns(PrivacyLeakWarning):
             model.fit(X, y)
         check_is_fitted(model)
 
+        # Test parameters and attributes
         self.assertEqual(model.n_features_in_, 3)
         self.assertEqual(model.n_classes_, 2)
         self.assertEqual(set(model.classes_), set([0, 1]))
         self.assertEqual(model.max_depth, 3)
         self.assertTrue(model.estimators_)
         self.assertEqual(len(model.estimators_), 5)
+
         self.assertTrue(model.predict(np.array([[12, 3, 14]])))
         self.assertTrue(np.all(model.bounds[0] == [2, 3, 4]))
         self.assertTrue(np.all(model.bounds[1] == [12, 13, 15]))
@@ -79,13 +78,15 @@ def test_with_bounds(self):
         X = np.array([[12, 3, 14], [12, 3, 4], [12, 3, 4], [2, 13, 4], [2, 13, 14], [2, 3, 14], [3, 5, 15]] * 3)
         y = np.array([1, 1, 1, 0, 0, 0, 1] * 3)
         model = RandomForestClassifier(epsilon=5, n_estimators=5, bounds=([2, 3, 4], [12, 13, 15]), classes=[0, 1],
-                                       max_depth=3)
+                                       max_depth=3, random_state=1)
         with self.assertRaises(NotFittedError):
             check_is_fitted(model)
 
         model.fit(X, y)
 
         check_is_fitted(model)
+
+        # Test parameters and attributes
         self.assertEqual(model.n_features_in_, 3)
         self.assertEqual(model.n_classes_, 2)
         self.assertEqual(set(model.classes_), set([0, 1]))
@@ -113,16 +114,41 @@ def test_non_binary_classes(self):
     def test_with_not_enough_bounds(self):
         X = np.array([[12, 3, 14], [12, 3, 4], [12, 3, 4], [2, 13, 4], [2, 13, 14], [2, 3, 14], [3, 5, 15]])
         y = np.array([1, 1, 1, 0, 0, 0, 1])
-        model = RandomForestClassifier(epsilon=2, n_estimators=5, bounds=([2, 3], [12, 13]))
+        model = RandomForestClassifier(epsilon=2, n_estimators=5, bounds=([2, 3], [12, 13]), classes=[0, 1])
+
         with self.assertRaises(NotFittedError):
             check_is_fitted(model)
+
         with self.assertRaises(ValueError):
             model.fit(X, y)
 
+    def test_random_state(self):
+        model0 = RandomForestClassifier(20, epsilon=0.1, bounds=(-2, 2), classes=[-1, 1], max_depth=3,
+                                        random_state=0)
+        model1 = RandomForestClassifier(20, epsilon=0.1, bounds=(-2, 2), classes=[-1, 1], max_depth=3,
+                                        random_state=2)
+        model0.fit(X, y)
+        model1.fit(X, y)
+
+        self.assertFalse(np.isclose(model0.predict_proba(np.array([[-2, -1]])),
+                                    model1.predict_proba(np.array([[-2, -1]]))).any(),
+                         (model0.predict_proba(np.array([[-2, -1]])), model1.predict_proba(np.array([[-2, -1]]))))
+
+        model1 = RandomForestClassifier(20, epsilon=0.1, bounds=(-2, 2), classes=[-1, 1], max_depth=3,
+                                        random_state=0)
+        model1.fit(X, y)
+        for i in range(20):
+            self.assertIsNotNone(model0.estimators_[i].random_state)
+            self.assertEqual(model0.estimators_[i].random_state, model1.estimators_[i].random_state)
+
+        self.assertTrue(np.isclose(model0.predict_proba(np.array([[-2, -1]])),
+                                   model1.predict_proba(np.array([[-2, -1]]))).any(),
+                        (model0.predict_proba(np.array([[-2, -1]])), model1.predict_proba(np.array([[-2, -1]]))))
+
     def test_sklearn_methods(self):
         depth = 3
         n_estimators = 4
-        clf = RandomForestClassifier(n_estimators, epsilon=5, bounds=(-2, 2), classes=np.array((-1, 1)),
+        clf = RandomForestClassifier(n_estimators, epsilon=5, bounds=(-2, 2), classes=[-1, 1],
                                      max_depth=depth)
         clf.fit(X, y)
 
@@ -155,8 +181,8 @@ def test_sklearn_methods(self):
     def test_warm_start(self):
         depth = 3
         n_estimators = 4
-        clf = RandomForestClassifier(n_estimators, epsilon=1, bounds=(-2, 2), classes=np.array((-1, 1)),
-                                     max_depth=depth, warm_start=True)
+        clf = RandomForestClassifier(n_estimators, epsilon=1, bounds=(-2, 2), classes=[-1, 1],
+                                     max_depth=depth, warm_start=True, random_state=0)
         clf.fit(X, y)
         first_estimator = clf[0].__dict__.copy()
 
@@ -174,10 +200,18 @@ def test_warm_start(self):
         with self.assertRaises(ValueError):
             clf.fit(X, y)
 
+        # Estimators should still have the same seed
+        clf2 = RandomForestClassifier(n_estimators * 2, epsilon=1, bounds=(-2, 2), classes=[-1, 1],
+                                      max_depth=depth, warm_start=True, random_state=0)
+        clf2.fit(X, y)
+
+        for i in range(n_estimators * 2):
+            self.assertEqual(clf.estimators_[i].random_state, clf2.estimators_[i].random_state)
+
     def test_parallel(self):
         depth = 3
         n_estimators = 4
-        clf = RandomForestClassifier(n_estimators, epsilon=5, bounds=(-2, 2), classes=np.array((-1, 1)),
+        clf = RandomForestClassifier(n_estimators, epsilon=5, bounds=(-2, 2), classes=[-1, 1],
                                      max_depth=depth, n_jobs=-1)
         clf.fit(X, y)
 
@@ -186,16 +220,27 @@ def test_parallel(self):
     def test_shuffle(self):
         depth = 3
         n_estimators = 4
-        clf = RandomForestClassifier(n_estimators, epsilon=5, bounds=(-2, 2), classes=np.array((-1, 1)),
-                                     max_depth=depth, shuffle=True)
+        clf = RandomForestClassifier(n_estimators, epsilon=5, bounds=(-2, 2), classes=[-1, 1],
+                                     max_depth=depth, random_state=0, shuffle=False)
         clf.fit(X, y)
 
-        self.assertIsNone(check_is_fitted(clf))
+        clf_shuffle = RandomForestClassifier(n_estimators, epsilon=5, bounds=(-2, 2), classes=[-1, 1],
+                                             max_depth=depth, random_state=0, shuffle=True)
+        clf_shuffle.fit(X, y)
+        self.assertIsNone(check_is_fitted(clf_shuffle))
+
+        # Estimators should have the same seed
+        for i in range(n_estimators):
+            self.assertEqual(clf.estimators_[i].random_state, clf_shuffle.estimators_[i].random_state)
+
+        # But predictions should be different, since trained on different data
+        self.assertFalse(np.isclose(clf.predict_proba(T), clf_shuffle.predict_proba(T)).all(),
+                         (clf.predict_proba(T), clf_shuffle.predict_proba(T)))
 
     def test_more_estimators_than_data(self):
         depth = 3
         n_estimators = 20
-        clf = RandomForestClassifier(n_estimators, epsilon=5, bounds=(-2, 2), classes=np.array((-1, 1)),
+        clf = RandomForestClassifier(n_estimators, epsilon=5, bounds=(-2, 2), classes=[-1, 1],
                                      max_depth=depth)
         clf.fit(X, y)
 
diff --git a/tests/models/test_StandardScaler.py b/tests/models/test_StandardScaler.py
index 72c473ec..ea371d66 100644
--- a/tests/models/test_StandardScaler.py
+++ b/tests/models/test_StandardScaler.py
@@ -4,7 +4,7 @@
 import sklearn.preprocessing as sk_pp
 
 from diffprivlib.models.standard_scaler import StandardScaler
-from diffprivlib.utils import PrivacyLeakWarning, DiffprivlibCompatibilityWarning, global_seed, BudgetError
+from diffprivlib.utils import PrivacyLeakWarning, DiffprivlibCompatibilityWarning, BudgetError
 
 
 class TestStandardScaler(TestCase):
@@ -88,11 +88,10 @@ def test_functionality(self):
         self.assertIsNotNone(ss.fit_transform(X))
 
     def test_similar_results(self):
-        global_seed(314159)
+        rng = np.random.RandomState(0)
+        X = rng.rand(100000, 5)
 
-        X = np.random.rand(100000, 5)
-
-        dp_ss = StandardScaler(bounds=(0, 1), epsilon=float("inf"))
+        dp_ss = StandardScaler(bounds=(0, 1), epsilon=float("inf"), random_state=rng)
         dp_ss.fit(X)
 
         sk_ss = sk_pp.StandardScaler()
@@ -104,6 +103,22 @@ def test_similar_results(self):
                         (dp_ss.var_, sk_ss.var_))
         self.assertTrue(np.all(dp_ss.n_samples_seen_ == sk_ss.n_samples_seen_))
 
+    def test_random_state(self):
+        rng = np.random.RandomState(0)
+        X = rng.rand(100000, 5)
+
+        ss0 = StandardScaler(bounds=(0, 1), epsilon=1, random_state=0)
+        ss1 = StandardScaler(bounds=(0, 1), epsilon=1, random_state=1)
+        ss0.fit(X)
+        ss1.fit(X)
+        self.assertFalse(np.any(ss0.mean_ == ss1.mean_))
+        self.assertFalse(np.any(ss0.var_ == ss1.var_))
+
+        ss1 = StandardScaler(bounds=(0, 1), epsilon=1, random_state=0)
+        ss1.fit(X)
+        self.assertTrue(np.all(ss0.mean_ == ss1.mean_), (ss0.mean_, ss1.mean_))
+        self.assertTrue(np.all(ss0.var_ == ss1.var_))
+
     def test_accountant(self):
         from diffprivlib.accountant import BudgetAccountant
         acc = BudgetAccountant()
diff --git a/tests/models/test_logistic_regression_path.py b/tests/models/test_logistic_regression_path.py
index 8b723080..b7b01457 100644
--- a/tests/models/test_logistic_regression_path.py
+++ b/tests/models/test_logistic_regression_path.py
@@ -10,7 +10,7 @@ def setup_class(self):
         X = np.array(
             [0.50, 0.75, 1.00, 1.25, 1.50, 1.75, 1.75, 2.00, 2.25, 2.50, 2.75, 3.00, 3.25, 3.50, 4.00, 4.25, 4.50, 4.75,
              5.00, 5.50])
-        y = np.array([0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1])
+        y = np.array([0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1], dtype=X.dtype)
         X = X[:, np.newaxis]
 
         self.X = X
@@ -35,7 +35,7 @@ def test_with_dataset(self):
     def test_pos_class(self):
         X = np.array(
             [0.50, 0.75, 1.00])
-        y = np.array([0, 1, 2])
+        y = np.array([0, 1, 2], dtype=X.dtype)
         X = X[:, np.newaxis]
 
         with self.assertRaises(ValueError):
@@ -53,3 +53,11 @@ def test_Cs(self):
         self.assertEqual(len(coefs), 3)
         self.assertEqual(len(Cs), 3)
         self.assertEqual(len(n_iter), 3)
+
+    def test_random_state(self):
+        out0, _, _ = _logistic_regression_path(self.X, self.y, epsilon=1, data_norm=1, Cs=[1e5], random_state=0)
+        out1, _, _ = _logistic_regression_path(self.X, self.y, epsilon=1, data_norm=1, Cs=[1e5], random_state=1)
+        self.assertFalse(np.any(out0 == out1))
+
+        out1, _, _ = _logistic_regression_path(self.X, self.y, epsilon=1, data_norm=1, Cs=[1e5], random_state=0)
+        self.assertTrue(np.all(out0 == out1))
diff --git a/tests/test_check_random_state.py b/tests/test_check_random_state.py
new file mode 100644
index 00000000..2cb51fcb
--- /dev/null
+++ b/tests/test_check_random_state.py
@@ -0,0 +1,29 @@
+import secrets
+from unittest import TestCase
+
+import numpy as np
+
+from diffprivlib.utils import check_random_state
+
+
+class TestCheckRandomState(TestCase):
+    def test_secure(self):
+        rng = secrets.SystemRandom()
+        self.assertIs(rng, check_random_state(rng, True))
+        self.assertIsNot(rng, check_random_state(None, True))
+        self.assertIsInstance(check_random_state(None, True), secrets.SystemRandom)
+
+        self.assertRaises(ValueError, check_random_state, rng, False)
+
+    def test_not_secure(self):
+        rng = np.random.mtrand._rand
+        self.assertIs(check_random_state(None), rng)
+        self.assertIsInstance(check_random_state(1), np.random.RandomState)
+
+        self.assertIsNot(check_random_state(1), rng)
+
+        rng = np.random.RandomState(1)
+        self.assertIs(check_random_state(rng), rng)
+        self.assertIsNot(check_random_state(1), rng)
+        self.assertEqual(rng.random(), check_random_state(1).random())
+        self.assertNotEqual(rng.random(), check_random_state(1).random())
diff --git a/tests/tools/test_histogram.py b/tests/tools/test_histogram.py
index acb6cea9..d015ccab 100644
--- a/tests/tools/test_histogram.py
+++ b/tests/tools/test_histogram.py
@@ -3,7 +3,7 @@
 
 from diffprivlib.accountant import BudgetAccountant
 from diffprivlib.tools.histograms import histogram
-from diffprivlib.utils import global_seed, BudgetError
+from diffprivlib.utils import BudgetError
 
 
 class TestHistogram(TestCase):
@@ -26,7 +26,6 @@ def test_same_edges(self):
         self.assertTrue((edges == dp_edges).all())
 
     def test_different_result(self):
-        global_seed(3141592653)
         a = np.array([1, 2, 3, 4, 5])
         hist, _ = np.histogram(a, bins=3, range=(0, 10))
         dp_hist, _ = histogram(a, epsilon=0.01, bins=3, range=(0, 10))
@@ -36,7 +35,6 @@ def test_different_result(self):
         self.assertTrue((hist != dp_hist).any())
 
     def test_density(self):
-        global_seed(3141592653)
         a = np.array([1, 2, 3, 4, 5])
         dp_hist, _ = histogram(a, epsilon=10, bins=3, range=(0, 10), density=True)
 
diff --git a/tests/tools/test_histogram2d.py b/tests/tools/test_histogram2d.py
index 5be40c75..bd4efcce 100644
--- a/tests/tools/test_histogram2d.py
+++ b/tests/tools/test_histogram2d.py
@@ -3,7 +3,7 @@
 
 from diffprivlib.accountant import BudgetAccountant
 from diffprivlib.tools.histograms import histogram2d
-from diffprivlib.utils import global_seed, PrivacyLeakWarning, BudgetError
+from diffprivlib.utils import PrivacyLeakWarning, BudgetError
 
 
 class TestHistogram2d(TestCase):
@@ -50,7 +50,6 @@ def test_same_edges(self):
         self.assertTrue((edges_y == dp_edges_y).all())
 
     def test_different_result(self):
-        global_seed(3141592653)
         x = np.array([1, 2, 3, 4, 5])
         y = np.array([5, 7, 1, 5, 9])
         hist, _, _ = np.histogram2d(x, y, bins=3, range=[(0, 10), (0, 10)])
@@ -61,7 +60,6 @@ def test_different_result(self):
         self.assertTrue((hist != dp_hist).any())
 
     def test_density(self):
-        global_seed(3141592653)
         x = np.array([1, 2, 3, 4, 5])
         y = np.array([5, 7, 1, 5, 9])
         dp_hist, _, _ = histogram2d(x, y, epsilon=1, bins=3, range=[(0, 10), (0, 10)], density=True)
diff --git a/tests/tools/test_histogramdd.py b/tests/tools/test_histogramdd.py
index f737a556..bb57a756 100644
--- a/tests/tools/test_histogramdd.py
+++ b/tests/tools/test_histogramdd.py
@@ -3,7 +3,7 @@
 
 from diffprivlib.accountant import BudgetAccountant
 from diffprivlib.tools.histograms import histogramdd
-from diffprivlib.utils import global_seed, PrivacyLeakWarning, BudgetError
+from diffprivlib.utils import PrivacyLeakWarning, BudgetError
 
 
 class TestHistogramdd(TestCase):
@@ -33,7 +33,6 @@ def test_same_edges(self):
             self.assertTrue((edges[i] == dp_edges[i]).all())
 
     def test_different_result(self):
-        global_seed(3141592653)
         a = np.array([1, 2, 3, 4, 5])
         hist, _ = np.histogramdd(a, bins=3, range=[(0, 10)])
         dp_hist, _ = histogramdd(a, epsilon=0.1, bins=3, range=[(0, 10)])
@@ -43,7 +42,6 @@ def test_different_result(self):
         self.assertTrue((hist != dp_hist).any())
 
     def test_density_1d(self):
-        global_seed(3141592653)
         a = np.array([1, 2, 3, 4, 5])
         dp_hist, _ = histogramdd(a, epsilon=10, bins=3, range=[(0, 10)], density=True)
 
@@ -52,7 +50,6 @@ def test_density_1d(self):
         self.assertAlmostEqual(dp_hist.sum(), 1.0 * 3 / 10)
 
     def test_density_2d(self):
-        global_seed(3141592653)
         a = np.array([[1, 2, 3, 4, 5], [6, 7, 8, 9, 10]]).T
         dp_hist, _ = histogramdd(a, epsilon=10, bins=3, range=[(0, 10), (0, 10)], density=True)
 
diff --git a/tests/tools/test_nanstd.py b/tests/tools/test_nanstd.py
index 4589528b..a300db79 100644
--- a/tests/tools/test_nanstd.py
+++ b/tests/tools/test_nanstd.py
@@ -3,7 +3,7 @@
 import numpy as np
 
 from diffprivlib.tools.utils import nanstd
-from diffprivlib.utils import PrivacyLeakWarning, global_seed, BudgetError
+from diffprivlib.utils import PrivacyLeakWarning, BudgetError, check_random_state
 
 
 class TestNanStd(TestCase):
@@ -38,7 +38,6 @@ def test_missing_bounds(self):
         self.assertIsNotNone(res)
 
     def test_large_epsilon(self):
-        global_seed(12345)
         a = np.random.random(1000)
         res = float(np.std(a))
         res_dp = nanstd(a, epsilon=5, bounds=(0, 1))
@@ -46,7 +45,6 @@ def test_large_epsilon(self):
         self.assertAlmostEqual(res, res_dp, delta=0.01)
 
     def test_large_epsilon_axis(self):
-        global_seed(12345)
         a = np.random.random((1000, 5))
         res = np.std(a, axis=0)
         res_dp = nanstd(a, epsilon=15, bounds=(0, 1), axis=0)
@@ -60,9 +58,10 @@ def test_array_like(self):
 
     def test_clipped_output(self):
         a = np.random.random((10,))
+        rng = check_random_state(0)
 
         for i in range(100):
-            self.assertTrue(0 <= nanstd(a, epsilon=1e-5, bounds=(0, 1)) <= 1)
+            self.assertTrue(0 <= nanstd(a, epsilon=1e-3, bounds=(0, 1), random_state=rng) <= 1)
 
     def test_nan(self):
         a = np.random.random((5, 5))
diff --git a/tests/tools/test_nanvar.py b/tests/tools/test_nanvar.py
index d429d7b9..ea8d23d8 100644
--- a/tests/tools/test_nanvar.py
+++ b/tests/tools/test_nanvar.py
@@ -3,7 +3,7 @@
 import numpy as np
 
 from diffprivlib.tools.utils import nanvar
-from diffprivlib.utils import PrivacyLeakWarning, BudgetError
+from diffprivlib.utils import PrivacyLeakWarning, BudgetError, check_random_state
 
 
 class TestNanVar(TestCase):
@@ -58,9 +58,10 @@ def test_array_like(self):
 
     def test_clipped_output(self):
         a = np.random.random((10,))
+        rng = check_random_state(0)
 
         for i in range(100):
-            self.assertTrue(0 <= nanvar(a, epsilon=1e-5, bounds=(0, 1)) <= 1)
+            self.assertTrue(0 <= nanvar(a, epsilon=1e-3, bounds=(0, 1), random_state=rng) <= 1)
 
     def test_nan(self):
         a = np.random.random((5, 5))
diff --git a/tests/tools/test_percentile.py b/tests/tools/test_percentile.py
index c32329db..1220fb35 100644
--- a/tests/tools/test_percentile.py
+++ b/tests/tools/test_percentile.py
@@ -4,6 +4,7 @@
 import pytest
 
 from diffprivlib.tools.quantiles import percentile
+from diffprivlib.utils import check_random_state
 
 
 class TestPercentile(TestCase):
@@ -19,9 +20,10 @@ def test_bad_percents(self):
         self.assertRaises(ValueError, percentile, a, [50] * 3 + [-1])
 
     def test_simple(self):
-        a = np.random.random(1000)
+        rng = check_random_state(0)
+        a = rng.random(1000)
 
-        res = percentile(a, 50, epsilon=5, bounds=(0, 1))
+        res = percentile(a, 50, epsilon=5, bounds=(0, 1), random_state=rng)
         self.assertAlmostEqual(res, 0.5, delta=0.05)
 
     @pytest.mark.filterwarnings("ignore:Bounds have not been specified")
diff --git a/tests/tools/test_quantile.py b/tests/tools/test_quantile.py
index 6d35092d..32f1da0b 100644
--- a/tests/tools/test_quantile.py
+++ b/tests/tools/test_quantile.py
@@ -4,7 +4,7 @@
 import pytest
 
 from diffprivlib.tools.quantiles import quantile
-from diffprivlib.utils import PrivacyLeakWarning, BudgetError
+from diffprivlib.utils import PrivacyLeakWarning, BudgetError, check_random_state
 
 
 class TestQuantile(TestCase):
@@ -47,9 +47,10 @@ def test_missing_bounds(self):
         self.assertIsNotNone(res)
 
     def test_large_epsilon(self):
-        a = np.random.random(1000)
+        rng = check_random_state(0)
+        a = rng.random(1000)
         res = np.quantile(a, 0.5)
-        res_dp = quantile(a, 0.5, epsilon=5, bounds=(0, 1))
+        res_dp = quantile(a, 0.5, epsilon=5, bounds=(0, 1), random_state=rng)
 
         self.assertAlmostEqual(float(res), float(res_dp), delta=0.01)
 
@@ -65,26 +66,29 @@ def test_output_type(self):
         self.assertTrue(isinstance(res, float))
 
     def test_simple(self):
-        a = np.random.random(1000)
+        rng = check_random_state(0)
+        a = rng.random(1000)
         for q in [0.1, 0.25, 0.5, 0.75, 0.9]:
-            res = quantile(a, q, epsilon=5, bounds=(0, 1))
+            res = quantile(a, q, epsilon=5, bounds=(0, 1), random_state=rng)
             self.assertAlmostEqual(res, q, delta=0.05)
 
     def test_normal(self):
-        a = np.random.normal(size=2000)
-        res = quantile(a, 0.5, epsilon=3, bounds=(-3, 3))
+        rng = check_random_state(0)
+        a = rng.normal(size=2000)
+        res = quantile(a, 0.5, epsilon=3, bounds=(-3, 3), random_state=rng)
         self.assertAlmostEqual(res, 0, delta=0.1)
 
     @pytest.mark.filterwarnings("ignore:Bounds have not been specified")
     def test_uniform_array(self):
         a = np.array([1] * 10)
-        res = quantile(a, 0.5, epsilon=1)
+        res = quantile(a, 0.5, epsilon=1, random_state=0)
         self.assertTrue(0 <= res <= 2)
 
     def test_multiple_q(self):
-        a = np.random.random(1000)
+        rng = check_random_state(0)
+        a = rng.random(1000)
         q = [0.1, 0.25, 0.5, 0.75, 0.9]
-        res = quantile(a, q, epsilon=5, bounds=(0, 1))
+        res = quantile(a, q, epsilon=5, bounds=(0, 1), random_state=rng)
 
         self.assertEqual(len(q), len(res))
         self.assertTrue(np.isclose(res, q, atol=0.05).all())
diff --git a/tests/tools/test_std.py b/tests/tools/test_std.py
index 0ab3d6a5..7a54d10f 100644
--- a/tests/tools/test_std.py
+++ b/tests/tools/test_std.py
@@ -3,7 +3,7 @@
 import numpy as np
 
 from diffprivlib.tools.utils import std
-from diffprivlib.utils import PrivacyLeakWarning, global_seed, BudgetError
+from diffprivlib.utils import PrivacyLeakWarning, BudgetError, check_random_state
 
 
 class TestStd(TestCase):
@@ -38,7 +38,6 @@ def test_missing_bounds(self):
         self.assertIsNotNone(res)
 
     def test_large_epsilon(self):
-        global_seed(12345)
         a = np.random.random(1000)
         res = float(np.std(a))
         res_dp = std(a, epsilon=5, bounds=(0, 1))
@@ -46,7 +45,6 @@ def test_large_epsilon(self):
         self.assertAlmostEqual(res, res_dp, delta=0.01)
 
     def test_large_epsilon_axis(self):
-        global_seed(12345)
         a = np.random.random((1000, 5))
         res = np.std(a, axis=0)
         res_dp = std(a, epsilon=15, bounds=(0, 1), axis=0)
@@ -60,9 +58,10 @@ def test_array_like(self):
 
     def test_clipped_output(self):
         a = np.random.random((10,))
+        rng = check_random_state(0)
 
         for i in range(100):
-            self.assertTrue(0 <= std(a, epsilon=1e-5, bounds=(0, 1)) <= 1)
+            self.assertTrue(0 <= std(a, epsilon=1e-3, bounds=(0, 1), random_state=rng) <= 1)
 
     def test_nan(self):
         a = np.random.random((5, 5))
diff --git a/tests/tools/test_var.py b/tests/tools/test_var.py
index f8e3218f..b2fccdb8 100644
--- a/tests/tools/test_var.py
+++ b/tests/tools/test_var.py
@@ -3,7 +3,7 @@
 import numpy as np
 
 from diffprivlib.tools.utils import var
-from diffprivlib.utils import PrivacyLeakWarning, BudgetError
+from diffprivlib.utils import PrivacyLeakWarning, BudgetError, check_random_state
 
 
 class TestVar(TestCase):
@@ -58,9 +58,10 @@ def test_array_like(self):
 
     def test_clipped_output(self):
         a = np.random.random((10,))
+        rng = check_random_state(0)
 
         for i in range(100):
-            self.assertTrue(0 <= var(a, epsilon=1e-5, bounds=(0, 1)) <= 1)
+            self.assertTrue(0 <= var(a, epsilon=1e-3, bounds=(0, 1), random_state=rng) <= 1)
 
     def test_nan(self):
         a = np.random.random((5, 5))