Repository: sympy/sympy

Issue:
Bspline error with symbolic knots
The following code uses a symbolic value in the knots, and results in an error
```
from sympy import *
x = Symbol('x')

Delta = Symbol('Delta',positive=True)
knots = [i*Delta for i in range(4)]

basis = bspline_basis_set(1, knots, x)
print(basis)
```
The error is
```
Traceback (most recent call last):
  File "bspline_sym.py", line 11, in <module>
    basis = bspline_basis_set(1, knots, x)
  File "/home/sympy/git2/sympy/sympy/functions/special/bsplines.py", line 210, in bspline_basis_set
    return [bspline_basis(d, tuple(knots), i, x) for i in range(n_splines)]
  File "/home/sympy/git2/sympy/sympy/functions/special/bsplines.py", line 210, in <listcomp>
    return [bspline_basis(d, tuple(knots), i, x) for i in range(n_splines)]
  File "/home/sympy/git2/sympy/sympy/functions/special/bsplines.py", line 166, in bspline_basis
    result = _add_splines(A, b1, B, b2)
  File "/home/sympy/git2/sympy/sympy/functions/special/bsplines.py", line 50, in _add_splines
    elif lower_2 < lower and upper_2 <= lower:
  File "/home/sympy/git2/sympy/sympy/core/relational.py", line 384, in __nonzero__
    raise TypeError("cannot determine truth value of Relational")
TypeError: cannot determine truth value of Relational
```

The code in bsplines.py decomposes the conditionals in the Piecewise expression, and expects each side of the conditional to be in particular order - the rhs is expected to the numerical (or fixed) value.
This order is respected in cases where the knots are numbers, or where the knots are expressions.   When a knot is an individual symbol the order changes and break expectations.

The issue appears because the bare symbol `Delta` is a knot.   If the range is changed to something like `range(2,5)`, then the issue does not appear because every knot is a multiplicative expression.


A workaround is to multiply the knots by a dummy symbol and remove the symbol later.
```
dummy = Symbol('dummy',positive=True)
knots = [i*Delta*dummy for i in range(4)]
print(knots)
basis = bspline_basis_set(1, knots, x)
basis = [b.subs(dummy, 1) for b in basis]
print(basis)
```

Discussion:
sympy-bot: :white_check_mark:

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot) (v158). I'm here to help you write a release notes entry. Please read the [guide on how to write release notes](https://github.com/sympy/sympy/wiki/Writing-Release-Notes).



Your release notes are in good order.

Here is what the release notes will look like:
* functions
  * fixed bug which arose when using symbolic knots ([#19272](https://github.com/sympy/sympy/pull/19272) by [@smichr](https://github.com/smichr))

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.7.

Note: This comment will be updated with the latest check if you edit the pull request. You need to reload the page to see it. <details><summary>Click here to see the pull request description that was parsed.</summary>

    <!-- Your title above should be a short description of what
    was changed. Do not include the issue number in the title. -->

    #### References to other Issues or PRs
    <!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
    format, e.g. "Fixes #1234" (see
    https://tinyurl.com/auto-closing for more information). Also, please
    write a comment on that issue linking back to this pull request once it is
    open. -->

    fixes #19262

    #### Brief description of what is fixed or changed


    #### Other comments


    #### Release Notes

    <!-- Write the release notes for this release below. See
    https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more information
    on how to write release notes. The bot will check your release notes
    automatically to see if they are formatted correctly. -->

    <!-- BEGIN RELEASE NOTES -->
    * functions
        * fixed bug which arose when using symbolic knots
    <!-- END RELEASE NOTES -->

</details><p>

codecov[bot]: # [Codecov](https://codecov.io/gh/sympy/sympy/pull/19272?src=pr&el=h1) Report
> Merging [#19272](https://codecov.io/gh/sympy/sympy/pull/19272?src=pr&el=desc) into [master](https://codecov.io/gh/sympy/sympy/commit/b561fee8a4810924f9481a5d78a83436b544a672&el=desc) will **decrease** coverage by `0.035%`.
> The diff coverage is `84.615%`.

```diff
@@              Coverage Diff              @@
##            master    #19272       +/-   ##
=============================================
- Coverage   75.649%   75.613%   -0.036%     
=============================================
  Files          651       651               
  Lines       169487    169466       -21     
  Branches     40021     40015        -6     
=============================================
- Hits        128216    128140       -76     
- Misses       35654     35721       +67     
+ Partials      5617      5605       -12     
```

jksuom: How about using a dummy with no assumptions when constructing the spline polynomials and replacing that with x only in the output?

Patch:
diff --git a/.travis.yml b/.travis.yml
index 753a1ecb79..346eaaadfd 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -53,7 +53,7 @@ matrix:
       # all of the dependencies are supported on 3.8.
       env:
         - TEST_ASCII="true"
-        - TEST_OPT_DEPENDENCY="matchpy numpy scipy gmpy2 matplotlib theano llvmlite autowrap cython wurlitzer python-symengine=0.5.1 tensorflow numexpr ipython antlr-python-runtime>=4.7,<4.8 antlr>=4.7,<4.8 cloudpickle pyglet pycosat lfortran python-clang lxml"
+        - TEST_OPT_DEPENDENCY="matchpy numpy scipy gmpy2 matplotlib<3.2 theano llvmlite autowrap cython wurlitzer python-symengine=0.5.1 tensorflow numexpr ipython antlr-python-runtime>=4.7,<4.8 antlr>=4.7,<4.8 cloudpickle pyglet pycosat lfortran python-clang lxml"
         - TEST_SAGE="true"
         - SYMPY_STRICT_COMPILER_CHECKS=1
       addons:
@@ -239,7 +239,7 @@ before_install:
         pip list | grep "numpy" && pip uninstall -y numpy;
     fi
   - if [[ "${TEST_SPHINX}" == "true" ]]; then
-      pip install "sphinx" "docutils" doctr matplotlib sphinx-math-dollar;
+      pip install "sphinx" "docutils" doctr "matplotlib<3.2" sphinx-math-dollar;
     fi
   - if [[ "${BENCHMARK}" == "true" ]]; then
       pip install asv;
diff --git a/doc/src/modules/stats.rst b/doc/src/modules/stats.rst
index ac03008709..05874ec7ec 100644
--- a/doc/src/modules/stats.rst
+++ b/doc/src/modules/stats.rst
@@ -117,6 +117,8 @@ Interface
 .. autofunction:: std
 .. autofunction:: sample
 .. autofunction:: sample_iter
+.. autofunction:: sympy.stats.rv.sample_iter_lambdify
+.. autofunction:: sympy.stats.rv.sample_iter_subs
 .. autofunction:: sympy.stats.rv.sampling_density
 .. autofunction:: sympy.stats.rv.sampling_P
 .. autofunction:: sympy.stats.rv.sampling_E
diff --git a/sympy/functions/special/bsplines.py b/sympy/functions/special/bsplines.py
index 97d6c1e270..24d7f7fd01 100644
--- a/sympy/functions/special/bsplines.py
+++ b/sympy/functions/special/bsplines.py
@@ -5,8 +5,25 @@
 from sympy.core.cache import lru_cache
 
 
-def _add_splines(c, b1, d, b2):
+def _ivl(cond, x):
+    """return the interval corresponding to the condition
+
+    Conditions in spline's Piecewise give the range over
+    which an expression is valid like (lo <= x) & (x <= hi).
+    This function returns (lo, hi).
+    """
+    from sympy.logic.boolalg import And
+    if isinstance(cond, And) and len(cond.args) == 2:
+        a, b = cond.args
+        if a.lts == x:
+            a, b = b, a
+        return a.lts, b.gts
+    raise TypeError('unexpected cond type: %s' % cond)
+
+
+def _add_splines(c, b1, d, b2, x):
     """Construct c*b1 + d*b2."""
+
     if b1 == S.Zero or c == S.Zero:
         rv = piecewise_fold(d * b2)
     elif b2 == S.Zero or d == S.Zero:
@@ -23,23 +40,17 @@ def _add_splines(c, b1, d, b2):
         # This merging algorithm assumes the conditions in
         # p1 and p2 are sorted
         for arg in p1.args[:-1]:
-            # Conditional of Piecewise are And objects
-            # the args of the And object is a tuple of two
-            # Relational objects the numerical value is in the .rhs
-            # of the Relational object
             expr = arg.expr
             cond = arg.cond
 
-            lower = cond.args[0].rhs
+            lower = _ivl(cond, x)[0]
 
             # Check p2 for matching conditions that can be merged
             for i, arg2 in enumerate(p2args):
                 expr2 = arg2.expr
                 cond2 = arg2.cond
 
-                lower_2 = cond2.args[0].rhs
-                upper_2 = cond2.args[1].rhs
-
+                lower_2, upper_2 = _ivl(cond2, x)
                 if cond2 == cond:
                     # Conditions match, join expressions
                     expr += expr2
@@ -64,7 +75,7 @@ def _add_splines(c, b1, d, b2):
         # Add final (0, True)
         new_args.append((0, True))
 
-        rv = Piecewise(*new_args)
+        rv = Piecewise(*new_args, evaluate=False)
 
     return rv.expand()
 
@@ -137,6 +148,11 @@ def bspline_basis(d, knots, n, x):
     .. [1] https://en.wikipedia.org/wiki/B-spline
 
     """
+    from sympy.core.symbol import Dummy
+    # make sure x has no assumptions so conditions don't evaluate
+    xvar = x
+    x = Dummy()
+
     knots = tuple(sympify(k) for k in knots)
     d = int(d)
     n = int(n)
@@ -163,10 +179,12 @@ def bspline_basis(d, knots, n, x):
         else:
             b1 = A = S.Zero
 
-        result = _add_splines(A, b1, B, b2)
+        result = _add_splines(A, b1, B, b2, x)
     else:
         raise ValueError("degree must be non-negative: %r" % n)
-    return result
+
+    # return result with user-given x
+    return result.xreplace({x: xvar})
 
 
 def bspline_basis_set(d, knots, x):
@@ -242,7 +260,7 @@ def interpolating_spline(d, x, X, Y):
     bspline_basis_set, interpolating_poly
 
     """
-    from sympy import symbols, Number, Dummy, Rational
+    from sympy import symbols, Dummy
     from sympy.solvers.solveset import linsolve
     from sympy.matrices.dense import Matrix
 
@@ -256,6 +274,7 @@ def interpolating_spline(d, x, X, Y):
         raise ValueError("Degree must be less than the number of control points.")
     if not all(a < b for a, b in zip(X, X[1:])):
         raise ValueError("The x-coordinates must be strictly increasing.")
+    X = [sympify(i) for i in X]
 
     # Evaluating knots value
     if d.is_odd:
@@ -264,7 +283,7 @@ def interpolating_spline(d, x, X, Y):
     else:
         j = d // 2
         interior_knots = [
-            Rational(a + b, 2) for a, b in zip(X[j : -j - 1], X[j + 1 : -j])
+            (a + b)/2 for a, b in zip(X[j : -j - 1], X[j + 1 : -j])
         ]
 
     knots = [X[0]] * (d + 1) + list(interior_knots) + [X[-1]] * (d + 1)
@@ -279,8 +298,7 @@ def interpolating_spline(d, x, X, Y):
 
     # Sorting the intervals
     #  ival contains the end-points of each interval
-    ival = [e.atoms(Number) for e in intervals]
-    ival = [list(sorted(e))[0] for e in ival]
+    ival = [_ivl(c, x) for c in intervals]
     com = zip(ival, intervals)
     com = sorted(com, key=lambda x: x[0])
     intervals = [y for x, y in com]
diff --git a/sympy/polys/factortools.py b/sympy/polys/factortools.py
index cfa59ae9b8..93298f47e6 100644
--- a/sympy/polys/factortools.py
+++ b/sympy/polys/factortools.py
@@ -124,64 +124,13 @@ def dmp_trial_division(f, factors, u, K):
 
 
 def dup_zz_mignotte_bound(f, K):
-    """
-    The Knuth-Cohen variant of Mignotte bound for
-    univariate polynomials in `K[x]`.
-
-    Examples
-    ========
-
-    >>> from sympy.polys import ring, ZZ
-    >>> R, x = ring("x", ZZ)
-
-    >>> f = x**3 + 14*x**2 + 56*x + 64
-    >>> R.dup_zz_mignotte_bound(f)
-    152
-
-    By checking `factor(f)` we can see that max coeff is 8
-
-    Also consider a case that `f` is irreducible for example `f = 2*x**2 + 3*x + 4`
-    To avoid a bug for these cases, we return the bound plus the max coefficient of `f`
-
-    >>> f = 2*x**2 + 3*x + 4
-    >>> R.dup_zz_mignotte_bound(f)
-    6
-
-    Lastly,To see the difference between the new and the old Mignotte bound
-    consider the irreducible polynomial::
-
-    >>> f = 87*x**7 + 4*x**6 + 80*x**5 + 17*x**4 + 9*x**3 + 12*x**2 + 49*x + 26
-    >>> R.dup_zz_mignotte_bound(f)
-    744
-
-    The new Mignotte bound is 744 whereas the old one (SymPy 1.5.1) is 1937664.
-
-
-    References
-    ==========
-
-    ..[1] [Abbott2013]_
-
-    """
-    from sympy import binomial
-
-    d = dup_degree(f)
-    delta = _ceil(d / 2)
-    delta2 = _ceil(delta / 2)
-
-    # euclidean-norm
-    eucl_norm = K.sqrt( sum( [cf**2 for cf in f] ) )
-
-    # biggest values of binomial coefficients (p. 538 of reference)
-    t1 = binomial(delta - 1, delta2)
-    t2 = binomial(delta - 1, delta2 - 1)
+    """Mignotte bound for univariate polynomials in `K[x]`. """
+    a = dup_max_norm(f, K)
+    b = abs(dup_LC(f, K))
+    n = dup_degree(f)
 
-    lc = K.abs(dup_LC(f, K))   # leading coefficient
-    bound = t1 * eucl_norm + t2 * lc   # (p. 538 of reference)
-    bound += dup_max_norm(f, K) # add max coeff for irreducible polys
-    bound = _ceil(bound / 2) * 2   # round up to even integer
+    return K.sqrt(K(n + 1))*2**n*a*b
 
-    return bound
 
 def dmp_zz_mignotte_bound(f, u, K):
     """Mignotte bound for multivariate polynomials in `K[X]`. """
diff --git a/sympy/polys/polytools.py b/sympy/polys/polytools.py
index f066202e69..6824c7f410 100644
--- a/sympy/polys/polytools.py
+++ b/sympy/polys/polytools.py
@@ -132,7 +132,7 @@ class Poly(Basic):
 
     """
 
-    __slots__ = ('rep', 'gens')
+    __slots__ = ('rep',)
 
     is_commutative = True
     is_Poly = True
@@ -158,15 +158,6 @@ def __new__(cls, rep, *gens, **args):
             else:
                 return cls._from_expr(rep, opt)
 
-    # Poly does not pass its args to Basic.__new__ to be stored in _args so we
-    # have to emulate them here with an args property that derives from rep
-    # and gens which are instance attributes. This also means we need to
-    # define _hashable_content. The _hashable_content is rep and gens but args
-    # uses expr instead of rep (expr is the Basic version of rep). Passing
-    # expr in args means that Basic methods like subs should work. Using rep
-    # otherwise means that Poly can remain more efficient than Basic by
-    # avoiding creating a Basic instance just to be hashable.
-
     @classmethod
     def new(cls, rep, *gens):
         """Construct :class:`Poly` instance from raw representation. """
@@ -176,22 +167,20 @@ def new(cls, rep, *gens):
         elif rep.lev != len(gens) - 1:
             raise PolynomialError("invalid arguments: %s, %s" % (rep, gens))
 
-        obj = Basic.__new__(cls)
+        expr = basic_from_dict(rep.to_sympy_dict(), *gens)
+
+        obj = Basic.__new__(cls, expr, *gens)
         obj.rep = rep
-        obj.gens = gens
 
         return obj
 
     @property
     def expr(self):
-        return basic_from_dict(self.rep.to_sympy_dict(), *self.gens)
+        return self.args[0]
 
     @property
-    def args(self):
-        return (self.expr,) + self.gens
-
-    def _hashable_content(self):
-        return (self.rep,) + self.gens
+    def gens(self):
+        return self.args[1:]
 
     @classmethod
     def from_dict(cls, rep, *gens, **args):
diff --git a/sympy/stats/crv.py b/sympy/stats/crv.py
index 61a8fb7a19..f684a4f21f 100644
--- a/sympy/stats/crv.py
+++ b/sympy/stats/crv.py
@@ -12,21 +12,17 @@
 
 from sympy import (Interval, Intersection, symbols, sympify, Dummy, nan,
         Integral, And, Or, Piecewise, cacheit, integrate, oo, Lambda,
-        Basic, S, exp, I, FiniteSet, Ne, Eq, Union, poly, series, factorial,
-        lambdify)
+        Basic, S, exp, I, FiniteSet, Ne, Eq, Union, poly, series, factorial)
 from sympy.core.function import PoleError
 from sympy.functions.special.delta_functions import DiracDelta
 from sympy.polys.polyerrors import PolynomialError
 from sympy.solvers.solveset import solveset
 from sympy.solvers.inequalities import reduce_rational_inequalities
 from sympy.core.sympify import _sympify
-from sympy.external import import_module
 from sympy.stats.rv import (RandomDomain, SingleDomain, ConditionalDomain,
         ProductDomain, PSpace, SinglePSpace, random_symbols, NamedArgsMixin)
+import random
 
-scipy = import_module('scipy')
-numpy = import_module('numpy')
-pymc3 = import_module('pymc3')
 
 class ContinuousDomain(RandomDomain):
     """
@@ -149,132 +145,6 @@ def __call__(self, *args):
         return self.pdf(*args)
 
 
-class SampleExternalContinuous:
-    """Class consisting of the methods that are used to sample values of random
-    variables from external libraries."""
-
-
-    scipy_rv_map = {
-        'BetaDistribution': lambda dist, size: scipy.stats.beta.rvs(a=float(dist.alpha),
-            b=float(dist.beta), size=size),
-        'BetaPrimeDistribution':lambda dist, size: scipy.stats.betaprime.rvs(a=float(dist.alpha),
-            b=float(dist.beta), size=size),
-        'CauchyDistribution': lambda dist, size: scipy.stats.cauchy.rvs(loc=float(dist.x0),
-            scale=float(dist.gamma), size=size),
-        'ChiDistribution': lambda dist, size: scipy.stats.chi.rvs(df=float(dist.k),
-            size=size),
-        'ChiSquaredDistribution': lambda dist, size: scipy.stats.chi2.rvs(df=float(dist.k),
-            size=size),
-        'ExponentialDistribution': lambda dist, size: scipy.stats.expon.rvs(loc=0,
-            scale=1/float(dist.rate), size=size),
-        'GammaDistribution': lambda dist, size: scipy.stats.gamma.rvs(a=float(dist.k), loc=0,
-            scale=float(dist.theta), size=size),
-        'GammaInverseDistribution': lambda dist, size: scipy.stats.invgamma.rvs(a=float(dist.a),
-            loc=0, scale=float(dist.b), size=size),
-        'LogNormalDistribution': lambda dist, size: scipy.stats.lognorm.rvs(s=float(dist.std),
-            loc=0, scale=exp(float(dist.mean)), size=size),
-        'NormalDistribution': lambda dist, size: scipy.stats.norm.rvs(float(dist.mean),
-            float(dist.std), size=size),
-        'GaussianInverseDistribution': lambda dist, size: scipy.stats.invgauss.rvs(
-            mu=float(dist.mean)/float(dist.shape), scale=float(dist.shape), size=size),
-        'ParetoDistribution': lambda dist, size: scipy.stats.pareto.rvs(b=float(dist.alpha),
-            scale=float(dist.xm), size=size),
-        'StudentTDistribution': lambda dist, size: scipy.stats.t.rvs(df=float(dist.nu),
-            size=size),
-        'UniformDistribution': lambda dist, size: scipy.stats.uniform.rvs(loc=float(dist.left),
-            scale=float(dist.right)-float(dist.left), size=size),
-        'WeibullDistribution': lambda dist, size: scipy.stats.weibull_min.rvs(loc=0,
-            c=float(dist.beta), scale=float(dist.alpha), size=size)
-        }
-
-    numpy_rv_map = {
-        'BetaDistribution': lambda dist, size: numpy.random.beta(a=float(dist.alpha),
-            b=float(dist.beta), size=size),
-        'ChiSquaredDistribution': lambda dist, size: numpy.random.chisquare(
-            df=float(dist.k), size=size),
-        'ExponentialDistribution': lambda dist, size: numpy.random.exponential(
-            1/float(dist.rate), size=size),
-        'GammaDistribution': lambda dist, size: numpy.random.gamma(float(dist.k),
-            float(dist.theta), size=size),
-        'LogNormalDistribution': lambda dist, size: numpy.random.lognormal(
-            float(dist.mean), float(dist.std), size=size),
-        'NormalDistribution': lambda dist, size: numpy.random.normal(
-            float(dist.mean), float(dist.std), size=size),
-        'ParetoDistribution': lambda dist, size: (numpy.random.pareto(
-            a=float(dist.alpha), size=size) + 1) * float(dist.xm),
-        'UniformDistribution': lambda dist, size: numpy.random.uniform(
-            low=float(dist.left), high=float(dist.right), size=size)
-    }
-
-    pymc3_rv_map = {
-        'BetaDistribution': lambda dist:
-            pymc3.Beta('X', alpha=float(dist.alpha), beta=float(dist.beta)),
-        'CauchyDistribution': lambda dist:
-            pymc3.Cauchy('X', alpha=float(dist.x0), beta=float(dist.gamma)),
-        'ChiSquaredDistribution': lambda dist:
-            pymc3.ChiSquared('X', nu=float(dist.k)),
-        'ExponentialDistribution': lambda dist:
-            pymc3.Exponential('X', lam=float(dist.rate)),
-        'GammaDistribution': lambda dist:
-            pymc3.Gamma('X', alpha=float(dist.k), beta=1/float(dist.theta)),
-        'LogNormalDistribution': lambda dist:
-            pymc3.Lognormal('X', mu=float(dist.mean), sigma=float(dist.std)),
-        'NormalDistribution': lambda dist:
-            pymc3.Normal('X', float(dist.mean), float(dist.std)),
-        'GaussianInverseDistribution': lambda dist:
-            pymc3.Wald('X', mu=float(dist.mean), lam=float(dist.shape)),
-        'ParetoDistribution': lambda dist:
-            pymc3.Pareto('X', alpha=float(dist.alpha), m=float(dist.xm)),
-        'UniformDistribution': lambda dist:
-            pymc3.Uniform('X', lower=float(dist.left), upper=float(dist.right))
-    }
-
-    @classmethod
-    def _sample_scipy(cls, dist, size):
-        """Sample from SciPy."""
-
-        dist_list = cls.scipy_rv_map.keys()
-
-        if dist.__class__.__name__ == 'ContinuousDistributionHandmade':
-            from scipy.stats import rv_continuous
-            z = Dummy('z')
-            handmade_pdf = lambdify(z, dist.pdf(z), 'scipy')
-            class scipy_pdf(rv_continuous):
-                def _pdf(self, x):
-                    return handmade_pdf(x)
-            scipy_rv = scipy_pdf(a=dist.set._inf, b=dist.set._sup, name='scipy_pdf')
-            return scipy_rv.rvs(size=size)
-
-        if dist.__class__.__name__ not in dist_list:
-            return None
-
-        return cls.scipy_rv_map[dist.__class__.__name__](dist, size)
-
-    @classmethod
-    def _sample_numpy(cls, dist, size):
-        """Sample from NumPy."""
-
-        dist_list = cls.numpy_rv_map.keys()
-
-        if dist.__class__.__name__ not in dist_list:
-            return None
-
-        return cls.numpy_rv_map[dist.__class__.__name__](dist, size)
-
-    @classmethod
-    def _sample_pymc3(cls, dist, size):
-        """Sample from PyMC3."""
-
-        dist_list = cls.pymc3_rv_map.keys()
-
-        if dist.__class__.__name__ not in dist_list:
-            return None
-
-        with pymc3.Model():
-            cls.pymc3_rv_map[dist.__class__.__name__](dist)
-            return pymc3.sample(size, chains=1, progressbar=False)[:]['X']
-
-
 class SingleContinuousDistribution(ContinuousDistribution, NamedArgsMixin):
     """ Continuous distribution of a single variable
 
@@ -301,25 +171,34 @@ def __new__(cls, *args):
     def check(*args):
         pass
 
-    def sample(self, size=1, library='scipy'):
+    def sample(self, size=()):
         """ A random realization from the distribution """
+        icdf = self._inverse_cdf_expression()
+        if not size:
+            return icdf(random.uniform(0, 1))
+        else:
+            return [icdf(random.uniform(0, 1))]*size
 
-        libraries = ['scipy', 'numpy', 'pymc3']
-        if library not in libraries:
-            raise NotImplementedError("Sampling from %s is not supported yet."
-                                        % str(library))
-        if not import_module(library):
-            raise ValueError("Failed to import %s" % library)
+    @cacheit
+    def _inverse_cdf_expression(self):
+        """ Inverse of the CDF
 
-        samps = getattr(SampleExternalContinuous, '_sample_' + library)(self, size)
+        Used by sample
+        """
+        x, z = symbols('x, z', positive=True, cls=Dummy)
+        # Invert CDF
+        try:
+            inverse_cdf = solveset(self.cdf(x) - z, x, S.Reals)
+            if isinstance(inverse_cdf, Intersection) and S.Reals in inverse_cdf.args:
+                inverse_cdf = list(inverse_cdf.args[1])
+        except NotImplementedError:
+            inverse_cdf = None
+        if not inverse_cdf or len(inverse_cdf) != 1:
+            raise NotImplementedError("Could not invert CDF")
 
-        if samps is not None:
-            return samps
-        raise NotImplementedError(
-                "Sampling for %s is not currently implemented from %s"
-                % (self.__class__.__name__, library)
-                )
+        (icdf,) = inverse_cdf
 
+        return Lambda(z, icdf)
 
     @cacheit
     def compute_cdf(self, **kwargs):
@@ -616,13 +495,13 @@ def set(self):
     def domain(self):
         return SingleContinuousDomain(sympify(self.symbol), self.set)
 
-    def sample(self, size=1, library='scipy'):
+    def sample(self, size=()):
         """
         Internal sample method
 
         Returns dictionary mapping RandomSymbol to realization value.
         """
-        return {self.value: self.distribution.sample(size, library=library)}
+        return {self.value: self.distribution.sample(size)}
 
     def compute_expectation(self, expr, rvs=None, evaluate=False, **kwargs):
         rvs = rvs or (self.value,)
diff --git a/sympy/stats/crv_types.py b/sympy/stats/crv_types.py
index 576ad128b9..a12502e2ac 100644
--- a/sympy/stats/crv_types.py
+++ b/sympy/stats/crv_types.py
@@ -54,6 +54,7 @@
 
 from __future__ import print_function, division
 
+import random
 
 from sympy import beta as beta_fn
 from sympy import cos, sin, tan, atan, exp, besseli, besselj, besselk
@@ -61,6 +62,7 @@
                    Piecewise, And, Eq, binomial, factorial, Sum, floor, Abs,
                    Lambda, Basic, lowergamma, erf, erfc, erfi, erfinv, I, asin,
                    hyper, uppergamma, sinh, Ne, expint, Rational, integrate)
+from sympy.external import import_module
 from sympy.matrices import MatrixBase, MatrixExpr
 from sympy.stats.crv import SingleContinuousPSpace, SingleContinuousDistribution
 from sympy.stats.joint_rv import JointPSpace, CompoundDistribution
@@ -377,6 +379,12 @@ def pdf(self, x):
         alpha, beta = self.alpha, self.beta
         return x**(alpha - 1) * (1 - x)**(beta - 1) / beta_fn(alpha, beta)
 
+    def sample(self, size=()):
+        if not size:
+            return random.betavariate(self.alpha, self.beta)
+        else:
+            return [random.betavariate(self.alpha, self.beta)]*size
+
     def _characteristic_function(self, t):
         return hyper((self.alpha,), (self.alpha + self.beta,), I*t)
 
@@ -1183,6 +1191,12 @@ def check(rate):
     def pdf(self, x):
         return self.rate * exp(-self.rate*x)
 
+    def sample(self, size=()):
+        if not size:
+            return random.expovariate(self.rate)
+        else:
+            return [random.expovariate(self.rate)]*size
+
     def _cdf(self, x):
         return Piecewise(
                 (S.One - exp(-self.rate*x), x >= 0),
@@ -1612,6 +1626,12 @@ def pdf(self, x):
         k, theta = self.k, self.theta
         return x**(k - 1) * exp(-x/theta) / (gamma(k)*theta**k)
 
+    def sample(self, size=()):
+        if not size:
+            return random.gammavariate(self.k, self.theta)
+        else:
+            return [random.gammavariate(self.k, self.theta)]*size
+
     def _cdf(self, x):
         k, theta = self.k, self.theta
         return Piecewise(
@@ -1719,6 +1739,14 @@ def _cdf(self, x):
         return Piecewise((uppergamma(a,b/x)/gamma(a), x > 0),
                         (S.Zero, True))
 
+    def sample(self, size=()):
+        scipy = import_module('scipy')
+        if scipy:
+            from scipy.stats import invgamma
+            return invgamma.rvs(float(self.a), 0, float(self.b), size=size)
+        else:
+            raise NotImplementedError('Sampling the Inverse Gamma Distribution requires Scipy.')
+
     def _characteristic_function(self, t):
         a, b = self.a, self.b
         return 2 * (-I*b*t)**(a/2) * besselk(a, sqrt(-4*I*b*t)) / gamma(a)
@@ -2386,6 +2414,12 @@ def pdf(self, x):
         mean, std = self.mean, self.std
         return exp(-(log(x) - mean)**2 / (2*std**2)) / (x*sqrt(2*pi)*std)
 
+    def sample(self, size=()):
+        if not size:
+            return random.lognormvariate(self.mean, self.std)
+        else:
+            return [random.lognormvariate(self.mean, self.std)]*size
+
     def _cdf(self, x):
         mean, std = self.mean, self.std
         return Piecewise(
@@ -2716,6 +2750,12 @@ def check(mean, std):
     def pdf(self, x):
         return exp(-(x - self.mean)**2 / (2*self.std**2)) / (sqrt(2*pi)*self.std)
 
+    def sample(self, size=()):
+        if not size:
+            return random.normalvariate(self.mean, self.std)
+        else:
+            return [random.normalvariate(self.mean, self.std)]*size
+
     def _cdf(self, x):
         mean, std = self.mean, self.std
         return erf(sqrt(2)*(-mean + x)/(2*std))/2 + S.Half
@@ -2844,6 +2884,15 @@ def pdf(self, x):
         mu, s = self.mean, self.shape
         return exp(-s*(x - mu)**2 / (2*x*mu**2)) * sqrt(s/((2*pi*x**3)))
 
+    def sample(self, size=()):
+        scipy = import_module('scipy')
+        if scipy:
+            from scipy.stats import invgauss
+            return invgauss.rvs(float(self.mean/self.shape), 0, float(self.shape), size=size)
+        else:
+            raise NotImplementedError(
+                'Sampling the Inverse Gaussian Distribution requires Scipy.')
+
     def _cdf(self, x):
         from sympy.stats import cdf
         mu, s = self.mean, self.shape
@@ -2948,6 +2997,12 @@ def pdf(self, x):
         xm, alpha = self.xm, self.alpha
         return alpha * xm**alpha / x**(alpha + 1)
 
+    def sample(self, size=()):
+        if not size:
+            return random.paretovariate(self.alpha)
+        else:
+            return [random.paretovariate(self.alpha)]*size
+
     def _cdf(self, x):
         xm, alpha = self.xm, self.alpha
         return Piecewise(
@@ -3797,6 +3852,12 @@ def expectation(self, expr, var, **kwargs):
                               Min(self.left, self.right): self.left})
         return result
 
+    def sample(self, size=()):
+        if not size:
+            return random.uniform(self.left, self.right)
+        else:
+            return [random.uniform(self.left, self.right)]*size
+
 
 def Uniform(name, left, right):
     r"""
@@ -4049,6 +4110,11 @@ def pdf(self, x):
         alpha, beta = self.alpha, self.beta
         return beta * (x/alpha)**(beta - 1) * exp(-(x/alpha)**beta) / alpha
 
+    def sample(self, size=()):
+        if not size:
+            return random.weibullvariate(self.alpha, self.beta)
+        else:
+            return [random.weibullvariate(self.alpha, self.beta)]*size
 
 def Weibull(name, alpha, beta):
     r"""
diff --git a/sympy/stats/drv.py b/sympy/stats/drv.py
index 97a1ac5a4e..53f60e9812 100644
--- a/sympy/stats/drv.py
+++ b/sympy/stats/drv.py
@@ -47,16 +47,19 @@ def __new__(cls, *args):
     def check(*args):
         pass
 
-    def sample(self, size=(1,), library='scipy'):
+    def sample(self, size=()):
         """ A random realization from the distribution"""
-        if hasattr(self,'_sample_scipy') and import_module('scipy'):
-            return self._sample_scipy(size=size)
+        if getattr(self,'_sample_scipy', None) and import_module('scipy'):
+            return self._sample_scipy(size)
         icdf = self._inverse_cdf_expression()
         samp_list = []
         while True:
             sample_ = floor(list(icdf(random.uniform(0, 1)))[0])
             if sample_ >= self.set.inf:
-                samp_list.append(sample_)
+                if not size:
+                    return sample_
+                else:
+                    samp_list.append(sample_)
             if len(samp_list) == size:
                 return samp_list
 
@@ -315,7 +318,7 @@ def set(self):
     def domain(self):
         return SingleDiscreteDomain(self.symbol, self.set)
 
-    def sample(self, size=(1,), library='scipy'):
+    def sample(self, size=()):
         """
         Internal sample method
 
diff --git a/sympy/stats/frv.py b/sympy/stats/frv.py
index e82fabba33..a71b2b6964 100644
--- a/sympy/stats/frv.py
+++ b/sympy/stats/frv.py
@@ -339,7 +339,7 @@ def conditional_space(self, condition):
                 for key, val in self._density.items() if domain._test(key))
         return FinitePSpace(domain, density)
 
-    def sample(self, size=(1,), library='scipy'):
+    def sample(self, size=()):
         """
         Internal sample method
 
diff --git a/sympy/stats/rv.py b/sympy/stats/rv.py
index 2a427915cd..fadb8c5db0 100644
--- a/sympy/stats/rv.py
+++ b/sympy/stats/rv.py
@@ -25,10 +25,6 @@
 from sympy.logic.boolalg import Boolean
 from sympy.sets.sets import FiniteSet, ProductSet, Intersection
 from sympy.solvers.solveset import solveset
-from sympy.external import import_module
-from sympy.utilities.misc import filldedent
-import warnings
-
 
 x = Symbol('x')
 
@@ -390,10 +386,9 @@ def domain(self):
     def density(self):
         raise NotImplementedError("Density not available for ProductSpaces")
 
-    def sample(self, size=(1,), library='scipy'):
+    def sample(self, size=()):
         return {k: v for space in self.spaces
-            for k, v in space.sample(size=size, library=library).items()}
-
+            for k, v in space.sample().items()}
 
     def probability(self, condition, **kwargs):
         cond_inv = False
@@ -721,8 +716,7 @@ def expectation(expr, condition=None, numsamples=None, evaluate=True, **kwargs):
     if not random_symbols(expr):  # expr isn't random?
         return expr
     if numsamples:  # Computing by monte carlo sampling?
-        evalf = kwargs.get('evalf', True)
-        return sampling_E(expr, condition, numsamples=numsamples, evalf=evalf)
+        return sampling_E(expr, condition, numsamples=numsamples)
 
     if expr.has(RandomIndexedSymbol):
         return pspace(expr).compute_expectation(expr, condition, evaluate, **kwargs)
@@ -1023,62 +1017,69 @@ def where(condition, given_condition=None, **kwargs):
     return pspace(condition).where(condition, **kwargs)
 
 
-def sample(expr, condition=None, size=(1,), library='scipy', numsamples=1,
-                                                                    **kwargs):
+def sample(expr, condition=None, size=(), **kwargs):
     """
     A realization of the random expression
 
-    Parameters
-    ==========
-
-    expr : Expression of random variables
-        Expression from which sample is extracted
-    condition : Expr containing RandomSymbols
-        A conditional expression
-    size : int, tuple
-        Represents size of each sample in numsamples
-    library : str
-        - 'scipy' : Sample using scipy
-        - 'numpy' : Sample using numpy
-        - 'pymc3' : Sample using PyMC3
-
-        Choose any of the available options to sample from as string,
-        by default is 'scipy'
-    numsamples : int
-        Number of samples, each with size as ``size``
-
     Examples
     ========
 
     >>> from sympy.stats import Die, sample, Normal
     >>> X, Y, Z = Die('X', 6), Die('Y', 6), Die('Z', 6)
 
-    >>> die_roll = sample(X + Y + Z) # doctest: +SKIP
+    >>> die_roll = sample(X + Y + Z) # A random realization of three dice
     >>> N = Normal('N', 3, 4)
-    >>> samp = next(sample(N))[0] # doctest: +SKIP
-    >>> samp in N.pspace.domain.set # doctest: +SKIP
-    True
-    >>> samp = next(sample(N, N>0))[0] # doctest: +SKIP
-    >>> samp > 0 # doctest: +SKIP
+    >>> samp = sample(N)
+    >>> samp in N.pspace.domain.set
     True
-    >>> samp_list = next(sample(N, size=4)) # doctest: +SKIP
-    >>> [sam in N.pspace.domain.set for sam in samp_list] # doctest: +SKIP
+    >>> samp_list = sample(N, size=4)
+    >>> [sam in N.pspace.domain.set for sam in samp_list]
     [True, True, True, True]
 
-    Returns
-    =======
+    """
+    return next(sample_iter(expr, condition, size=size, numsamples=1))
 
-    sample: iterator object
-        iterator object containing the sample/samples of given expr
 
+def sample_iter(expr, condition=None, size=(), numsamples=S.Infinity, **kwargs):
     """
-    message = ("The return type of sample has been changed to return an "
-                  "iterator object since version 1.7. For more information see "
-                  "https://github.com/sympy/sympy/issues/19061")
-    warnings.warn(filldedent(message))
-    return sample_iter(expr, condition, size=size, library=library,
-                                                        numsamples=numsamples)
+    Returns an iterator of realizations from the expression given a condition
+
+    Parameters
+    ==========
+
+    expr: Expr
+        Random expression to be realized
+    condition: Expr, optional
+        A conditional expression
+    numsamples: integer, optional
+        Length of the iterator (defaults to infinity)
+
+    Examples
+    ========
+
+    >>> from sympy.stats import Normal, sample_iter
+    >>> X = Normal('X', 0, 1)
+    >>> expr = X*X + 3
+    >>> iterator = sample_iter(expr, numsamples=3)
+    >>> list(iterator) # doctest: +SKIP
+    [12, 4, 7]
+
+    See Also
+    ========
+
+    sample
+    sampling_P
+    sampling_E
+    sample_iter_lambdify
+    sample_iter_subs
 
+    """
+    # lambdify is much faster but not as robust
+    try:
+        return sample_iter_lambdify(expr, condition, size=size, numsamples=numsamples, **kwargs)
+    # use subs when lambdify fails
+    except TypeError:
+        return sample_iter_subs(expr, condition, size=size, numsamples=numsamples, **kwargs)
 
 def quantile(expr, evaluate=True, **kwargs):
     r"""
@@ -1126,73 +1127,41 @@ def quantile(expr, evaluate=True, **kwargs):
     else:
         return result
 
-def sample_iter(expr, condition=None, size=(1,), library='scipy',
-                    numsamples=S.Infinity, **kwargs):
-
+def sample_iter_lambdify(expr, condition=None, size=(), numsamples=S.Infinity, **kwargs):
     """
-    Returns an iterator of realizations from the expression given a condition
-
-    Parameters
-    ==========
-
-    expr: Expr
-        Random expression to be realized
-    condition: Expr, optional
-        A conditional expression
-    size : int, tuple
-        Represents size of each sample in numsamples
-    numsamples: integer, optional
-        Length of the iterator (defaults to infinity)
-
-    Examples
-    ========
-
-    >>> from sympy.stats import Normal, sample_iter
-    >>> X = Normal('X', 0, 1)
-    >>> expr = X*X + 3
-    >>> iterator = sample_iter(expr, numsamples=3) # doctest: +SKIP
-    >>> list(iterator) # doctest: +SKIP
-    [12, 4, 7]
-
-    Returns
-    =======
-
-    sample_iter: iterator object
-        iterator object containing the sample/samples of given expr
+    Uses lambdify for computation. This is fast but does not always work.
 
     See Also
     ========
 
-    sample
-    sampling_P
-    sampling_E
+    sample_iter
 
     """
-    if not import_module(library):
-        raise ValueError("Failed to import %s" % library)
-
-    if condition is not None:
+    if condition:
         ps = pspace(Tuple(expr, condition))
     else:
         ps = pspace(expr)
 
     rvs = list(ps.values)
-    if library == 'pymc3':
-        # Currently unable to lambdify in pymc3
-        # TODO : Remove 'pymc3' when lambdify accepts 'pymc3' as module
-        fn = lambdify(rvs, expr, **kwargs)
-    else:
-        fn = lambdify(rvs, expr, modules=library, **kwargs)
-    if condition is not None:
+    fn = lambdify(rvs, expr, **kwargs)
+    if condition:
         given_fn = lambdify(rvs, condition, **kwargs)
 
+    # Check that lambdify can handle the expression
+    # Some operations like Sum can prove difficult
+    d = ps.sample(size)  # a dictionary that maps RVs to values
+    args = [d[rv] for rv in rvs]
+    fn(*args)
+    if condition:
+        given_fn(*args)
+
     def return_generator():
         count = 0
         while count < numsamples:
-            d = ps.sample(size=size, library=library)  # a dictionary that maps RVs to values
+            d = ps.sample(size)  # a dictionary that maps RVs to values
             args = [d[rv] for rv in rvs]
 
-            if condition is not None:  # Check that these values satisfy the condition
+            if condition:  # Check that these values satisfy the condition
                 gd = given_fn(*args)
                 if gd != True and gd != False:
                     raise ValueError(
@@ -1204,21 +1173,39 @@ def return_generator():
             count += 1
     return return_generator()
 
-def sample_iter_lambdify(expr, condition=None, size=(1,), numsamples=S.Infinity,
-                                                                    **kwargs):
 
-    return sample_iter(expr, condition=condition, size=size, numsamples=numsamples,
-                                                                        **kwargs)
+def sample_iter_subs(expr, condition=None, size=(), numsamples=S.Infinity, **kwargs):
+    """
+    Uses subs for computation. This is slow but almost always works.
+
+    See Also
+    ========
 
-def sample_iter_subs(expr, condition=None, size=(1,), numsamples=S.Infinity,
-                                                                    **kwargs):
+    sample_iter
 
-    return sample_iter(expr, condition=condition, size=size, numsamples=numsamples,
-                                                                        **kwargs)
+    """
+    if condition is not None:
+        ps = pspace(Tuple(expr, condition))
+    else:
+        ps = pspace(expr)
 
+    count = 0
+    while count < numsamples:
+        d = ps.sample(size)  # a dictionary that maps RVs to values
 
-def sampling_P(condition, given_condition=None, library='scipy', numsamples=1,
-                                                    evalf=True, **kwargs):
+        if condition is not None:  # Check that these values satisfy the condition
+            gd = condition.xreplace(d)
+            if gd != True and gd != False:
+                raise ValueError("Conditions must not contain free symbols")
+            if not gd:  # If the values don't satisfy then try again
+                continue
+
+        yield expr.xreplace(d)
+        count += 1
+
+
+def sampling_P(condition, given_condition=None, numsamples=1,
+               evalf=True, **kwargs):
     """
     Sampling version of P
 
@@ -1234,10 +1221,13 @@ def sampling_P(condition, given_condition=None, library='scipy', numsamples=1,
     count_true = 0
     count_false = 0
 
-    samples = sample_iter(condition, given_condition, library=library,
+    samples = sample_iter(condition, given_condition,
                           numsamples=numsamples, **kwargs)
 
     for sample in samples:
+        if sample != True and sample != False:
+            raise ValueError("Conditions must not contain free symbols")
+
         if sample:
             count_true += 1
         else:
@@ -1250,7 +1240,7 @@ def sampling_P(condition, given_condition=None, library='scipy', numsamples=1,
         return result
 
 
-def sampling_E(expr, given_condition=None, library='scipy', numsamples=1,
+def sampling_E(expr, given_condition=None, numsamples=1,
                evalf=True, **kwargs):
     """
     Sampling version of E
@@ -1262,20 +1252,17 @@ def sampling_E(expr, given_condition=None, library='scipy', numsamples=1,
     sampling_P
     sampling_density
     """
-    samples = list(sample_iter(expr, given_condition, library=library,
-                          numsamples=numsamples, **kwargs))
-    try:
-        result = Add(*[samp[0] for samp in samples]) / numsamples
-    except TypeError:
-        result = Add(*[samp for samp in samples]) / numsamples
 
+    samples = sample_iter(expr, given_condition,
+                          numsamples=numsamples, **kwargs)
+
+    result = Add(*list(samples)) / numsamples
     if evalf:
         return result.evalf()
     else:
         return result
 
-def sampling_density(expr, given_condition=None, library='scipy',
-                    numsamples=1, **kwargs):
+def sampling_density(expr, given_condition=None, numsamples=1, **kwargs):
     """
     Sampling version of density
 
@@ -1287,13 +1274,9 @@ def sampling_density(expr, given_condition=None, library='scipy',
     """
 
     results = {}
-    for result in sample_iter(expr, given_condition, library=library,
+    for result in sample_iter(expr, given_condition,
                               numsamples=numsamples, **kwargs):
-        try:
-            results[result[0]] = results.get(result[0], 0) + 1
-        except TypeError:
-            results[result] = results.get(result, 0) + 1
-
+        results[result] = results.get(result, 0) + 1
     return results
 
 


Test patch:
diff --git a/sympy/functions/special/tests/test_bsplines.py b/sympy/functions/special/tests/test_bsplines.py
index 11473fe367..6a838e3849 100644
--- a/sympy/functions/special/tests/test_bsplines.py
+++ b/sympy/functions/special/tests/test_bsplines.py
@@ -1,6 +1,5 @@
 from sympy.functions import bspline_basis_set, interpolating_spline
-from sympy import Piecewise, Interval, And
-from sympy import symbols, Rational, S
+from sympy import (Piecewise, Interval, And, symbols, Rational, S)
 from sympy.testing.pytest import slow
 
 x, y = symbols('x,y')
@@ -147,3 +146,17 @@ def test_6_points_degree_3():
     assert spline == Piecewise((6058*x**3/5301 - 18427*x**2/5301 + x*Rational(12622, 5301) + 3, (x >= -1) & (x <= 2)),
                                (-8327*x**3/5301 + 67883*x**2/5301 - x*Rational(159998, 5301) + Rational(43661, 1767), (x >= 2) & (x <= 3)),
                                (5414*x**3/47709 - 1386*x**2/589 + x*Rational(4267, 279) - Rational(12232, 589), (x >= 3) & (x <= 12)))
+
+
+def test_issue_19262():
+    Delta = symbols('Delta', positive=True)
+    knots = [i*Delta for i in range(4)]
+    basis = bspline_basis_set(1, knots, x)
+    y = symbols('y', nonnegative=True)
+    basis2 = bspline_basis_set(1, knots, y)
+    assert basis[0].subs(x, y) == basis2[0]
+    assert interpolating_spline(1, x,
+        [Delta*i for i in [1, 2, 4, 7]], [3, 6, 5, 7]
+        )  == Piecewise((3*x/Delta, (Delta <= x) & (x <= 2*Delta)),
+        (7 - x/(2*Delta), (x >= 2*Delta) & (x <= 4*Delta)),
+        (Rational(7, 3) + 2*x/(3*Delta), (x >= 4*Delta) & (x <= 7*Delta)))
diff --git a/sympy/polys/tests/test_factortools.py b/sympy/polys/tests/test_factortools.py
index 34fbe88267..7bcd935756 100644
--- a/sympy/polys/tests/test_factortools.py
+++ b/sympy/polys/tests/test_factortools.py
@@ -27,8 +27,7 @@ def test_dmp_trial_division():
 
 def test_dup_zz_mignotte_bound():
     R, x = ring("x", ZZ)
-    assert R.dup_zz_mignotte_bound(2*x**2 + 3*x + 4) == 6
-    assert R.dup_zz_mignotte_bound(x**3 + 14*x**2 + 56*x + 64) == 152
+    assert R.dup_zz_mignotte_bound(2*x**2 + 3*x + 4) == 32
 
 
 def test_dmp_zz_mignotte_bound():
diff --git a/sympy/stats/tests/test_continuous_rv.py b/sympy/stats/tests/test_continuous_rv.py
index 53391ebbed..35d17d1e29 100644
--- a/sympy/stats/tests/test_continuous_rv.py
+++ b/sympy/stats/tests/test_continuous_rv.py
@@ -24,7 +24,7 @@
 from sympy.stats.joint_rv_types import MultivariateLaplaceDistribution, MultivariateNormalDistribution
 from sympy.stats.crv import SingleContinuousPSpace, SingleContinuousDomain
 from sympy.stats.joint_rv import JointPSpace
-from sympy.testing.pytest import raises, XFAIL, slow, skip, ignore_warnings
+from sympy.testing.pytest import raises, XFAIL, slow, skip
 from sympy.testing.randtest import verify_numerically as tn
 
 oo = S.Infinity
@@ -321,15 +321,10 @@ def test_moment_generating_function():
 
 def test_sample_continuous():
     Z = ContinuousRV(z, exp(-z), set=Interval(0, oo))
-    assert density(Z)(-1) == 0
-
-    scipy = import_module('scipy')
-    if not scipy:
-        skip('Scipy is not installed. Abort tests')
-    with ignore_warnings(UserWarning):
-        assert next(sample(Z))[0] in Z.pspace.domain.set
+    assert sample(Z) in Z.pspace.domain.set
     sym, val = list(Z.pspace.sample().items())[0]
-    assert sym == Z and val[0] in Interval(0, oo)
+    assert sym == Z and val in Interval(0, oo)
+    assert density(Z)(-1) == 0
 
 
 def test_ContinuousRV():
@@ -736,8 +731,7 @@ def test_sampling_gamma_inverse():
     if not scipy:
         skip('Scipy not installed. Abort tests for sampling of gamma inverse.')
     X = GammaInverse("x", 1, 1)
-    with ignore_warnings(UserWarning):
-        assert next(sample(X))[0] in X.pspace.domain.set
+    assert sample(X) in X.pspace.domain.set
 
 def test_gompertz():
     b = Symbol("b", positive=True)
@@ -860,19 +854,15 @@ def test_lognormal():
     #assert variance(X) == (exp(std**2)-1) * exp(2*mean + std**2)
 
     # Right now, only density function and sampling works
-    scipy = import_module('scipy')
-    if not scipy:
-        skip('Scipy is not installed. Abort tests')
-    with ignore_warnings(UserWarning):
-        for i in range(3):
-            X = LogNormal('x', i, 1)
-            assert next(sample(X))[0] in X.pspace.domain.set
+
+    for i in range(3):
+        X = LogNormal('x', i, 1)
+        assert sample(X) in X.pspace.domain.set
 
     size = 5
-    with ignore_warnings(UserWarning):
-        samps = next(sample(X, size=size))
-        for samp in samps:
-            assert samp in X.pspace.domain.set
+    samps = sample(X, size=size)
+    for samp in samps:
+        assert samp in X.pspace.domain.set
     # The sympy integrator can't do this too well
     #assert E(X) ==
     raises(NotImplementedError, lambda: moment_generating_function(X))
@@ -975,8 +965,7 @@ def test_sampling_gaussian_inverse():
     if not scipy:
         skip('Scipy not installed. Abort tests for sampling of Gaussian inverse.')
     X = GaussianInverse("x", 1, 1)
-    with ignore_warnings(UserWarning):
-        assert next(sample(X, library='scipy'))[0] in X.pspace.domain.set
+    assert sample(X) in X.pspace.domain.set
 
 def test_pareto():
     xm, beta = symbols('xm beta', positive=True)
@@ -1262,9 +1251,6 @@ def test_wignersemicircle():
 
 
 def test_prefab_sampling():
-    scipy = import_module('scipy')
-    if not scipy:
-        skip('Scipy is not installed. Abort tests')
     N = Normal('X', 0, 1)
     L = LogNormal('L', 0, 1)
     E = Exponential('Ex', 1)
@@ -1277,13 +1263,12 @@ def test_prefab_sampling():
     variables = [N, L, E, P, W, U, B, G]
     niter = 10
     size = 5
-    with ignore_warnings(UserWarning):
-        for var in variables:
-            for i in range(niter):
-                assert next(sample(var))[0] in var.pspace.domain.set
-                samps = next(sample(var, size=size))
-                for samp in samps:
-                    assert samp in var.pspace.domain.set
+    for var in variables:
+        for i in range(niter):
+            assert sample(var) in var.pspace.domain.set
+            samps = sample(var, size=size)
+            for samp in samps:
+                assert samp in var.pspace.domain.set
 
 def test_input_value_assertions():
     a, b = symbols('a b')
@@ -1330,6 +1315,7 @@ def test_NormalDistribution():
     nd = NormalDistribution(0, 1)
     x = Symbol('x')
     assert nd.cdf(x) == erf(sqrt(2)*x/2)/2 + S.Half
+    assert isinstance(nd.sample(), float) or nd.sample().is_Number
     assert nd.expectation(1, x) == 1
     assert nd.expectation(x, x) == 0
     assert nd.expectation(x**2, x) == 1
@@ -1528,90 +1514,6 @@ def test_ContinuousDistributionHandmade():
     assert E(space.value) == Rational(3, 2)
     assert variance(space.value) == Rational(13, 12)
 
-
-def test_sample_numpy():
-    distribs_numpy = [
-        Beta("B", 1, 1),
-        Normal("N", 0, 1),
-        Gamma("G", 2, 7),
-        Exponential("E", 2),
-        LogNormal("LN", 0, 1),
-        Pareto("P", 1, 1),
-        ChiSquared("CS", 2),
-        Uniform("U", 0, 1)
-    ]
-    size = 3
-    numpy = import_module('numpy')
-    if not numpy:
-        skip('Numpy is not installed. Abort tests for _sample_numpy.')
-    else:
-        with ignore_warnings(UserWarning):
-            for X in distribs_numpy:
-                samps = next(sample(X, size=size, library='numpy'))
-                for sam in samps:
-                    assert sam in X.pspace.domain.set
-
-
-def test_sample_scipy():
-    distribs_scipy = [
-        Beta("B", 1, 1),
-        BetaPrime("BP", 1, 1),
-        Cauchy("C", 1, 1),
-        Chi("C", 1),
-        Normal("N", 0, 1),
-        Gamma("G", 2, 7),
-        GammaInverse("GI", 1, 1),
-        GaussianInverse("GUI", 1, 1),
-        Exponential("E", 2),
-        LogNormal("LN", 0, 1),
-        Pareto("P", 1, 1),
-        StudentT("S", 2),
-        ChiSquared("CS", 2),
-        Uniform("U", 0, 1)
-    ]
-    size = 3
-    numsamples = 5
-    scipy = import_module('scipy')
-    if not scipy:
-        skip('Scipy is not installed. Abort tests for _sample_scipy.')
-    else:
-        with ignore_warnings(UserWarning):
-            g_sample = list(sample(Gamma("G", 2, 7), size=size, numsamples=numsamples))
-            assert len(g_sample) == numsamples
-            for X in distribs_scipy:
-                samps = next(sample(X, size=size, library='scipy'))
-                samps2 = next(sample(X, size=(2, 2), library='scipy'))
-                for sam in samps:
-                    assert sam in X.pspace.domain.set
-                for i in range(2):
-                    for j in range(2):
-                        assert samps2[i][j] in X.pspace.domain.set
-
-def test_sample_pymc3():
-    distribs_pymc3 = [
-        Beta("B", 1, 1),
-        Cauchy("C", 1, 1),
-        Normal("N", 0, 1),
-        Gamma("G", 2, 7),
-        GaussianInverse("GI", 1, 1),
-        Exponential("E", 2),
-        LogNormal("LN", 0, 1),
-        Pareto("P", 1, 1),
-        ChiSquared("CS", 2),
-        Uniform("U", 0, 1)
-    ]
-    size = 3
-    pymc3 = import_module('pymc3')
-    if not pymc3:
-        skip('PyMC3 is not installed. Abort tests for _sample_pymc3.')
-    else:
-        with ignore_warnings(UserWarning):
-            for X in distribs_pymc3:
-                samps = next(sample(X, size=size, library='pymc3'))
-                for sam in samps:
-                    assert sam in X.pspace.domain.set
-
-
 def test_issue_16318():
     #test compute_expectation function of the SingleContinuousDomain
     N = SingleContinuousDomain(x, Interval(0, 1))
diff --git a/sympy/stats/tests/test_discrete_rv.py b/sympy/stats/tests/test_discrete_rv.py
index 2e495fe30c..ffef04fffb 100644
--- a/sympy/stats/tests/test_discrete_rv.py
+++ b/sympy/stats/tests/test_discrete_rv.py
@@ -13,7 +13,7 @@
                                     NegativeBinomial, Skellam, YuleSimon, Zeta,
                                     DiscreteRV)
 from sympy.stats.rv import sample
-from sympy.testing.pytest import slow, nocache_fail, raises, skip, ignore_warnings
+from sympy.testing.pytest import slow, nocache_fail, raises, skip
 from sympy.external import import_module
 
 x = Symbol('x')
@@ -137,14 +137,10 @@ def test_zeta():
 @slow
 def test_sample_discrete():
     X = Geometric('X', S.Half)
-    scipy = import_module('scipy')
-    if not scipy:
-        skip('Scipy not installed. Abort tests')
-    with ignore_warnings(UserWarning):
-        assert next(sample(X))[0] in X.pspace.domain.set
-        samps = next(sample(X, size=2)) # This takes long time if ran without scipy
-        for samp in samps:
-            assert samp in X.pspace.domain.set
+    assert sample(X) in X.pspace.domain.set
+    samps = sample(X, size=4)
+    for samp in samps:
+        assert samp in X.pspace.domain.set
 
 def test_discrete_probability():
     X = Geometric('X', Rational(1, 5))
@@ -307,26 +303,23 @@ def test_sampling_methods():
     if not numpy:
         skip('Numpy is not installed. Abort tests for _sample_numpy.')
     else:
-        with ignore_warnings(UserWarning):
-            for X in distribs_numpy:
-                samps = X.pspace.distribution._sample_numpy(size)
-                for samp in samps:
-                    assert samp in X.pspace.domain.set
+        for X in distribs_numpy:
+            samps = X.pspace.distribution._sample_numpy(size)
+            for samp in samps:
+                assert samp in X.pspace.domain.set
     scipy = import_module('scipy')
     if not scipy:
         skip('Scipy is not installed. Abort tests for _sample_scipy.')
     else:
-        with ignore_warnings(UserWarning):
-            for X in distribs_scipy:
-                samps = next(sample(X, size=size))
-                for samp in samps:
-                    assert samp in X.pspace.domain.set
+        for X in distribs_scipy:
+            samps = sample(X, size=size)
+            for samp in samps:
+                assert samp in X.pspace.domain.set
     pymc3 = import_module('pymc3')
     if not pymc3:
         skip('PyMC3 is not installed. Abort tests for _sample_pymc3.')
     else:
-        with ignore_warnings(UserWarning):
-            for X in distribs_pymc3:
-                samps = X.pspace.distribution._sample_pymc3(size)
-                for samp in samps:
-                    assert samp in X.pspace.domain.set
+        for X in distribs_pymc3:
+            samps = X.pspace.distribution._sample_pymc3(size)
+            for samp in samps:
+                assert samp in X.pspace.domain.set
diff --git a/sympy/stats/tests/test_finite_rv.py b/sympy/stats/tests/test_finite_rv.py
index 6c868c4774..178506c0d1 100644
--- a/sympy/stats/tests/test_finite_rv.py
+++ b/sympy/stats/tests/test_finite_rv.py
@@ -11,7 +11,7 @@
 from sympy.stats.frv_types import DieDistribution, BinomialDistribution, \
     HypergeometricDistribution
 from sympy.stats.rv import Density
-from sympy.testing.pytest import raises, skip, ignore_warnings
+from sympy.testing.pytest import raises, skip
 
 
 def BayesTest(A, B):
@@ -146,11 +146,7 @@ def test_given():
     X = Die('X', 6)
     assert density(X, X > 5) == {S(6): S.One}
     assert where(X > 2, X > 5).as_boolean() == Eq(X.symbol, 6)
-    scipy = import_module('scipy')
-    if not scipy:
-        skip('Scipy is not installed. Abort tests')
-    with ignore_warnings(UserWarning):
-        assert next(sample(X, X > 5)) == 6
+    assert sample(X, X > 5) == 6
 
 
 def test_domains():
@@ -449,27 +445,25 @@ def test_sampling_methods():
     distribs_pymc3 = [BetaBinomial("B", 1, 1, 1)]
 
     size = 5
-    with ignore_warnings(UserWarning):
-        for X in distribs_random:
-            sam = X.pspace.distribution._sample_random(size)
-            for i in range(size):
-                assert sam[i] in X.pspace.domain.set
+
+    for X in distribs_random:
+        sam = X.pspace.distribution._sample_random(size)
+        for i in range(size):
+            assert sam[i] in X.pspace.domain.set
 
     scipy = import_module('scipy')
     if not scipy:
         skip('Scipy not installed. Abort tests for _sample_scipy.')
     else:
-        with ignore_warnings(UserWarning):
-            for X in distribs_scipy:
-                sam = X.pspace.distribution._sample_scipy(size)
-                for i in range(size):
-                    assert sam[i] in X.pspace.domain.set
+        for X in distribs_scipy:
+            sam = X.pspace.distribution._sample_scipy(size)
+            for i in range(size):
+                assert sam[i] in X.pspace.domain.set
     pymc3 = import_module('pymc3')
     if not pymc3:
         skip('PyMC3 not installed. Abort tests for _sample_pymc3.')
     else:
-        with ignore_warnings(UserWarning):
-            for X in distribs_pymc3:
-                sam = X.pspace.distribution._sample_pymc3(size)
-                for i in range(size):
-                    assert sam[i] in X.pspace.domain.set
+        for X in distribs_pymc3:
+            sam = X.pspace.distribution._sample_pymc3(size)
+            for i in range(size):
+                assert sam[i] in X.pspace.domain.set
diff --git a/sympy/stats/tests/test_rv.py b/sympy/stats/tests/test_rv.py
index 6d91e3097d..2770ef570b 100644
--- a/sympy/stats/tests/test_rv.py
+++ b/sympy/stats/tests/test_rv.py
@@ -8,8 +8,7 @@
         DiscreteUniform, Poisson, characteristic_function, moment_generating_function)
 from sympy.stats.rv import (IndependentProductPSpace, rs_swap, Density, NamedArgsMixin,
         RandomSymbol, sample_iter, PSpace)
-from sympy.testing.pytest import raises, skip, XFAIL, ignore_warnings
-from sympy.external import import_module
+from sympy.testing.pytest import raises
 from sympy.core.numbers import comp
 from sympy.stats.frv_types import BernoulliDistribution
 
@@ -81,9 +80,6 @@ def test_sample_iter():
     Y = DiscreteUniform('Y', [1,2,7])
     Z = Poisson('Z', 2)
 
-    scipy = import_module('scipy')
-    if not scipy:
-        skip('Scipy is not installed. Abort tests')
     expr = X**2 + 3
     iterator = sample_iter(expr)
 
@@ -104,6 +100,7 @@ def is_iterator(obj):
             return True
         else:
             return False
+
     assert is_iterator(iterator)
     assert is_iterator(iterator2)
     assert is_iterator(iterator3)
@@ -183,38 +180,26 @@ def test_H():
 def test_Sample():
     X = Die('X', 6)
     Y = Normal('Y', 0, 1)
-    z = Symbol('z', integer=True)
+    z = Symbol('z')
 
-    scipy = import_module('scipy')
-    if not scipy:
-        skip('Scipy is not installed. Abort tests')
-    with ignore_warnings(UserWarning):
-        assert next(sample(X)) in [1, 2, 3, 4, 5, 6]
-        assert next(sample(X + Y))[0].is_Float
+    assert sample(X) in [1, 2, 3, 4, 5, 6]
+    assert sample(X + Y).is_Float
 
-    assert P(X + Y > 0, Y < 0, numsamples=10).is_number
+    P(X + Y > 0, Y < 0, numsamples=10).is_number
     assert E(X + Y, numsamples=10).is_number
-    assert E(X**2 + Y, numsamples=10).is_number
-    assert E((X + Y)**2, numsamples=10).is_number
     assert variance(X + Y, numsamples=10).is_number
 
-    raises(TypeError, lambda: P(Y > z, numsamples=5))
+    raises(ValueError, lambda: P(Y > z, numsamples=5))
 
     assert P(sin(Y) <= 1, numsamples=10) == 1
     assert P(sin(Y) <= 1, cos(Y) < 1, numsamples=10) == 1
 
-    assert all(i in range(1, 7) for i in density(X, numsamples=10))
-    assert all(i in range(4, 7) for i in density(X, X>3, numsamples=10))
+    # Make sure this doesn't raise an error
+    E(Sum(1/z**Y, (z, 1, oo)), Y > 2, numsamples=3)
 
 
-@XFAIL
-def test_samplingE():
-    scipy = import_module('scipy')
-    if not scipy:
-        skip('Scipy is not installed. Abort tests')
-    Y = Normal('Y', 0, 1)
-    z = Symbol('z', integer=True)
-    assert E(Sum(1/z**Y, (z, 1, oo)), Y > 2, numsamples=3).is_number
+    assert all(i in range(1, 7) for i in density(X, numsamples=10))
+    assert all(i in range(4, 7) for i in density(X, X>3, numsamples=10))
 
 
 def test_given():