diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 75353212..1aade107 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,6 +1,7 @@
 Master Branch
 =============
 
+
 Version 4.1.0 (2020-11-05)
 ==========================
 
diff --git a/chaospy/descriptives/conditional.py b/chaospy/descriptives/conditional.py
index 634a4d95..67d84288 100644
--- a/chaospy/descriptives/conditional.py
+++ b/chaospy/descriptives/conditional.py
@@ -2,6 +2,7 @@
 from itertools import product
 
 import numpy
+import chaospy
 import numpoly
 
 from . import expected
@@ -17,7 +18,7 @@ def E_cond(poly, freeze, dist, **kws):
     Args:
         poly (numpoly.ndpoly):
             Polynomial to find conditional expected value on.
-        freeze (numpy.ndarray):
+        freeze (numpy.ndpoly):
             Boolean values defining the conditional variables. True values
             implies that the value is conditioned on, e.g. frozen during the
             expected value calculation.
@@ -35,27 +36,31 @@ def E_cond(poly, freeze, dist, **kws):
         >>> poly
         polynomial([1, q0, q1, 10*q0*q1-1])
         >>> dist = chaospy.J(chaospy.Gamma(1, 1), chaospy.Normal(0, 2))
-        >>> chaospy.E_cond(poly, [1, 0], dist)
+        >>> chaospy.E_cond(poly, q0, dist)
         polynomial([1.0, q0, 0.0, -1.0])
-        >>> chaospy.E_cond(poly, [0, 1], dist)
+        >>> chaospy.E_cond(poly, q1, dist)
         polynomial([1.0, 1.0, q1, 10.0*q1-1.0])
-        >>> chaospy.E_cond(poly, [1, 1], dist)
+        >>> chaospy.E_cond(poly, [q0, q1], dist)
         polynomial([1, q0, q1, 10*q0*q1-1])
-        >>> chaospy.E_cond(poly, [0, 0], dist)
+        >>> chaospy.E_cond(poly, [], dist)
         polynomial([1.0, 1.0, 0.0, -1.0])
+        >>> chaospy.E_cond(4, [], dist)
+        array(4)
 
     """
     poly = numpoly.set_dimensions(poly, len(dist))
-    if not poly.isconstant:
+    if poly.isconstant():
         return poly.tonumpy()
     assert not dist.stochastic_dependent, dist
 
-    freeze = numpoly.polynomial(freeze)
-    if freeze.isconstant():
-        freeze = freeze.tonumpy().astype(bool)
+    freeze = numpoly.aspolynomial(freeze)
+    if not freeze.size:
+        return numpoly.polynomial(chaospy.E(poly, dist))
+    if not freeze.isconstant():
+        freeze = [name in freeze.names for name in poly.names]
     else:
-        poly, freeze = numpoly.align_exponents(poly, freeze)
-        freeze = numpy.isin(poly.keys, freeze.keys)
+        freeze = freeze.tonumpy()
+    freeze = numpy.asarray(freeze, dtype=bool)
 
     # decompose into frozen and unfrozen part
     poly = numpoly.decompose(poly)
@@ -71,5 +76,4 @@ def E_cond(poly, freeze, dist, **kws):
     # Remove frozen coefficients, such that poly == sum(frozen*unfrozen) holds
     for key in unfrozen.keys:
         unfrozen[key] = unfrozen[key] != 0
-
     return numpoly.sum(frozen*expected.E(unfrozen, dist), 0)
diff --git a/chaospy/descriptives/correlation/pearson.py b/chaospy/descriptives/correlation/pearson.py
index af2d521e..eaf90d14 100644
--- a/chaospy/descriptives/correlation/pearson.py
+++ b/chaospy/descriptives/correlation/pearson.py
@@ -40,9 +40,6 @@ def Corr(poly, dist=None, **kws):
     else:
         poly = numpoly.polynomial(poly)
 
-    if not poly.shape:
-        return numpy.ones((1, 1))
-
     cov = chaospy.Cov(poly, dist, **kws)
     var = numpy.diag(cov)
     vvar = numpy.sqrt(numpy.outer(var, var))
diff --git a/chaospy/descriptives/correlation/spearman.py b/chaospy/descriptives/correlation/spearman.py
index 076bb538..009feaab 100644
--- a/chaospy/descriptives/correlation/spearman.py
+++ b/chaospy/descriptives/correlation/spearman.py
@@ -1,8 +1,12 @@
 """Spearman's correlation coefficient."""
+import numpy
 from scipy.stats import spearmanr
 
+import numpoly
+import chaospy
 
-def Spearman(poly, dist, sample=10000, retall=False, **kws):
+
+def Spearman(poly, dist=None, sample=10000, retall=False, **kws):
     """
     Calculate Spearman's rank-order correlation coefficient.
 
@@ -24,10 +28,35 @@ def Spearman(poly, dist, sample=10000, retall=False, **kws):
             The two-sided p-value for a hypothesis test whose null hypothesis
             is that two sets of data are uncorrelated, has same dimension as
             ``rho``.
+
+    Examples:
+        >>> distribution = chaospy.MvNormal(
+        ...     [3, 4], [[2, .5], [.5, 1]])
+        >>> corr, pvalue = chaospy.Spearman(distribution, sample=50, retall=True)
+        >>> corr.round(4)
+        array([[1.   , 0.603],
+               [0.603, 1.   ]])
+        >>> pvalue.round(8)
+        array([[0.00e+00, 3.58e-06],
+               [3.58e-06, 0.00e+00]])
+
     """
+    if isinstance(poly, chaospy.Distribution):
+        poly, dist = numpoly.variable(len(poly)), poly
+    else:
+        poly = numpoly.polynomial(poly)
     samples = dist.sample(sample, **kws)
-    poly = polynomials.flatten(poly)
-    Y = poly(*samples)
+    corr = numpy.eye(len(poly))
+    pval = numpy.zeros((len(poly), len(poly)))
+    evals = poly.ravel()(*samples)
+    assert len(poly) == len(evals)
+    for idx in range(len(poly)):
+        for idy in range(idx+1, len(poly)):
+            if idx == idy:
+                pass
+            spear = spearmanr(evals[idx], evals[idy])
+            pval[idx, idy] = pval[idy, idx] = spear.pvalue
+            corr[idx, idy] = corr[idy, idx] = spear.correlation
     if retall:
-        return spearmanr(Y.T)
-    return spearmanr(Y.T)[0]
+        return corr, pval
+    return corr
diff --git a/chaospy/descriptives/covariance.py b/chaospy/descriptives/covariance.py
index 42c20c13..e6ba9ef3 100644
--- a/chaospy/descriptives/covariance.py
+++ b/chaospy/descriptives/covariance.py
@@ -33,13 +33,17 @@ def Cov(poly, dist=None, **kws):
                [  0. ,   2. ,   0.5,   0. ],
                [  0. ,   0.5,   1. ,   0. ],
                [  0. ,   0. ,   0. , 225. ]])
+        >>> chaospy.Cov([1, 2, 3], dist)
+        array([[0., 0., 0.],
+               [0., 0., 0.],
+               [0., 0., 0.]])
 
     """
     if dist is None:
         dist, poly = poly, numpoly.variable(len(poly))
     poly = numpoly.set_dimensions(poly, len(dist))
-    if not poly.isconstant:
-        return poly.tonumpy()**2
+    if poly.isconstant():
+        return numpy.zeros((len(poly), len(poly)))
     poly = poly-E(poly, dist)
     poly = numpoly.outer(poly, poly)
     return E(poly, dist)
diff --git a/chaospy/descriptives/kurtosis.py b/chaospy/descriptives/kurtosis.py
index d2bdadcf..0f9f9a3c 100644
--- a/chaospy/descriptives/kurtosis.py
+++ b/chaospy/descriptives/kurtosis.py
@@ -38,14 +38,17 @@ def Kurt(poly, dist=None, fisher=True, **kws):
         >>> poly = chaospy.polynomial([1, q0, q1, 10*q0*q1-1])
         >>> chaospy.Kurt(poly, dist).round(4)
         array([nan,  6.,  0., 15.])
+        >>> chaospy.Kurt(4., dist)
+        array(nan)
+
     """
     adjust = 3 if fisher else 0
 
     if dist is None:
         dist, poly = poly, numpoly.variable(len(poly))
     poly = numpoly.set_dimensions(poly, len(dist))
-    if not poly.isconstant:
-        return poly.tonumpy()**4-adjust
+    if poly.isconstant():
+        return numpy.full(poly.shape, numpy.nan)
 
     poly = poly-E(poly, dist, **kws)
     poly = numpoly.true_divide(poly, Std(poly, dist, **kws))
diff --git a/chaospy/descriptives/percentile.py b/chaospy/descriptives/percentile.py
index 3700dd4a..e563930c 100644
--- a/chaospy/descriptives/percentile.py
+++ b/chaospy/descriptives/percentile.py
@@ -39,14 +39,11 @@ def Perc(poly, q, dist, sample=10000, **kws):
     shape = poly.shape
     poly = poly.flatten()
 
-    q = numpy.array(q)/100.
+    q = numpy.asarray(q).ravel()/100.
     dim = len(dist)
 
     # Interior
-    Z = dist.sample(sample, **kws)
-    if dim==1:
-        Z = (Z,)
-        q = numpy.array([q])
+    Z = dist.sample(sample, **kws).reshape(len(dist), sample)
     poly1 = poly(*Z)
 
     # Min/max
diff --git a/chaospy/descriptives/skewness.py b/chaospy/descriptives/skewness.py
index 81bda767..b90dea1f 100644
--- a/chaospy/descriptives/skewness.py
+++ b/chaospy/descriptives/skewness.py
@@ -1,4 +1,5 @@
 """Skewness operator."""
+import numpy
 import numpoly
 
 from .expected import E
@@ -31,12 +32,15 @@ def Skew(poly, dist=None, **kws):
         >>> poly = chaospy.polynomial([1, q0, q1, 10*q0*q1-1])
         >>> chaospy.Skew(poly, dist)
         array([nan,  2.,  0.,  0.])
+        >>> chaospy.Skew(2., dist)
+        array(nan)
+
     """
     if dist is None:
         dist, poly = poly, numpoly.variable(len(poly))
     poly = numpoly.set_dimensions(poly, len(dist))
-    if not poly.isconstant:
-        return poly.tonumpy()**3
+    if poly.isconstant():
+        return numpy.full(poly.shape, numpy.nan)
 
     poly = poly-E(poly, dist, **kws)
     poly = numpoly.true_divide(poly, Std(poly, dist, **kws))
diff --git a/chaospy/descriptives/variance.py b/chaospy/descriptives/variance.py
index 9a04f846..b7e5a456 100644
--- a/chaospy/descriptives/variance.py
+++ b/chaospy/descriptives/variance.py
@@ -29,13 +29,15 @@ def Var(poly, dist=None, **kws):
         >>> poly = chaospy.polynomial([1, q0, q1, 10*q0*q1])
         >>> chaospy.Var(poly, dist)
         array([  0.,   1.,   4., 800.])
+        >>> chaospy.Var(2., dist)
+        array(0.)
+
     """
     if dist is None:
         dist, poly = poly, numpoly.variable(len(poly))
     poly = numpoly.set_dimensions(poly, len(dist))
-    if not poly.isconstant:
-
-        return poly.tonumpy()**2
+    if poly.isconstant():
+        return numpy.zeros(poly.shape)
     poly = poly-E(poly, dist, **kws)
     poly = numpoly.square(poly)
     return E(poly, dist, **kws)
diff --git a/chaospy/quadrature/gauss_patterson.py b/chaospy/quadrature/gauss_patterson.py
index f258606f..f93f71d6 100644
--- a/chaospy/quadrature/gauss_patterson.py
+++ b/chaospy/quadrature/gauss_patterson.py
@@ -60,11 +60,13 @@ def quad_gauss_patterson(order, domain):
 
     Example:
         >>> abscissas, weights = chaospy.quad_gauss_patterson(
-        ...     2, chaospy.Uniform(0, 1))
-        >>> abscissas.round(4)
-        array([[0.0198, 0.1127, 0.2829, 0.5   , 0.7171, 0.8873, 0.9802]])
-        >>> weights.round(4)
-        array([0.0523, 0.1342, 0.2007, 0.2255, 0.2007, 0.1342, 0.0523])
+        ...     1, chaospy.Iid(chaospy.Uniform(0, 1), 2))
+        >>> abscissas.round(3)
+        array([[0.113, 0.113, 0.113, 0.5  , 0.5  , 0.5  , 0.887, 0.887, 0.887],
+               [0.113, 0.5  , 0.887, 0.113, 0.5  , 0.887, 0.113, 0.5  , 0.887]])
+        >>> weights.round(3)
+        array([0.077, 0.123, 0.077, 0.123, 0.198, 0.123, 0.077, 0.123, 0.077])
+
     """
     if isinstance(domain, chaospy.Distribution):
         abscissas, weights = quad_gauss_patterson(
@@ -75,24 +77,22 @@ def quad_gauss_patterson(order, domain):
         weights /= numpy.sum(weights)
         return abscissas, weights
 
-    lower, upper = numpy.array(domain)
-    lower = numpy.asarray(lower).flatten()
-    upper = numpy.asarray(upper).flatten()
-
-    if len(lower) > 1:
+    domain = numpy.array(numpy.broadcast_arrays(*domain))
+    lower, upper = numpy.atleast_2d(domain.T).T
+    lower, upper, order = numpy.broadcast_arrays(lower, upper, order)
 
-        values = [quad_gauss_patterson(order[i], (domain[0][i], domain[1][i]))
-                  for i in range(len(domain[0]))]
+    if lower.size > 1:
+        values = [quad_gauss_patterson(order_, (lower_, upper_))
+                  for order_, lower_, upper_ in zip(order, lower, upper)]
 
-        abscissas = [_[0][0] for _ in values]
-        weights = [_[1] for _ in values]
+        abscissas = [value[0][0] for value in values]
+        weights = [value[1] for value in values]
         abscissas = combine(abscissas).T
         weights = numpy.prod(combine(weights), -1)
-
         return abscissas, weights
 
-    order = sorted(PATTERSON_VALUES.keys())[order]
-    abscissas, weights = PATTERSON_VALUES[order]
+    order = sorted(PATTERSON_VALUES.keys())[int(order)]
+    abscissas, weights = PATTERSON_VALUES[int(order)]
 
     abscissas = .5*(abscissas*(upper-lower)+upper+lower)
     abscissas = abscissas.reshape(1, abscissas.size)
diff --git a/chaospy/quadrature/genz_keister/frontend.py b/chaospy/quadrature/genz_keister/frontend.py
index df7b879c..25269faf 100644
--- a/chaospy/quadrature/genz_keister/frontend.py
+++ b/chaospy/quadrature/genz_keister/frontend.py
@@ -22,13 +22,15 @@ def quad_genz_keister(order, dist, rule=24):
     """
     Genz-Keister quadrature rule.
 
-    Eabsicassample:
+    Examples:
         >>> abscissas, weights = quad_genz_keister(
-        ...         order=1, dist=chaospy.Uniform(0, 1))
-        >>> abscissas.round(4)
-        array([[0.0416, 0.5   , 0.9584]])
-        >>> weights.round(4)
-        array([0.1667, 0.6667, 0.1667])
+        ...         order=1, dist=chaospy.Iid(chaospy.Uniform(0, 1), 2))
+        >>> abscissas.round(2)
+        array([[0.04, 0.04, 0.04, 0.5 , 0.5 , 0.5 , 0.96, 0.96, 0.96],
+               [0.04, 0.5 , 0.96, 0.04, 0.5 , 0.96, 0.04, 0.5 , 0.96]])
+        >>> weights.round(2)
+        array([0.03, 0.11, 0.03, 0.11, 0.44, 0.11, 0.03, 0.11, 0.03])
+
     """
     assert isinstance(rule, int)
 
@@ -51,5 +53,4 @@ def quad_genz_keister(order, dist, rule=24):
     abscissas, weights = foo(order)
     abscissas = dist.inv(scipy.special.ndtr(abscissas))
     abscissas = abscissas.reshape(1, abscissas.size)
-
     return abscissas, weights
diff --git a/chaospy/quadrature/genz_keister/gk16.py b/chaospy/quadrature/genz_keister/gk16.py
index 56712560..2b615bcc 100644
--- a/chaospy/quadrature/genz_keister/gk16.py
+++ b/chaospy/quadrature/genz_keister/gk16.py
@@ -372,7 +372,3 @@ def quad_genz_keister_16(order):
         1.8684014894510604E-18,
     )),
 }
-
-if __name__ == "__main__":
-    import doctest
-    doctest.testmod()
diff --git a/chaospy/quadrature/leja.py b/chaospy/quadrature/leja.py
index 28245954..b4f61434 100644
--- a/chaospy/quadrature/leja.py
+++ b/chaospy/quadrature/leja.py
@@ -55,21 +55,21 @@ def quad_leja(
                 The quadrature weights with ``weights.shape == (N,)``.
 
     Example:
-        >>> abscissas, weights = quad_leja(3, chaospy.Normal(0, 1))
-        >>> abscissas.round(4)
-        array([[-2.7173, -1.4142,  0.    ,  1.7635]])
-        >>> weights.round(4)
-        array([0.022 , 0.1629, 0.6506, 0.1645])
+        >>> abscissas, weights = quad_leja(
+        ...     2, chaospy.Iid(chaospy.Normal(0, 1), 2))
+        >>> abscissas.round(2)
+        array([[-1.41, -1.41, -1.41,  0.  ,  0.  ,  0.  ,  1.76,  1.76,  1.76],
+               [-1.41,  0.  ,  1.76, -1.41,  0.  ,  1.76, -1.41,  0.  ,  1.76]])
+        >>> weights.round(3)
+        array([0.05 , 0.133, 0.04 , 0.133, 0.359, 0.107, 0.04 , 0.107, 0.032])
+
     """
     if len(dist) > 1:
-        if dist.stochastic_depedent:
+        if dist.stochastic_dependent:
             raise chaospy.StochasticallyDependentError(
                 "Leja quadrature do not supper distribution with dependencies.")
-        if isinstance(order, int):
-            out = [quad_leja(order, _) for _ in dist]
-        else:
-            out = [quad_leja(order[_], dist[_]) for _ in range(len(dist))]
-
+        order = numpy.broadcast_to(order, len(dist))
+        out = [quad_leja(order[_], dist[_]) for _ in range(len(dist))]
         abscissas = [_[0][0] for _ in out]
         weights = [_[1] for _ in out]
         abscissas = combine(abscissas).T
@@ -90,11 +90,14 @@ def objective(abscissas_):
         def fmin(idx):
             """Bound minimization."""
             try:
-                x, fx = fminbound(objective, abscissas[idx], abscissas[idx+1], full_output=1)[:2]
-            except UnboundLocalError:
-                x = abscissas[idx] + 0.5*(3-5**0.5) * (abscissas[idx+1]-abscissas[idx])
-                fx = objective(x)
-            return x, fx
+                xopt, fval, _, _ = fminbound(objective, abscissas[idx],
+                                             abscissas[idx+1], full_output=True)
+            # Hard coded solution to scipy/scipy#11207 for scipy < 1.5.0.
+            except UnboundLocalError:  # pragma: no cover
+                xopt = abscissas[idx]+0.5*(3-5**0.5)*(
+                    abscissas[idx+1]-abscissas[idx])
+                fx = objective(xopt)
+            return xopt, fval
 
         opts, vals = zip(*[fmin(idx) for idx in range(len(abscissas)-1)])
         index = numpy.argmin(vals)
@@ -110,10 +113,10 @@ def fmin(idx):
 def create_weights(
         nodes,
         dist,
-        rule="fejer",
+        rule="clenshaw_curtis",
 ):
     """Create weights for the Laja method."""
-    _, poly, _ = chaospy.stieltjes(len(nodes)-1, dist)
-    poly = poly.flatten()
+    _, poly, _ = chaospy.stieltjes(len(nodes)-1, dist, rule=rule)
+    poly = poly.ravel()
     weights = numpy.linalg.inv(poly(nodes))
     return weights[:, 0]
diff --git a/chaospy/regression.py b/chaospy/regression.py
index d62dc27e..f7124fae 100644
--- a/chaospy/regression.py
+++ b/chaospy/regression.py
@@ -114,7 +114,7 @@ def fit_regression(
         evals = evals.reshape(len(evals), -1)
 
     if model is None:
-        uhat, _, _, _ = numpy.linalg.lstsq(poly_evals, evals)
+        uhat, _, _, _ = numpy.linalg.lstsq(poly_evals, evals, rcond=None)
 
     else:
         try: