ENH: np.random.default_gen() (#13840)

* ENH: Rename seed_seq argument to seed and replace Generator() with default_gen()

10 files changed, 303 insertions, 203 deletions

diff --git a/numpy/random/__init__.py b/numpy/random/__init__.py
index 6543161f2..c7df6eb50 100644
--- a/numpy/random/__init__.py
+++ b/numpy/random/__init__.py
@@ -3,16 +3,42 @@
 Random Number Generation
 ========================
 
-Instantiate a BitGenerator and wrap it in a Generator
-which will convert the uniform stream to a number of distributions. The "bare"
-functions are kept for legacy code, they should be called with the newer API
-via ``np.random.Generator().function`` instead
+Use ``default_gen()`` to create a `Generator` and call its methods.
+
+=============== =========================================================
+Generator
+--------------- ---------------------------------------------------------
+Generator       Class implementing all of the random number distributions
+default_gen     Default constructor for ``Generator``
+=============== =========================================================
+
+============================================= ===
+BitGenerator Streams that work with Generator
+--------------------------------------------- ---
+MT19937
+PCG64
+Philox
+SFC64
+============================================= ===
+
+============================================= ===
+Getting entropy to initialize a BitGenerator
+--------------------------------------------- ---
+SeedSequence
+============================================= ===
+
+
+Legacy
+------
+
+For backwards compatibility with previous versions of numpy before 1.17, the
+various aliases to the global `RandomState` methods are left alone and do not
+use the new `Generator` API.
 
 ==================== =========================================================
 Utility functions
 -------------------- ---------------------------------------------------------
 random               Uniformly distributed floats over ``[0, 1)``
-integers             Uniformly distributed integers, replaces ``randint``
 bytes                Uniformly distributed random bytes.
 permutation          Randomly permute a sequence / generate a random sequence.
 shuffle              Randomly permute a sequence in place.
@@ -94,20 +120,6 @@ get_state            Get tuple representing internal state of generator.
 set_state            Set state of generator.
 ==================== =========================================================
 
-============================================= ===
-BitGenerator Streams that work with Generator
---------------------------------------------- ---
-MT19937
-PCG64
-Philox
-SFC64
-============================================= ===
-
-============================================= ===
-Getting entropy to initialize a BitGenerator
---------------------------------------------- ---
-SeedSequence
-============================================= ===
 
 """
 from __future__ import division, absolute_import, print_function
@@ -167,7 +179,7 @@ __all__ = [
 
 from . import mtrand
 from .mtrand import *
-from .generator import Generator
+from .generator import Generator, default_gen
 from .bit_generator import SeedSequence
 from .mt19937 import MT19937
 from .pcg64 import PCG64
@@ -176,7 +188,7 @@ from .sfc64 import SFC64
 from .mtrand import RandomState
 
 __all__ += ['Generator', 'RandomState', 'SeedSequence', 'MT19937',
-            'Philox', 'PCG64', 'SFC64']
+            'Philox', 'PCG64', 'SFC64', 'default_gen']
 
 
 def __RandomState_ctor():
diff --git a/numpy/random/bit_generator.pxd b/numpy/random/bit_generator.pxd
index 846e3d30a..79fe69275 100644
--- a/numpy/random/bit_generator.pxd
+++ b/numpy/random/bit_generator.pxd
@@ -13,7 +13,6 @@ cdef class BitGenerator():
 
 cdef class SeedSequence():
     cdef readonly object entropy
-    cdef readonly object program_entropy
     cdef readonly tuple spawn_key
     cdef readonly int pool_size
     cdef readonly object pool
diff --git a/numpy/random/bit_generator.pyx b/numpy/random/bit_generator.pyx
index 8886d2e05..cd1b628eb 100644
--- a/numpy/random/bit_generator.pyx
+++ b/numpy/random/bit_generator.pyx
@@ -116,7 +116,7 @@ def _coerce_to_uint32_array(x):
     Examples
     --------
     >>> import numpy as np
-    >>> from seed_seq import _coerce_to_uint32_array
+    >>> from np.random.bit_generator import _coerce_to_uint32_array
     >>> _coerce_to_uint32_array(12345)
     array([12345], dtype=uint32)
     >>> _coerce_to_uint32_array('12345')
@@ -473,17 +473,20 @@ ISpawnableSeedSequence.register(SeedSequence)
 
 cdef class BitGenerator():
     """
-    BitGenerator(seed_seq=None)
+    BitGenerator(seed=None)
 
     Base Class for generic BitGenerators, which provide a stream
     of random bits based on different algorithms. Must be overridden.
 
     Parameters
     ----------
-    seed_seq : {None, ISeedSequence, int, sequence[int]}, optional
-        A ISeedSequence to initialize the BitGenerator. If None, one will be
-        created. If an int or a sequence of ints, it will be used as the
-        entropy for creating a SeedSequence.
+    seed : {None, int, array_like[ints], ISeedSequence}, optional
+        A seed to initialize the `BitGenerator`. If None, then fresh,
+        unpredictable entropy will be pulled from the OS. If an ``int`` or
+        ``array_like[ints]`` is passed, then it will be passed to
+        `SeedSequence` to derive the initial `BitGenerator` state. One may also
+        pass in an implementor of the `ISeedSequence` interface like
+        `SeedSequence`.
 
     Attributes
     ----------
@@ -498,7 +501,7 @@ cdef class BitGenerator():
     SeedSequence
     """
 
-    def __init__(self, seed_seq=None):
+    def __init__(self, seed=None):
         self.lock = Lock()
         self._bitgen.state = <void *>0
         if type(self) is BitGenerator:
@@ -509,9 +512,9 @@ cdef class BitGenerator():
 
         cdef const char *name = "BitGenerator"
         self.capsule = PyCapsule_New(<void *>&self._bitgen, name, NULL)
-        if not isinstance(seed_seq, ISeedSequence):
-            seed_seq = SeedSequence(seed_seq)
-        self._seed_seq = seed_seq
+        if not isinstance(seed, ISeedSequence):
+            seed = SeedSequence(seed)
+        self._seed_seq = seed
 
     # Pickling support:
     def __getstate__(self):
diff --git a/numpy/random/generator.pyx b/numpy/random/generator.pyx
index 05323c422..e6b70f176 100644
--- a/numpy/random/generator.pyx
+++ b/numpy/random/generator.pyx
@@ -34,9 +34,18 @@ __all__ = ['Generator', 'beta', 'binomial', 'bytes', 'chisquare', 'choice',
 np.import_array()
 
 
+cdef bint _check_bit_generator(object bitgen):
+    """Check if an object satisfies the BitGenerator interface.
+    """
+    if not hasattr(bitgen, "capsule"):
+        return False
+    cdef const char *name = "BitGenerator"
+    return PyCapsule_IsValid(bitgen.capsule, name)
+
+
 cdef class Generator:
     """
-    Generator(bit_generator=None)
+    Generator(bit_generator)
 
     Container for the BitGenerators.
 
@@ -48,6 +57,9 @@ cdef class Generator:
     array filled with generated values is returned. If `size` is a tuple,
     then an array with that shape is filled and returned.
 
+    The function :func:`numpy.random.default_gen` will instantiate
+    a `Generator` with numpy's default `BitGenerator`.
+
     **No Compatibility Guarantee**
 
     ``Generator`` does not provide a version compatibility guarantee. In
@@ -55,9 +67,8 @@ cdef class Generator:
 
     Parameters
     ----------
-    bit_generator : BitGenerator, optional
-        BitGenerator to use as the core generator. If none is provided, uses
-        Xoshiro256.
+    bit_generator : BitGenerator
+        BitGenerator to use as the core generator.
 
     Notes
     -----
@@ -70,16 +81,14 @@ cdef class Generator:
 
     Examples
     --------
-    >>> from numpy.random import Generator
-    >>> rg = Generator()
+    >>> from numpy.random import Generator, PCG64
+    >>> rg = Generator(PCG64())
     >>> rg.standard_normal()
     -0.203  # random
 
-    Using a specific generator
-
-    >>> from numpy.random import MT19937
-    >>> rg = Generator(MT19937())
-
+    See Also
+    --------
+    default_gen : Recommended constructor for `Generator`.
     """
     cdef public object _bit_generator
     cdef bitgen_t _bitgen
@@ -87,9 +96,7 @@ cdef class Generator:
     cdef object lock
     _poisson_lam_max = POISSON_LAM_MAX
 
-    def __init__(self, bit_generator=None):
-        if bit_generator is None:
-            bit_generator = PCG64()
+    def __init__(self, bit_generator):
         self._bit_generator = bit_generator
 
         capsule = bit_generator.capsule
@@ -166,16 +173,17 @@ cdef class Generator:
 
         Examples
         --------
-        >>> np.random.Generator().random()
+        >>> rng = np.random.default_gen()
+        >>> rng.random()
         0.47108547995356098 # random
-        >>> type(np.random.Generator().random())
+        >>> type(rng.random())
         <class 'float'>
-        >>> np.random.Generator().random((5,))
+        >>> rng.random((5,))
         array([ 0.30220482,  0.86820401,  0.1654503 ,  0.11659149,  0.54323428]) # random
 
         Three-by-two array of random numbers from [-5, 0):
 
-        >>> 5 * np.random.Generator().random((3, 2)) - 5
+        >>> 5 * rng.random((3, 2)) - 5
         array([[-3.99149989, -0.52338984], # random
                [-2.99091858, -0.79479508],
                [-1.23204345, -1.75224494]])
@@ -321,7 +329,7 @@ cdef class Generator:
         --------
         Output a 3x8000 array:
 
-        >>> n = np.random.Generator().standard_exponential((3, 8000))
+        >>> n = np.random.default_gen().standard_exponential((3, 8000))
 
         """
         key = np.dtype(dtype).name
@@ -389,30 +397,31 @@ cdef class Generator:
 
         Examples
         --------
-        >>> np.random.Generator().integers(2, size=10)
+        >>> rng = np.random.default_gen()
+        >>> rng.integers(2, size=10)
         array([1, 0, 0, 0, 1, 1, 0, 0, 1, 0])  # random
-        >>> np.random.Generator().integers(1, size=10)
+        >>> rng.integers(1, size=10)
         array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
 
         Generate a 2 x 4 array of ints between 0 and 4, inclusive:
 
-        >>> np.random.Generator().integers(5, size=(2, 4))
+        >>> rng.integers(5, size=(2, 4))
         array([[4, 0, 2, 1],
                [3, 2, 2, 0]])  # random
 
         Generate a 1 x 3 array with 3 different upper bounds
 
-        >>> np.random.Generator().integers(1, [3, 5, 10])
+        >>> rng.integers(1, [3, 5, 10])
         array([2, 2, 9])  # random
 
         Generate a 1 by 3 array with 3 different lower bounds
 
-        >>> np.random.Generator().integers([1, 5, 7], 10)
+        >>> rng.integers([1, 5, 7], 10)
         array([9, 8, 7])  # random
 
         Generate a 2 by 4 array using broadcasting with dtype of uint8
 
-        >>> np.random.Generator().integers([1, 3, 5, 7], [[10], [20]], dtype=np.uint8)
+        >>> rng.integers([1, 3, 5, 7], [[10], [20]], dtype=np.uint8)
         array([[ 8,  6,  9,  7],
                [ 1, 16,  9, 12]], dtype=uint8)  # random
 
@@ -484,7 +493,7 @@ cdef class Generator:
 
         Examples
         --------
-        >>> np.random.Generator().bytes(10)
+        >>> np.random.default_gen().bytes(10)
         ' eh\\x85\\x022SZ\\xbf\\xa4' #random
 
         """
@@ -550,33 +559,34 @@ cdef class Generator:
         --------
         Generate a uniform random sample from np.arange(5) of size 3:
 
-        >>> np.random.Generator().choice(5, 3)
+        >>> rng = np.random.default_gen()
+        >>> rng.choice(5, 3)
         array([0, 3, 4]) # random
-        >>> #This is equivalent to np.random.Generator().integers(0,5,3)
+        >>> #This is equivalent to rng.integers(0,5,3)
 
         Generate a non-uniform random sample from np.arange(5) of size 3:
 
-        >>> np.random.Generator().choice(5, 3, p=[0.1, 0, 0.3, 0.6, 0])
+        >>> rng.choice(5, 3, p=[0.1, 0, 0.3, 0.6, 0])
         array([3, 3, 0]) # random
 
         Generate a uniform random sample from np.arange(5) of size 3 without
         replacement:
 
-        >>> np.random.Generator().choice(5, 3, replace=False)
+        >>> rng.choice(5, 3, replace=False)
         array([3,1,0]) # random
-        >>> #This is equivalent to np.random.Generator().permutation(np.arange(5))[:3]
+        >>> #This is equivalent to rng.permutation(np.arange(5))[:3]
 
         Generate a non-uniform random sample from np.arange(5) of size
         3 without replacement:
 
-        >>> np.random.Generator().choice(5, 3, replace=False, p=[0.1, 0, 0.3, 0.6, 0])
+        >>> rng.choice(5, 3, replace=False, p=[0.1, 0, 0.3, 0.6, 0])
         array([2, 3, 0]) # random
 
         Any of the above can be repeated with an arbitrary array-like
         instead of just integers. For instance:
 
         >>> aa_milne_arr = ['pooh', 'rabbit', 'piglet', 'Christopher']
-        >>> np.random.Generator().choice(aa_milne_arr, 5, p=[0.5, 0.1, 0.1, 0.3])
+        >>> rng.choice(aa_milne_arr, 5, p=[0.5, 0.1, 0.1, 0.3])
         array(['pooh', 'pooh', 'pooh', 'Christopher', 'piglet'], # random
               dtype='<U11')
 
@@ -784,7 +794,7 @@ cdef class Generator:
         --------
         Draw samples from the distribution:
 
-        >>> s = np.random.Generator().uniform(-1,0,1000)
+        >>> s = np.random.default_gen().uniform(-1,0,1000)
 
         All values are within the given interval:
 
@@ -869,8 +879,8 @@ cdef class Generator:
         -----
         For random samples from :math:`N(\\mu, \\sigma^2)`, use one of::
 
-            mu + sigma * np.random.Generator().standard_normal(size=...)
-            np.random.Generator().normal(mu, sigma, size=...)
+            mu + sigma * gen.standard_normal(size=...)
+            gen.normal(mu, sigma, size=...)
 
         See Also
         --------
@@ -880,22 +890,23 @@ cdef class Generator:
 
         Examples
         --------
-        >>> np.random.Generator().standard_normal()
+        >>> rng = np.random.default_gen()
+        >>> rng.standard_normal()
         2.1923875335537315 #random
 
-        >>> s = np.random.Generator().standard_normal(8000)
+        >>> s = rng.standard_normal(8000)
         >>> s
         array([ 0.6888893 ,  0.78096262, -0.89086505, ...,  0.49876311,  # random
                -0.38672696, -0.4685006 ])                                # random
         >>> s.shape
         (8000,)
-        >>> s = np.random.Generator().standard_normal(size=(3, 4, 2))
+        >>> s = rng.standard_normal(size=(3, 4, 2))
         >>> s.shape
         (3, 4, 2)
 
         Two-by-four array of samples from :math:`N(3, 6.25)`:
 
-        >>> 3 + 2.5 * np.random.Generator().standard_normal(size=(2, 4))
+        >>> 3 + 2.5 * rng.standard_normal(size=(2, 4))
         array([[-4.49401501,  4.00950034, -1.81814867,  7.29718677],   # random
                [ 0.39924804,  4.68456316,  4.99394529,  4.84057254]])  # random
 
@@ -962,8 +973,8 @@ cdef class Generator:
         The function has its peak at the mean, and its "spread" increases with
         the standard deviation (the function reaches 0.607 times its maximum at
         :math:`x + \\sigma` and :math:`x - \\sigma` [2]_).  This implies that
-        `numpy.random.Generator().normal` is more likely to return
-        samples lying close to the mean, rather than those far away.
+        :meth:`normal` is more likely to return samples lying close to the
+        mean, rather than those far away.
 
         References
         ----------
@@ -978,7 +989,7 @@ cdef class Generator:
         Draw samples from the distribution:
 
         >>> mu, sigma = 0, 0.1 # mean and standard deviation
-        >>> s = np.random.Generator().normal(mu, sigma, 1000)
+        >>> s = np.random.default_gen().normal(mu, sigma, 1000)
 
         Verify the mean and the variance:
 
@@ -1000,7 +1011,7 @@ cdef class Generator:
 
         Two-by-four array of samples from N(3, 6.25):
 
-        >>> np.random.Generator().normal(3, 2.5, size=(2, 4))
+        >>> np.random.default_gen().normal(3, 2.5, size=(2, 4))
         array([[-4.49401501,  4.00950034, -1.81814867,  7.29718677],   # random
                [ 0.39924804,  4.68456316,  4.99394529,  4.84057254]])  # random
 
@@ -1074,7 +1085,7 @@ cdef class Generator:
         Draw samples from the distribution:
 
         >>> shape, scale = 2., 1. # mean and width
-        >>> s = np.random.Generator().standard_gamma(shape, 1000000)
+        >>> s = np.random.default_gen().standard_gamma(shape, 1000000)
 
         Display the histogram of the samples, along with
         the probability density function:
@@ -1162,7 +1173,7 @@ cdef class Generator:
         Draw samples from the distribution:
 
         >>> shape, scale = 2., 2.  # mean=4, std=2*sqrt(2)
-        >>> s = np.random.Generator().gamma(shape, scale, 1000)
+        >>> s = np.random.default_gen().gamma(shape, scale, 1000)
 
         Display the histogram of the samples, along with
         the probability density function:
@@ -1252,7 +1263,7 @@ cdef class Generator:
 
         >>> dfnum = 1. # between group degrees of freedom
         >>> dfden = 48. # within groups degrees of freedom
-        >>> s = np.random.Generator().f(dfnum, dfden, 1000)
+        >>> s = np.random.default_gen().f(dfnum, dfden, 1000)
 
         The lower bound for the top 1% of the samples is :
 
@@ -1328,12 +1339,13 @@ cdef class Generator:
         distribution for the null hypothesis.  We'll plot the two probability
         distributions for comparison.
 
+        >>> rng = np.random.default_gen()
         >>> dfnum = 3 # between group deg of freedom
         >>> dfden = 20 # within groups degrees of freedom
         >>> nonc = 3.0
-        >>> nc_vals = np.random.Generator().noncentral_f(dfnum, dfden, nonc, 1000000)
+        >>> nc_vals = rng.noncentral_f(dfnum, dfden, nonc, 1000000)
         >>> NF = np.histogram(nc_vals, bins=50, density=True)
-        >>> c_vals = np.random.Generator().f(dfnum, dfden, 1000000)
+        >>> c_vals = rng.f(dfnum, dfden, 1000000)
         >>> F = np.histogram(c_vals, bins=50, density=True)
         >>> import matplotlib.pyplot as plt
         >>> plt.plot(F[1][1:], F[0])
@@ -1405,7 +1417,7 @@ cdef class Generator:
 
         Examples
         --------
-        >>> np.random.Generator().chisquare(2,4)
+        >>> np.random.default_gen().chisquare(2,4)
         array([ 1.89920014,  9.00867716,  3.13710533,  5.62318272]) # random
 
         """
@@ -1463,8 +1475,9 @@ cdef class Generator:
         --------
         Draw values from the distribution and plot the histogram
 
+        >>> rng = np.random.default_gen()
         >>> import matplotlib.pyplot as plt
-        >>> values = plt.hist(np.random.Generator().noncentral_chisquare(3, 20, 100000),
+        >>> values = plt.hist(rng.noncentral_chisquare(3, 20, 100000),
         ...                   bins=200, density=True)
         >>> plt.show()
 
@@ -1472,9 +1485,9 @@ cdef class Generator:
         and compare to a chisquare.
 
         >>> plt.figure()
-        >>> values = plt.hist(np.random.Generator().noncentral_chisquare(3, .0000001, 100000),
+        >>> values = plt.hist(rng.noncentral_chisquare(3, .0000001, 100000),
         ...                   bins=np.arange(0., 25, .1), density=True)
-        >>> values2 = plt.hist(np.random.Generator().chisquare(3, 100000),
+        >>> values2 = plt.hist(rng.chisquare(3, 100000),
         ...                    bins=np.arange(0., 25, .1), density=True)
         >>> plt.plot(values[1][0:-1], values[0]-values2[0], 'ob')
         >>> plt.show()
@@ -1483,7 +1496,7 @@ cdef class Generator:
         distribution.
 
         >>> plt.figure()
-        >>> values = plt.hist(np.random.Generator().noncentral_chisquare(3, 20, 100000),
+        >>> values = plt.hist(rng.noncentral_chisquare(3, 20, 100000),
         ...                   bins=200, density=True)
         >>> plt.show()
 
@@ -1549,7 +1562,7 @@ cdef class Generator:
         Draw samples and plot the distribution:
 
         >>> import matplotlib.pyplot as plt
-        >>> s = np.random.Generator().standard_cauchy(1000000)
+        >>> s = np.random.default_gen().standard_cauchy(1000000)
         >>> s = s[(s>-25) & (s<25)]  # truncate distribution so it plots well
         >>> plt.hist(s, bins=100)
         >>> plt.show()
@@ -1622,7 +1635,7 @@ cdef class Generator:
         We have 10 degrees of freedom, so is the sample mean within 95% of the
         recommended value?
 
-        >>> s = np.random.Generator().standard_t(10, size=100000)
+        >>> s = np.random.default_gen().standard_t(10, size=100000)
         >>> np.mean(intake)
         6753.636363636364
         >>> intake.std(ddof=1)
@@ -1716,7 +1729,7 @@ cdef class Generator:
         Draw samples from the distribution:
 
         >>> mu, kappa = 0.0, 4.0 # mean and dispersion
-        >>> s = np.random.Generator().vonmises(mu, kappa, 1000)
+        >>> s = np.random.default_gen().vonmises(mu, kappa, 1000)
 
         Display the histogram of the samples, along with
         the probability density function:
@@ -1816,7 +1829,7 @@ cdef class Generator:
         Draw samples from the distribution:
 
         >>> a, m = 3., 2.  # shape and mode
-        >>> s = (np.random.Generator().pareto(a, 1000) + 1) * m
+        >>> s = (np.random.default_gen().pareto(a, 1000) + 1) * m
 
         Display the histogram of the samples, along with the probability
         density function:
@@ -1908,8 +1921,9 @@ cdef class Generator:
         --------
         Draw samples from the distribution:
 
+        >>> rng = np.random.default_gen()
         >>> a = 5. # shape
-        >>> s = np.random.Generator().weibull(a, 1000)
+        >>> s = rng.weibull(a, 1000)
 
         Display the histogram of the samples, along with
         the probability density function:
@@ -1919,7 +1933,7 @@ cdef class Generator:
         >>> def weib(x,n,a):
         ...     return (a / n) * (x / n)**(a - 1) * np.exp(-(x / n)**a)
 
-        >>> count, bins, ignored = plt.hist(np.random.Generator().weibull(5.,1000))
+        >>> count, bins, ignored = plt.hist(rng.weibull(5.,1000))
         >>> x = np.arange(1,100.)/50.
         >>> scale = count.max()/weib(x, 1., 5.).max()
         >>> plt.plot(x, weib(x, 1., 5.)*scale)
@@ -1987,9 +2001,10 @@ cdef class Generator:
         --------
         Draw samples from the distribution:
 
+        >>> rng = np.random.default_gen()
         >>> a = 5. # shape
         >>> samples = 1000
-        >>> s = np.random.Generator().power(a, samples)
+        >>> s = rng.power(a, samples)
 
         Display the histogram of the samples, along with
         the probability density function:
@@ -2005,20 +2020,20 @@ cdef class Generator:
         Compare the power function distribution to the inverse of the Pareto.
 
         >>> from scipy import stats  # doctest: +SKIP
-        >>> rvs = np.random.Generator().power(5, 1000000)
-        >>> rvsp = np.random.Generator().pareto(5, 1000000)
+        >>> rvs = rng.power(5, 1000000)
+        >>> rvsp = rng.pareto(5, 1000000)
         >>> xx = np.linspace(0,1,100)
         >>> powpdf = stats.powerlaw.pdf(xx,5)  # doctest: +SKIP
 
         >>> plt.figure()
         >>> plt.hist(rvs, bins=50, density=True)
         >>> plt.plot(xx,powpdf,'r-')  # doctest: +SKIP
-        >>> plt.title('np.random.Generator().power(5)')
+        >>> plt.title('power(5)')
 
         >>> plt.figure()
         >>> plt.hist(1./(1.+rvsp), bins=50, density=True)
         >>> plt.plot(xx,powpdf,'r-')  # doctest: +SKIP
-        >>> plt.title('inverse of 1 + np.random.Generator().pareto(5)')
+        >>> plt.title('inverse of 1 + Generator.pareto(5)')
 
         >>> plt.figure()
         >>> plt.hist(1./(1.+rvsp), bins=50, density=True)
@@ -2093,7 +2108,7 @@ cdef class Generator:
         Draw samples from the distribution
 
         >>> loc, scale = 0., 1.
-        >>> s = np.random.Generator().laplace(loc, scale, 1000)
+        >>> s = np.random.default_gen().laplace(loc, scale, 1000)
 
         Display the histogram of the samples, along with
         the probability density function:
@@ -2195,8 +2210,9 @@ cdef class Generator:
         --------
         Draw samples from the distribution:
 
+        >>> rng = np.random.default_gen()
         >>> mu, beta = 0, 0.1 # location and scale
-        >>> s = np.random.Generator().gumbel(mu, beta, 1000)
+        >>> s = rng.gumbel(mu, beta, 1000)
 
         Display the histogram of the samples, along with
         the probability density function:
@@ -2214,7 +2230,7 @@ cdef class Generator:
         >>> means = []
         >>> maxima = []
         >>> for i in range(0,1000) :
-        ...    a = np.random.Generator().normal(mu, beta, 1000)
+        ...    a = rng.normal(mu, beta, 1000)
         ...    means.append(a.mean())
         ...    maxima.append(a.max())
         >>> count, bins, ignored = plt.hist(maxima, 30, density=True)
@@ -2296,7 +2312,7 @@ cdef class Generator:
         Draw samples from the distribution:
 
         >>> loc, scale = 10, 1
-        >>> s = np.random.Generator().logistic(loc, scale, 10000)
+        >>> s = np.random.default_gen().logistic(loc, scale, 10000)
         >>> import matplotlib.pyplot as plt
         >>> count, bins, ignored = plt.hist(s, bins=50)
 
@@ -2378,8 +2394,9 @@ cdef class Generator:
         --------
         Draw samples from the distribution:
 
+        >>> rng = np.random.default_gen()
         >>> mu, sigma = 3., 1. # mean and standard deviation
-        >>> s = np.random.Generator().lognormal(mu, sigma, 1000)
+        >>> s = rng.lognormal(mu, sigma, 1000)
 
         Display the histogram of the samples, along with
         the probability density function:
@@ -2401,9 +2418,10 @@ cdef class Generator:
 
         >>> # Generate a thousand samples: each is the product of 100 random
         >>> # values, drawn from a normal distribution.
+        >>> rng = rng
         >>> b = []
         >>> for i in range(1000):
-        ...    a = 10. + np.random.Generator().standard_normal(100)
+        ...    a = 10. + rng.standard_normal(100)
         ...    b.append(np.product(a))
 
         >>> b = np.array(b) / np.min(b) # scale values to be positive
@@ -2471,7 +2489,8 @@ cdef class Generator:
         Draw values from the distribution and plot the histogram
 
         >>> from matplotlib.pyplot import hist
-        >>> values = hist(np.random.Generator().rayleigh(3, 100000), bins=200, density=True)
+        >>> rng = np.random.default_gen()
+        >>> values = hist(rng.rayleigh(3, 100000), bins=200, density=True)
 
         Wave heights tend to follow a Rayleigh distribution. If the mean wave
         height is 1 meter, what fraction of waves are likely to be larger than 3
@@ -2479,7 +2498,7 @@ cdef class Generator:
 
         >>> meanvalue = 1
         >>> modevalue = np.sqrt(2 / np.pi) * meanvalue
-        >>> s = np.random.Generator().rayleigh(modevalue, 1000000)
+        >>> s = rng.rayleigh(modevalue, 1000000)
 
         The percentage of waves larger than 3 meters is:
 
@@ -2551,7 +2570,7 @@ cdef class Generator:
         Draw values from the distribution and plot the histogram:
 
         >>> import matplotlib.pyplot as plt
-        >>> h = plt.hist(np.random.Generator().wald(3, 2, 100000), bins=200, density=True)
+        >>> h = plt.hist(np.random.default_gen().wald(3, 2, 100000), bins=200, density=True)
         >>> plt.show()
 
         """
@@ -2618,7 +2637,7 @@ cdef class Generator:
         Draw values from the distribution and plot the histogram:
 
         >>> import matplotlib.pyplot as plt
-        >>> h = plt.hist(np.random.Generator().triangular(-3, 0, 8, 100000), bins=200,
+        >>> h = plt.hist(np.random.default_gen().triangular(-3, 0, 8, 100000), bins=200,
         ...              density=True)
         >>> plt.show()
 
@@ -2730,8 +2749,9 @@ cdef class Generator:
         --------
         Draw samples from the distribution:
 
+        >>> rng = np.random.default_gen()
         >>> n, p = 10, .5  # number of trials, probability of each trial
-        >>> s = np.random.Generator().binomial(n, p, 1000)
+        >>> s = rng.binomial(n, p, 1000)
         # result of flipping a coin 10 times, tested 1000 times.
 
         A real world example. A company drills 9 wild-cat oil exploration
@@ -2741,7 +2761,7 @@ cdef class Generator:
         Let's do 20,000 trials of the model, and count the number that
         generate zero positive results.
 
-        >>> sum(np.random.Generator().binomial(9, 0.1, 20000) == 0)/20000.
+        >>> sum(rng.binomial(9, 0.1, 20000) == 0)/20000.
         # answer = 0.38885, or 38%.
 
         """
@@ -2868,7 +2888,7 @@ cdef class Generator:
         for each successive well, that is what is the probability of a
         single success after drilling 5 wells, after 6 wells, etc.?
 
-        >>> s = np.random.Generator().negative_binomial(1, 0.1, 100000)
+        >>> s = np.random.default_gen().negative_binomial(1, 0.1, 100000)
         >>> for i in range(1, 11): # doctest: +SKIP
         ...    probability = sum(s<i) / 100000.
         ...    print(i, "wells drilled, probability of one success =", probability)
@@ -2932,7 +2952,8 @@ cdef class Generator:
         Draw samples from the distribution:
 
         >>> import numpy as np
-        >>> s = np.random.Generator().poisson(5, 10000)
+        >>> rng = np.random.default_gen()
+        >>> s = rng.poisson(5, 10000)
 
         Display histogram of the sample:
 
@@ -2942,7 +2963,7 @@ cdef class Generator:
 
         Draw each 100 values for lambda 100 and 500:
 
-        >>> s = np.random.Generator().poisson(lam=(100., 500.), size=(100, 2))
+        >>> s = rng.poisson(lam=(100., 500.), size=(100, 2))
 
         """
         return disc(&random_poisson, &self._bitgen, size, self.lock, 1, 0,
@@ -3007,7 +3028,7 @@ cdef class Generator:
         Draw samples from the distribution:
 
         >>> a = 2. # parameter
-        >>> s = np.random.Generator().zipf(a, 1000)
+        >>> s = np.random.default_gen().zipf(a, 1000)
 
         Display the histogram of the samples, along with
         the probability density function:
@@ -3068,7 +3089,7 @@ cdef class Generator:
         Draw ten thousand values from the geometric distribution,
         with the probability of an individual success equal to 0.35:
 
-        >>> z = np.random.Generator().geometric(p=0.35, size=10000)
+        >>> z = np.random.default_gen().geometric(p=0.35, size=10000)
 
         How many trials succeeded after a single run?
 
@@ -3168,9 +3189,10 @@ cdef class Generator:
         --------
         Draw samples from the distribution:
 
+        >>> rng = np.random.default_gen()
         >>> ngood, nbad, nsamp = 100, 2, 10
         # number of good, number of bad, and number of samples
-        >>> s = np.random.Generator().hypergeometric(ngood, nbad, nsamp, 1000)
+        >>> s = rng.hypergeometric(ngood, nbad, nsamp, 1000)
         >>> from matplotlib.pyplot import hist
         >>> hist(s)
         #   note that it is very unlikely to grab both bad items
@@ -3179,7 +3201,7 @@ cdef class Generator:
         If you pull 15 marbles at random, how likely is it that
         12 or more of them are one color?
 
-        >>> s = np.random.Generator().hypergeometric(15, 15, 15, 100000)
+        >>> s = rng.hypergeometric(15, 15, 15, 100000)
         >>> sum(s>=12)/100000. + sum(s<=3)/100000.
         #   answer = 0.003 ... pretty unlikely!
 
@@ -3283,7 +3305,7 @@ cdef class Generator:
         Draw samples from the distribution:
 
         >>> a = .6
-        >>> s = np.random.Generator().logseries(a, 10000)
+        >>> s = np.random.default_gen().logseries(a, 10000)
         >>> import matplotlib.pyplot as plt
         >>> count, bins, ignored = plt.hist(s)
 
@@ -3373,7 +3395,7 @@ cdef class Generator:
         Diagonal covariance means that points are oriented along x or y-axis:
 
         >>> import matplotlib.pyplot as plt
-        >>> x, y = np.random.Generator().multivariate_normal(mean, cov, 5000).T
+        >>> x, y = np.random.default_gen().multivariate_normal(mean, cov, 5000).T
         >>> plt.plot(x, y, 'x')
         >>> plt.axis('equal')
         >>> plt.show()
@@ -3393,7 +3415,7 @@ cdef class Generator:
         --------
         >>> mean = (1, 2)
         >>> cov = [[1, 0], [0, 1]]
-        >>> x = np.random.Generator().multivariate_normal(mean, cov, (3, 3))
+        >>> x = np.random.default_gen().multivariate_normal(mean, cov, (3, 3))
         >>> x.shape
         (3, 3, 2)
 
@@ -3507,14 +3529,15 @@ cdef class Generator:
         --------
         Throw a dice 20 times:
 
-        >>> np.random.Generator().multinomial(20, [1/6.]*6, size=1)
+        >>> rng = np.random.default_gen()
+        >>> rng.multinomial(20, [1/6.]*6, size=1)
         array([[4, 1, 7, 5, 2, 1]])  # random
 
         It landed 4 times on 1, once on 2, etc.
 
         Now, throw the dice 20 times, and 20 times again:
 
-        >>> np.random.Generator().multinomial(20, [1/6.]*6, size=2)
+        >>> rng.multinomial(20, [1/6.]*6, size=2)
         array([[3, 4, 3, 3, 4, 3],
                [2, 4, 3, 4, 0, 7]])  # random
 
@@ -3524,7 +3547,7 @@ cdef class Generator:
         Now, do one experiment throwing the dice 10 time, and 10 times again,
         and another throwing the dice 20 times, and 20 times again:
 
-        >>> np.random.Generator().multinomial([[10], [20]], [1/6.]*6, size=2)
+        >>> rng.multinomial([[10], [20]], [1/6.]*6, size=2)
         array([[[2, 4, 0, 1, 2, 1],
                 [1, 3, 0, 3, 1, 2]],
                [[1, 4, 4, 4, 4, 3],
@@ -3535,7 +3558,7 @@ cdef class Generator:
 
         A loaded die is more likely to land on number 6:
 
-        >>> np.random.Generator().multinomial(100, [1/7.]*5 + [2/7.])
+        >>> rng.multinomial(100, [1/7.]*5 + [2/7.])
         array([11, 16, 14, 17, 16, 26])  # random
 
         The probability inputs should be normalized. As an implementation
@@ -3544,12 +3567,12 @@ cdef class Generator:
         A biased coin which has twice as much weight on one side as on the
         other should be sampled like so:
 
-        >>> np.random.Generator().multinomial(100, [1.0 / 3, 2.0 / 3])  # RIGHT
+        >>> rng.multinomial(100, [1.0 / 3, 2.0 / 3])  # RIGHT
         array([38, 62])  # random
 
         not like:
 
-        >>> np.random.Generator().multinomial(100, [1.0, 2.0])  # WRONG
+        >>> rng.multinomial(100, [1.0, 2.0])  # WRONG
         Traceback (most recent call last):
         ValueError: pvals < 0, pvals > 1 or pvals contains NaNs
 
@@ -3683,7 +3706,7 @@ cdef class Generator:
         average length, but allowing some variation in the relative sizes of
         the pieces.
 
-        >>> s = np.random.Generator().dirichlet((10, 5, 3), 20).transpose()
+        >>> s = np.random.default_gen().dirichlet((10, 5, 3), 20).transpose()
 
         >>> import matplotlib.pyplot as plt
         >>> plt.barh(range(20), s[0])
@@ -3775,15 +3798,16 @@ cdef class Generator:
 
         Examples
         --------
+        >>> rng = np.random.default_gen()
         >>> arr = np.arange(10)
-        >>> np.random.Generator().shuffle(arr)
+        >>> rng.shuffle(arr)
         >>> arr
         [1 7 5 2 9 4 3 6 0 8] # random
 
         Multi-dimensional arrays are only shuffled along the first axis:
 
         >>> arr = np.arange(9).reshape((3, 3))
-        >>> np.random.Generator().shuffle(arr)
+        >>> rng.shuffle(arr)
         >>> arr
         array([[3, 4, 5], # random
                [6, 7, 8],
@@ -3885,14 +3909,15 @@ cdef class Generator:
 
         Examples
         --------
-        >>> np.random.Generator().permutation(10)
+        >>> rng = np.random.default_gen()
+        >>> rng.permutation(10)
         array([1, 7, 4, 3, 0, 9, 2, 5, 8, 6]) # random
 
-        >>> np.random.Generator().permutation([1, 4, 9, 12, 15])
+        >>> rng.permutation([1, 4, 9, 12, 15])
         array([15,  1,  9,  4, 12]) # random
 
         >>> arr = np.arange(9).reshape((3, 3))
-        >>> np.random.Generator().permutation(arr)
+        >>> rng.permutation(arr)
         array([[6, 7, 8], # random
                [0, 1, 2],
                [3, 4, 5]])
@@ -3918,47 +3943,34 @@ cdef class Generator:
         self.shuffle(idx)
         return arr[idx]
 
-_random_generator = Generator()
-
-beta = _random_generator.beta
-binomial = _random_generator.binomial
-bytes = _random_generator.bytes
-chisquare = _random_generator.chisquare
-choice = _random_generator.choice
-dirichlet = _random_generator.dirichlet
-exponential = _random_generator.exponential
-f = _random_generator.f
-gamma = _random_generator.gamma
-geometric = _random_generator.geometric
-gumbel = _random_generator.gumbel
-hypergeometric = _random_generator.hypergeometric
-integers = _random_generator.integers
-laplace = _random_generator.laplace
-logistic = _random_generator.logistic
-lognormal = _random_generator.lognormal
-logseries = _random_generator.logseries
-multinomial = _random_generator.multinomial
-multivariate_normal = _random_generator.multivariate_normal
-negative_binomial = _random_generator.negative_binomial
-noncentral_chisquare = _random_generator.noncentral_chisquare
-noncentral_f = _random_generator.noncentral_f
-normal = _random_generator.normal
-pareto = _random_generator.pareto
-permutation = _random_generator.permutation
-poisson = _random_generator.poisson
-power = _random_generator.power
-random = _random_generator.random
-rayleigh = _random_generator.rayleigh
-shuffle = _random_generator.shuffle
-standard_cauchy = _random_generator.standard_cauchy
-standard_exponential = _random_generator.standard_exponential
-standard_gamma = _random_generator.standard_gamma
-standard_normal = _random_generator.standard_normal
-standard_t = _random_generator.standard_t
-triangular = _random_generator.triangular
-uniform = _random_generator.uniform
-vonmises = _random_generator.vonmises
-wald = _random_generator.wald
-weibull = _random_generator.weibull
-zipf = _random_generator.zipf
 
+def default_gen(seed=None):
+    """Construct a new Generator with the default BitGenerator (PCG64).
+
+    Parameters
+    ----------
+    seed : {None, int, array_like[ints], ISeedSequence, BitGenerator, Generator}, optional
+        A seed to initialize the `BitGenerator`. If None, then fresh,
+        unpredictable entropy will be pulled from the OS. If an ``int`` or
+        ``array_like[ints]`` is passed, then it will be passed to
+        `SeedSequence` to derive the initial `BitGenerator` state. One may also
+        pass in an implementor of the `ISeedSequence` interface like
+        `SeedSequence`.
+        Additionally, when passed a `BitGenerator`, it will be wrapped by
+        `Generator`. If passed a `Generator`, it will be returned unaltered.
+
+    Notes
+    -----
+    When `seed` is omitted or ``None``, a new `BitGenerator` and `Generator` will
+    be instantiated each time. This function does not manage a default global
+    instance.
+    """
+    if _check_bit_generator(seed):
+        # We were passed a BitGenerator, so just wrap it up.
+        return Generator(seed)
+    elif isinstance(seed, Generator):
+        # Pass through a Generator.
+        return seed
+    # Otherwise we need to instantiate a new BitGenerator and Generator as
+    # normal.
+    return Generator(PCG64(seed))
diff --git a/numpy/random/mt19937.pyx b/numpy/random/mt19937.pyx
index db571cbef..49c3622f5 100644
--- a/numpy/random/mt19937.pyx
+++ b/numpy/random/mt19937.pyx
@@ -43,16 +43,19 @@ cdef uint64_t mt19937_raw(void *st) nogil:
 
 cdef class MT19937(BitGenerator):
     """
-    MT19937(seed_seq=None)
+    MT19937(seed=None)
 
     Container for the Mersenne Twister pseudo-random number generator.
 
     Parameters
     ----------
-    seed_seq : {None, SeedSequence, int, array_like[ints]}, optional
-        A SeedSequence to initialize the BitGenerator. If None, one will be
-        created. If an int or array_like[ints], it will be used as the entropy
-        for creating a SeedSequence.
+    seed : {None, int, array_like[ints], ISeedSequence}, optional
+        A seed to initialize the `BitGenerator`. If None, then fresh,
+        unpredictable entropy will be pulled from the OS. If an ``int`` or
+        ``array_like[ints]`` is passed, then it will be passed to
+        `SeedSequence` to derive the initial `BitGenerator` state. One may also
+        pass in an implementor of the `ISeedSequence` interface like
+        `SeedSequence`.
 
     Attributes
     ----------
diff --git a/numpy/random/pcg64.pyx b/numpy/random/pcg64.pyx
index 7cfa19771..c07ec09cb 100644
--- a/numpy/random/pcg64.pyx
+++ b/numpy/random/pcg64.pyx
@@ -43,10 +43,13 @@ cdef class PCG64(BitGenerator):
 
     Parameters
     ----------
-    seed_seq : {None, SeedSequence, int, array_like[ints]}, optional
-        A SeedSequence to initialize the BitGenerator. If None, one will be
-        created. If an int or array_like[ints], it will be used as the entropy
-        for creating a SeedSequence.
+    seed : {None, int, array_like[ints], ISeedSequence}, optional
+        A seed to initialize the `BitGenerator`. If None, then fresh,
+        unpredictable entropy will be pulled from the OS. If an ``int`` or
+        ``array_like[ints]`` is passed, then it will be passed to
+        `SeedSequence` to derive the initial `BitGenerator` state. One may also
+        pass in an implementor of the `ISeedSequence` interface like
+        `SeedSequence`.
 
     Notes
     -----
@@ -99,8 +102,8 @@ cdef class PCG64(BitGenerator):
     cdef pcg64_state rng_state
     cdef pcg64_random_t pcg64_random_state
 
-    def __init__(self, seed_seq=None):
-        BitGenerator.__init__(self, seed_seq)
+    def __init__(self, seed=None):
+        BitGenerator.__init__(self, seed)
         self.rng_state.pcg_state = &self.pcg64_random_state
 
         self._bitgen.state = <void *>&self.rng_state
diff --git a/numpy/random/philox.pyx b/numpy/random/philox.pyx
index 1594242d8..8b7683017 100644
--- a/numpy/random/philox.pyx
+++ b/numpy/random/philox.pyx
@@ -62,18 +62,21 @@ cdef class Philox(BitGenerator):
 
     Parameters
     ----------
-    seed_seq : {None, SeedSequence, int, array_like[ints]}, optional
-        A SeedSequence to initialize the BitGenerator. If None, one will be
-        created. If an int or array_like[ints], it will be used as the entropy
-        for creating a SeedSequence.
+    seed : {None, int, array_like[ints], ISeedSequence}, optional
+        A seed to initialize the `BitGenerator`. If None, then fresh,
+        unpredictable entropy will be pulled from the OS. If an ``int`` or
+        ``array_like[ints]`` is passed, then it will be passed to
+        `SeedSequence` to derive the initial `BitGenerator` state. One may also
+        pass in an implementor of the `ISeedSequence` interface like
+        `SeedSequence`.
     counter : {None, int, array_like}, optional
         Counter to use in the Philox state. Can be either
         a Python int (long in 2.x) in [0, 2**256) or a 4-element uint64 array.
         If not provided, the RNG is initialized at 0.
     key : {None, int, array_like}, optional
         Key to use in the Philox state.  Unlike seed, the value in key is
-        directly set. Can be either a Python int (long in 2.x) in [0, 2**128)
-        or a 2-element uint64 array. `key` and `seed` cannot both be used.
+        directly set. Can be either a Python int in [0, 2**128) or a 2-element
+        uint64 array. `key` and `seed` cannot both be used.
 
     Attributes
     ----------
@@ -143,7 +146,7 @@ cdef class Philox(BitGenerator):
 
     **Compatibility Guarantee**
 
-    ``Philox`` makes a guarantee that a fixed seed and will always produce
+    ``Philox`` makes a guarantee that a fixed seed will always produce
     the same random integer stream.
 
     Examples
@@ -164,16 +167,18 @@ cdef class Philox(BitGenerator):
     cdef philox4x64_key_t philox_key
     cdef philox4x64_ctr_t philox_ctr
 
-    def __init__(self, seed_seq=None, counter=None, key=None):
-        if seed_seq is not None and key is not None:
+    def __init__(self, seed=None, counter=None, key=None):
+        if seed is not None and key is not None:
             raise ValueError('seed and key cannot be both used')
-        BitGenerator.__init__(self, seed_seq)
+        BitGenerator.__init__(self, seed)
         self.rng_state.ctr = &self.philox_ctr
         self.rng_state.key = &self.philox_key
         if key is not None:
             key = int_to_array(key, 'key', 128, 64)
             for i in range(2):
                 self.rng_state.key.v[i] = key[i]
+            # The seed sequence is invalid.
+            self._seed_seq = None
         else:
             key = self._seed_seq.generate_state(2, np.uint64)
             for i in range(2):
diff --git a/numpy/random/sfc64.pyx b/numpy/random/sfc64.pyx
index 74a07da73..1afe3bd0d 100644
--- a/numpy/random/sfc64.pyx
+++ b/numpy/random/sfc64.pyx
@@ -32,16 +32,19 @@ cdef double sfc64_double(void* st) nogil:
 
 cdef class SFC64(BitGenerator):
     """
-    SFC64(seed_seq=None)
+    SFC64(seed=None)
 
     BitGenerator for Chris Doty-Humphrey's Small Fast Chaotic PRNG.
 
     Parameters
     ----------
-    seed_seq : {None, ISeedSequence, int, array_like[ints]}, optional
-        A SeedSequence to initialize the BitGenerator. If None, one will be
-        created. If an int or array_like[ints], it will be used as the entropy
-        for creating a SeedSequence.
+    seed : {None, int, array_like[ints], ISeedSequence}, optional
+        A seed to initialize the `BitGenerator`. If None, then fresh,
+        unpredictable entropy will be pulled from the OS. If an ``int`` or
+        ``array_like[ints]`` is passed, then it will be passed to
+        `SeedSequence` to derive the initial `BitGenerator` state. One may also
+        pass in an implementor of the `ISeedSequence` interface like
+        `SeedSequence`.
 
     Notes
     -----
@@ -81,8 +84,8 @@ cdef class SFC64(BitGenerator):
 
     cdef sfc64_state rng_state
 
-    def __init__(self, seed_seq=None):
-        BitGenerator.__init__(self, seed_seq)
+    def __init__(self, seed=None):
+        BitGenerator.__init__(self, seed)
         self._bitgen.state = <void *>&self.rng_state
         self._bitgen.next_uint64 = &sfc64_uint64
         self._bitgen.next_uint32 = &sfc64_uint32
diff --git a/numpy/random/tests/test_direct.py b/numpy/random/tests/test_direct.py
index 1e824fbcc..97cec6e26 100644
--- a/numpy/random/tests/test_direct.py
+++ b/numpy/random/tests/test_direct.py
@@ -7,7 +7,8 @@ from numpy.testing import (assert_equal, assert_allclose, assert_array_equal,
 import pytest
 
 from numpy.random import (
-  Generator, MT19937, PCG64, Philox, RandomState, SeedSequence, SFC64
+    Generator, MT19937, PCG64, Philox, RandomState, SeedSequence, SFC64,
+    default_gen
 )
 from numpy.random.common import interface
 
@@ -400,3 +401,18 @@ class TestSFC64(Base):
         cls.seed_error_type = (ValueError, TypeError)
         cls.invalid_init_types = [(3.2,), ([None],), (1, None)]
         cls.invalid_init_values = [(-1,)]
+
+
+class TestDefaultGen(object):
+    def test_seed(self):
+        for args in [(), (None,), (1234,), ([1234, 5678],)]:
+            rg = default_gen(*args)
+            assert isinstance(rg.bit_generator, PCG64)
+
+    def test_passthrough(self):
+        bg = Philox()
+        rg = default_gen(bg)
+        assert rg.bit_generator is bg
+        rg2 = default_gen(rg)
+        assert rg2 is rg
+        assert rg2.bit_generator is bg
diff --git a/numpy/random/tests/test_smoke.py b/numpy/random/tests/test_smoke.py
index 4263335f6..12feec098 100644
--- a/numpy/random/tests/test_smoke.py
+++ b/numpy/random/tests/test_smoke.py
@@ -764,6 +764,50 @@ class TestSFC64(RNG):
         cls._extra_setup()
 
 
+class TestPCG64(RNG):
+    @classmethod
+    def setup_class(cls):
+        cls.bit_generator = PCG64
+        cls.advance = 2**63 + 2**31 + 2**15 + 1
+        cls.seed = [12345]
+        cls.rg = Generator(cls.bit_generator(*cls.seed))
+        cls.initial_state = cls.rg.bit_generator.state
+        cls.seed_vector_bits = 64
+        cls._extra_setup()
+
+
+class TestDefaultGen(RNG):
+    @classmethod
+    def setup_class(cls):
+        # This will duplicate some tests that directly instantiate a fresh
+        # Generator(), but that's okay.
+        cls.bit_generator = PCG64
+        cls.advance = 2**63 + 2**31 + 2**15 + 1
+        cls.seed = [12345]
+        cls.rg = np.random.default_gen(*cls.seed)
+        cls.initial_state = cls.rg.bit_generator.state
+        cls.seed_vector_bits = 64
+        cls._extra_setup()
+
+    def test_default_is_pcg64(self):
+        # In order to change the default BitGenerator, we'll go through
+        # a deprecation cycle to move to a different function.
+        assert_(isinstance(self.rg.bit_generator, PCG64))
+
+    def test_seed(self):
+        np.random.default_gen()
+        np.random.default_gen(None)
+        np.random.default_gen(12345)
+        np.random.default_gen(0)
+        np.random.default_gen(43660444402423911716352051725018508569)
+        np.random.default_gen([43660444402423911716352051725018508569,
+                               279705150948142787361475340226491943209])
+        with pytest.raises(ValueError):
+            np.random.default_gen(-1)
+        with pytest.raises(ValueError):
+            np.random.default_gen([12345, -1])
+
+
 class TestEntropy(object):
     def test_entropy(self):
         e1 = entropy.random_entropy()