: tuple of ints, default: ()
    The shape of the sample to draw. If the parameters of the distribution
    underlying the random variable are arrays of shape ``param_shape``,
    the output array will be of shape ``shape + param_shape``.
method : {None, 'formula', 'inverse_transform'}
    The strategy used to produce the sample. By default (``None``),
    the infrastructure chooses between the following options,
    listed in order of precedence.

    - ``'formula'``: an implementation specific to the distribution
    - ``'inverse_transform'``: generate a uniformly distributed sample and
      return the inverse CDF at these arguments.

    Not all `method` options are available for all distributions.
    If the selected `method` is not available, a `NotImplementedError``
    will be raised.
rng : `numpy.random.Generator` or `scipy.stats.QMCEngine`, optional
    Pseudo- or quasi-random number generator state. When `rng` is None,
    a new `numpy.random.Generator` is created using entropy from the
    operating system. Types other than `numpy.random.Generator` and
    `scipy.stats.QMCEngine` are passed to `numpy.random.default_rng`
    to instantiate a ``Generator``.

    If `rng` is an instance of `scipy.stats.QMCEngine` configured to use
    scrambling and `shape` is not empty, then each slice along the zeroth
    axis of the result is a "quasi-independent", low-discrepancy sequence;
    that is, they are distinct sequences that can be treated as statistically
    independent for most practical purposes. Separate calls to `sample`
    produce new quasi-independent, low-discrepancy sequences.

References
----------
.. [1] Sampling (statistics), *Wikipedia*,
       https://en.wikipedia.org/wiki/Sampling_(statistics)

Examples
--------
Instantiate a distribution with the desired parameters:

>>> import numpy as np
>>> from scipy import stats
>>> X = stats.Uniform(a=0., b=1.)

Generate a pseudorandom sample:

>>> x = X.sample((1000, 1))
>>> octiles = (np.arange(8) + 1) / 8
>>> np.count_nonzero(x <= octiles, axis=0)
array([ 148,  263,  387,  516,  636,  751,  865, 1000])  # may vary

>>> X = stats.Uniform(a=np.zeros((3, 1)), b=np.ones(2))
>>> X.a.shape,
(3, 2)
>>> x = X.sample(shape=(5, 4))
>>> x.shape
(5, 4, 3, 2)

r