Repository: mwaskom/seaborn

Issue:
Remove scipy dependency?
scipy is a common core package, but it is a source of some [installation difficulties](https://stackoverflow.com/questions/43549269/seaborn-importerror-dll-load-failed-the-specified-module-could-not-be-found) and is a rather heavy install that may not be desired in all data science contexts.

It would not be too hard to remove seaborn's dependency on it. The main things that would need to change are:

- Internal implementation of kernel density estimation (scipy's implementation is pretty simple, and would not be too hard to copy/adapt)
- Scipy's `integrate.trapz` algorithm is used in a few places for testing the kernel density normalization; this could be approximated
- `clustergrid` uses `scipy` for clustering. This could change to an optional dependency that is required only when `clustermap` is used, but I also periodically consider spinning `clustermap` off into its own package.
- Testing/documenting the `fit=` argument in `distplot` (but this has not be implemented in the distribution upgrades, and is a candidate for deprecation/removal
- `scipy.distance` is used for binning in `regplot`, but just using 1D euclidean distance; easy to do in numpy
- A few utility functions (e.g. `IQR`) that are very easy to implement locally.

Maybe a few other uses. It wouldn't be a trivial change, but it's something to think about.

Discussion:
codecov[bot]: # [Codecov](https://codecov.io/gh/mwaskom/seaborn/pull/2398?src=pr&el=h1) Report
> Merging [#2398](https://codecov.io/gh/mwaskom/seaborn/pull/2398?src=pr&el=desc) (88c86ee) into [master](https://codecov.io/gh/mwaskom/seaborn/commit/0a24478a550132f1882e5be5f5dbc0fc446a8a6c?el=desc) (0a24478) will **decrease** coverage by `0.32%`.
> The diff coverage is `55.55%`.

[![Impacted file tree graph](https://codecov.io/gh/mwaskom/seaborn/pull/2398/graphs/tree.svg?width=650&height=150&src=pr&token=SLuzHnA9Bl)](https://codecov.io/gh/mwaskom/seaborn/pull/2398?src=pr&el=tree)

```diff
@@            Coverage Diff             @@
##           master    #2398      +/-   ##
==========================================
- Coverage   97.02%   96.69%   -0.33%     
==========================================
  Files          17       17              
  Lines        6049     6026      -23     
==========================================
- Hits         5869     5827      -42     
- Misses        180      199      +19     
```


| [Impacted Files](https://codecov.io/gh/mwaskom/seaborn/pull/2398?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [seaborn/\_statistics.py](https://codecov.io/gh/mwaskom/seaborn/pull/2398/diff?src=pr&el=tree#diff-c2VhYm9ybi9fc3RhdGlzdGljcy5weQ==) | `97.97% <55.55%> (-2.03%)` | :arrow_down: |
| [seaborn/\_docstrings.py](https://codecov.io/gh/mwaskom/seaborn/pull/2398/diff?src=pr&el=tree#diff-c2VhYm9ybi9fZG9jc3RyaW5ncy5weQ==) | `97.05% <0.00%> (-2.95%)` | :arrow_down: |
| [seaborn/regression.py](https://codecov.io/gh/mwaskom/seaborn/pull/2398/diff?src=pr&el=tree#diff-c2VhYm9ybi9yZWdyZXNzaW9uLnB5) | `96.52% <0.00%> (-1.28%)` | :arrow_down: |
| [seaborn/distributions.py](https://codecov.io/gh/mwaskom/seaborn/pull/2398/diff?src=pr&el=tree#diff-c2VhYm9ybi9kaXN0cmlidXRpb25zLnB5) | `95.65% <0.00%> (-0.78%)` | :arrow_down: |
| [seaborn/categorical.py](https://codecov.io/gh/mwaskom/seaborn/pull/2398/diff?src=pr&el=tree#diff-c2VhYm9ybi9jYXRlZ29yaWNhbC5weQ==) | `97.81% <0.00%> (-0.18%)` | :arrow_down: |
| [seaborn/\_core.py](https://codecov.io/gh/mwaskom/seaborn/pull/2398/diff?src=pr&el=tree#diff-c2VhYm9ybi9fY29yZS5weQ==) | `98.11% <0.00%> (-0.17%)` | :arrow_down: |
| [seaborn/axisgrid.py](https://codecov.io/gh/mwaskom/seaborn/pull/2398/diff?src=pr&el=tree#diff-c2VhYm9ybi9heGlzZ3JpZC5weQ==) | `95.47% <0.00%> (-0.05%)` | :arrow_down: |
| [seaborn/matrix.py](https://codecov.io/gh/mwaskom/seaborn/pull/2398/diff?src=pr&el=tree#diff-c2VhYm9ybi9tYXRyaXgucHk=) | `95.78% <0.00%> (-0.02%)` | :arrow_down: |
| [seaborn/rcmod.py](https://codecov.io/gh/mwaskom/seaborn/pull/2398/diff?src=pr&el=tree#diff-c2VhYm9ybi9yY21vZC5weQ==) | `100.00% <0.00%> (ø)` | |
| [seaborn/\_decorators.py](https://codecov.io/gh/mwaskom/seaborn/pull/2398/diff?src=pr&el=tree#diff-c2VhYm9ybi9fZGVjb3JhdG9ycy5weQ==) | `100.00% <0.00%> (ø)` | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/mwaskom/seaborn/pull/2398?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/mwaskom/seaborn/pull/2398?src=pr&el=footer). Last update [0a24478...88c86ee](https://codecov.io/gh/mwaskom/seaborn/pull/2398?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).


Patch:
diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml
index 1e3a0049..e431cf18 100644
--- a/.github/workflows/ci.yaml
+++ b/.github/workflows/ci.yaml
@@ -25,8 +25,7 @@ jobs:
       - name: Install seaborn
         run: |
           python -m pip install --upgrade pip
-          pip install `cat ci/deps_latest.txt ci/utils.txt`
-          pip install .
+          pip install .[all] -r ci/utils.txt
 
       - name: Install doc tools
         run: |
@@ -47,20 +46,24 @@ jobs:
 
         python: [3.7.x, 3.8.x, 3.9.x]
         target: [test]
+        install: [all]
         deps: [latest]
         backend: [agg]
 
         include:
           - python: 3.7.x
             target: unittests
+            install: all
             deps: pinned
             backend: agg
           - python: 3.9.x
             target: unittests
-            deps: minimal
+            install: light
+            deps: latest
             backend: agg
           - python: 3.9.x
             target: test
+            install: all
             deps: latest
             backend: tkagg
 
@@ -77,8 +80,9 @@ jobs:
       - name: Install seaborn
         run: |
           python -m pip install --upgrade pip
-          pip install `cat ci/deps_${{ matrix.deps }}.txt ci/utils.txt`
-          pip install .
+          if [[ ${{matrix.install}} == 'all' ]]; then EXTRAS='[all]'; fi
+          if [[ ${{matrix.deps }} == 'pinned' ]]; then DEPS='-r ci/deps_pinned.txt'; fi
+          pip install .$EXTRAS $DEPS -r ci/utils.txt
 
       - name: Cache datastes
         run: python ci/cache_test_datasets.py
diff --git a/README.md b/README.md
index c1435c80..6a79e937 100644
--- a/README.md
+++ b/README.md
@@ -29,20 +29,29 @@ Dependencies
 
 Seaborn supports Python 3.7+ and no longer supports Python 2.
 
-Installation requires [numpy](https://numpy.org/), [scipy](https://www.scipy.org/), [pandas](https://pandas.pydata.org/), and [matplotlib](https://matplotlib.org/). Some functions will optionally use [statsmodels](https://www.statsmodels.org/) if it is installed.
+Installation requires [numpy](https://numpy.org/), [pandas](https://pandas.pydata.org/), and [matplotlib](https://matplotlib.org/). Some functions will optionally use [scipy](https://www.scipy.org/) and/or [statsmodels](https://www.statsmodels.org/) if they are available.
 
 
 Installation
 ------------
 
-The latest stable release (and older versions) can be installed from PyPI:
+The latest stable release (and required dependencies) can be installed from PyPI:
 
     pip install seaborn
 
+It is also possible to include the optional dependencies:
+
+    pip install seaborn[all]
+
 You may instead want to use the development version from Github:
 
-    pip install git+https://github.com/mwaskom/seaborn.git#egg=seaborn
+    pip install git+https://github.com/mwaskom/seaborn.git
+
+Seaborn is also available from Anaconda and can be installed with conda:
+
+    conda install seaborn
 
+Note that the main anaconda repository typically lags PyPI in adding new releases.
 
 Testing
 -------
diff --git a/ci/deps_minimal.txt b/ci/deps_minimal.txt
deleted file mode 100644
index 8fed0418..00000000
--- a/ci/deps_minimal.txt
+++ /dev/null
@@ -1,4 +0,0 @@
-numpy
-scipy
-matplotlib
-pandas
diff --git a/ci/deps_pinned.txt b/ci/deps_pinned.txt
index f71e59a8..f2d7f145 100644
--- a/ci/deps_pinned.txt
+++ b/ci/deps_pinned.txt
@@ -1,5 +1,5 @@
 numpy~=1.16.0
-scipy~=1.2.0
 pandas~=0.23.0
 matplotlib~=3.0.0
+scipy~=1.2.0
 statsmodels~=0.9.0
diff --git a/doc/installing.rst b/doc/installing.rst
index 7768cae9..80ed1e70 100644
--- a/doc/installing.rst
+++ b/doc/installing.rst
@@ -9,12 +9,17 @@ Installing and getting started
 
    <div class="col-md-9">
 
-Official releases of seaborn can be installed from `PyPI
-<https://pypi.org/project/seaborn/>`_::
+Official releases of seaborn can be installed from `PyPI <https://pypi.org/project/seaborn/>`_:
 
     pip install seaborn
 
-The library is also included as part of the `Anaconda <https://repo.anaconda.com/>`_ distribution::
+The basic invocation of `pip` will install seaborn and, if necessary, its mandatory dependencies.
+It is possible to include optional dependencies that give access to a few advanced features:
+
+    pip install seaborn[all]
+
+The library is also included as part of the `Anaconda <https://repo.anaconda.com/>`_ distribution,
+and it can be installed with `conda`:
 
     conda install seaborn
 
@@ -26,15 +31,11 @@ Supported Python versions
 
 - Python 3.7+
 
-Required dependencies
+Mandatory dependencies
 ^^^^^^^^^^^^^^^^^^^^^^
 
-If not already present, these libraries will be downloaded when you install seaborn.
-
 - `numpy <https://numpy.org/>`__
 
-- `scipy <https://www.scipy.org/>`__
-
 - `pandas <https://pandas.pydata.org/>`__
 
 - `matplotlib <https://matplotlib.org>`__
@@ -44,7 +45,9 @@ Optional dependencies
 
 - `statsmodels <https://www.statsmodels.org/>`__, for advanced regression plots
 
-- `fastcluster <https://pypi.org/project/fastcluster/>`__, for clustering large matrices
+- `scipy <https://www.scipy.org/>`__, for clustering matrices and some advanced options
+
+- `fastcluster <https://pypi.org/project/fastcluster/>`__, faster clustering of large matrices
 
 Quickstart
 ~~~~~~~~~~
diff --git a/doc/releases/v0.12.0.txt b/doc/releases/v0.12.0.txt
index 6642658b..f628ec21 100644
--- a/doc/releases/v0.12.0.txt
+++ b/doc/releases/v0.12.0.txt
@@ -2,4 +2,8 @@
 v0.12.0 (Unreleased)
 --------------------
 
-- Following `NEP29 <https://numpy.org/neps/nep-0029-deprecation_policy.html>`_, dropped support for Python 3.6 and bumped the minimally-supported versions of the library dependencies.
\ No newline at end of file
+- Made `scipy` an optional dependency and added `pip install seaborn[all]` as a method for ensuring the availability of compatible `scipy` and `statsmodels` libraries. This has a few minor implications for existing code, which are explained in the Github pull request (:pr:`2398`).
+
+- Following `NEP29 <https://numpy.org/neps/nep-0029-deprecation_policy.html>`_, dropped support for Python 3.6 and bumped the minimally-supported versions of the library dependencies.
+
+- Removed several previously-deprecated utility functions (`iqr`, `percentiles`, `pmf_hist`, and `sort_df`).
\ No newline at end of file
diff --git a/requirements.txt b/requirements.txt
deleted file mode 100644
index 8fed0418..00000000
--- a/requirements.txt
+++ /dev/null
@@ -1,4 +0,0 @@
-numpy
-scipy
-matplotlib
-pandas
diff --git a/seaborn/_statistics.py b/seaborn/_statistics.py
index b76f397e..951375d7 100644
--- a/seaborn/_statistics.py
+++ b/seaborn/_statistics.py
@@ -26,7 +26,12 @@ The classes should behave roughly in the style of scikit-learn.
 """
 from numbers import Number
 import numpy as np
-from scipy import stats
+try:
+    from scipy.stats import gaussian_kde
+    _no_scipy = False
+except ImportError:
+    from .external.kde import gaussian_kde
+    _no_scipy = True
 
 from .utils import _check_argument
 
@@ -61,7 +66,7 @@ class KDE:
         clip : pair of numbers None, or a pair of such pairs
             Do not evaluate the density outside of these limits.
         cumulative : bool, optional
-            If True, estimate a cumulative distribution function.
+            If True, estimate a cumulative distribution function. Requires scipy.
 
         """
         if clip is None:
@@ -74,6 +79,9 @@ class KDE:
         self.clip = clip
         self.cumulative = cumulative
 
+        if cumulative and _no_scipy:
+            raise RuntimeError("Cumulative KDE evaluation requires scipy")
+
         self.support = None
 
     def _define_support_grid(self, x, bw, cut, clip, gridsize):
@@ -129,7 +137,7 @@ class KDE:
         if weights is not None:
             fit_kws["weights"] = weights
 
-        kde = stats.gaussian_kde(fit_data, **fit_kws)
+        kde = gaussian_kde(fit_data, **fit_kws)
         kde.set_bandwidth(kde.factor * self.bw_adjust)
 
         return kde
diff --git a/seaborn/categorical.py b/seaborn/categorical.py
index c060be4d..76b012e0 100644
--- a/seaborn/categorical.py
+++ b/seaborn/categorical.py
@@ -1,18 +1,26 @@
 from textwrap import dedent
 from numbers import Number
+import warnings
 import colorsys
+from functools import partial
+
 import numpy as np
-from scipy import stats
 import pandas as pd
+try:
+    from scipy.stats import gaussian_kde
+    _no_scipy = False
+except ImportError:
+    from .external.kde import gaussian_kde
+    _no_scipy = True
+
 import matplotlib as mpl
 from matplotlib.collections import PatchCollection
 import matplotlib.patches as Patches
 import matplotlib.pyplot as plt
-import warnings
 
 from ._core import variable_type, infer_orient, categorical_order
 from . import utils
-from .utils import remove_na
+from .utils import remove_na, _normal_quantile_func
 from .algorithms import bootstrap
 from .palettes import color_palette, husl_palette, light_palette, dark_palette
 from .axisgrid import FacetGrid, _facet_docs
@@ -662,7 +670,7 @@ class _ViolinPlotter(_CategoricalPlotter):
 
     def fit_kde(self, x, bw):
         """Estimate a KDE for a vector of data with flexible bandwidth."""
-        kde = stats.gaussian_kde(x, bw)
+        kde = gaussian_kde(x, bw)
 
         # Extract the numeric bandwidth from the KDE object
         bw_used = kde.factor
@@ -942,7 +950,7 @@ class _ViolinPlotter(_CategoricalPlotter):
         """Draw boxplot information at center of the density."""
         # Compute the boxplot statistics
         q25, q50, q75 = np.percentile(data, [25, 50, 75])
-        whisker_lim = 1.5 * stats.iqr(data)
+        whisker_lim = 1.5 * (q75 - q25)
         h1 = np.min(data[data >= (q25 - whisker_lim)])
         h2 = np.max(data[data <= (q75 + whisker_lim)])
 
@@ -1099,7 +1107,7 @@ class _StripPlotter(_CategoricalScatterPlotter):
             jlim = float(jitter)
         if self.hue_names is not None and dodge:
             jlim /= len(self.hue_names)
-        self.jitterer = stats.uniform(-jlim, jlim * 2).rvs
+        self.jitterer = partial(np.random.uniform, low=-jlim, high=+jlim)
 
     def draw_stripplot(self, ax, kws):
         """Draw the points onto `ax`."""
@@ -1120,7 +1128,7 @@ class _StripPlotter(_CategoricalScatterPlotter):
 
                 # Plot the points in centered positions
                 cat_pos = np.ones(strip_data.size) * i
-                cat_pos += self.jitterer(len(strip_data))
+                cat_pos += self.jitterer(size=len(strip_data))
                 kws.update(c=palette[point_colors])
                 if self.orient == "v":
                     ax.scatter(cat_pos, strip_data, **kws)
@@ -1138,7 +1146,7 @@ class _StripPlotter(_CategoricalScatterPlotter):
                     # Plot the points in centered positions
                     center = i + offsets[j]
                     cat_pos = np.ones(strip_data.size) * center
-                    cat_pos += self.jitterer(len(strip_data))
+                    cat_pos += self.jitterer(size=len(strip_data))
                     kws.update(c=palette[point_colors])
                     if self.orient == "v":
                         ax.scatter(cat_pos, strip_data, **kws)
@@ -1846,7 +1854,7 @@ class _LVPlotter(_CategoricalPlotter):
         elif self.k_depth == 'proportion':
             k = int(np.log2(n)) - int(np.log2(n * p)) + 1
         elif self.k_depth == 'trustworthy':
-            point_conf = 2 * stats.norm.ppf((1 - self.trust_alpha / 2)) ** 2
+            point_conf = 2 * _normal_quantile_func((1 - self.trust_alpha / 2)) ** 2
             k = int(np.log2(n / point_conf)) + 1
         else:
             k = int(self.k_depth)  # allow having k as input
diff --git a/seaborn/distributions.py b/seaborn/distributions.py
index 9e17f1bf..23772d51 100644
--- a/seaborn/distributions.py
+++ b/seaborn/distributions.py
@@ -11,7 +11,6 @@ import matplotlib.pyplot as plt
 import matplotlib.transforms as tx
 from matplotlib.colors import to_rgba
 from matplotlib.collections import LineCollection
-from scipy import stats
 
 from ._core import (
     VectorPlotter,
@@ -34,6 +33,7 @@ from .utils import (
 )
 from .palettes import color_palette
 from .external import husl
+from .external.kde import gaussian_kde
 from ._decorators import _deprecate_positional_args
 from ._docstrings import (
     DocstringComponents,
@@ -2395,7 +2395,8 @@ def _freedman_diaconis_bins(a):
     a = np.asarray(a)
     if len(a) < 2:
         return 1
-    h = 2 * stats.iqr(a) / (len(a) ** (1 / 3))
+    iqr = np.subtract.reduce(np.nanpercentile(a, [75, 25]))
+    h = 2 * iqr / (len(a) ** (1 / 3))
     # fall back to sqrt(a) bins if iqr is 0
     if h == 0:
         return int(np.sqrt(a.size))
@@ -2642,7 +2643,7 @@ def distplot(a=None, bins=None, hist=True, kde=True, rug=False, fit=None,
         gridsize = fit_kws.pop("gridsize", 200)
         cut = fit_kws.pop("cut", 3)
         clip = fit_kws.pop("clip", (-np.inf, np.inf))
-        bw = stats.gaussian_kde(a).scotts_factor() * a.std(ddof=1)
+        bw = gaussian_kde(a).scotts_factor() * a.std(ddof=1)
         x = _kde_support(a, bw, gridsize, cut, clip)
         params = fit.fit(a)
         y = pdf(x)
diff --git a/seaborn/external/kde.py b/seaborn/external/kde.py
new file mode 100644
index 00000000..1b1d04db
--- /dev/null
+++ b/seaborn/external/kde.py
@@ -0,0 +1,382 @@
+"""
+This module was copied from the scipy project.
+
+In the process of copying, some methods were removed because they depended on
+other parts of scipy (especially on compiled components), allowing seaborn to
+have a simple and pure Python implementation. These include:
+
+- integrate_gaussian
+- integrate_box
+- integrate_box_1d
+- integrate_kde
+- logpdf
+- resample
+
+Additionally, the numpy.linalg module was subsituted for scipy.linalg,
+and the examples section (with doctests) was removed from the docstring
+
+The original scipy license is copied below:
+
+Copyright (c) 2001-2002 Enthought, Inc.  2003-2019, SciPy Developers.
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+
+1. Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above
+   copyright notice, this list of conditions and the following
+   disclaimer in the documentation and/or other materials provided
+   with the distribution.
+
+3. Neither the name of the copyright holder nor the names of its
+   contributors may be used to endorse or promote products derived
+   from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+"""
+
+# -------------------------------------------------------------------------------
+#
+#  Define classes for (uni/multi)-variate kernel density estimation.
+#
+#  Currently, only Gaussian kernels are implemented.
+#
+#  Written by: Robert Kern
+#
+#  Date: 2004-08-09
+#
+#  Modified: 2005-02-10 by Robert Kern.
+#              Contributed to SciPy
+#            2005-10-07 by Robert Kern.
+#              Some fixes to match the new scipy_core
+#
+#  Copyright 2004-2005 by Enthought, Inc.
+#
+# -------------------------------------------------------------------------------
+
+import numpy as np
+from numpy import (asarray, atleast_2d, reshape, zeros, newaxis, dot, exp, pi,
+                   sqrt, ravel, power, atleast_1d, squeeze, sum, transpose,
+                   ones, cov)
+from numpy import linalg
+
+
+__all__ = ['gaussian_kde']
+
+
+class gaussian_kde(object):
+    """Representation of a kernel-density estimate using Gaussian kernels.
+
+    Kernel density estimation is a way to estimate the probability density
+    function (PDF) of a random variable in a non-parametric way.
+    `gaussian_kde` works for both uni-variate and multi-variate data.   It
+    includes automatic bandwidth determination.  The estimation works best for
+    a unimodal distribution; bimodal or multi-modal distributions tend to be
+    oversmoothed.
+
+    Parameters
+    ----------
+    dataset : array_like
+        Datapoints to estimate from. In case of univariate data this is a 1-D
+        array, otherwise a 2-D array with shape (# of dims, # of data).
+    bw_method : str, scalar or callable, optional
+        The method used to calculate the estimator bandwidth.  This can be
+        'scott', 'silverman', a scalar constant or a callable.  If a scalar,
+        this will be used directly as `kde.factor`.  If a callable, it should
+        take a `gaussian_kde` instance as only parameter and return a scalar.
+        If None (default), 'scott' is used.  See Notes for more details.
+    weights : array_like, optional
+        weights of datapoints. This must be the same shape as dataset.
+        If None (default), the samples are assumed to be equally weighted
+
+    Attributes
+    ----------
+    dataset : ndarray
+        The dataset with which `gaussian_kde` was initialized.
+    d : int
+        Number of dimensions.
+    n : int
+        Number of datapoints.
+    neff : int
+        Effective number of datapoints.
+
+        .. versionadded:: 1.2.0
+    factor : float
+        The bandwidth factor, obtained from `kde.covariance_factor`, with which
+        the covariance matrix is multiplied.
+    covariance : ndarray
+        The covariance matrix of `dataset`, scaled by the calculated bandwidth
+        (`kde.factor`).
+    inv_cov : ndarray
+        The inverse of `covariance`.
+
+    Methods
+    -------
+    evaluate
+    __call__
+    integrate_gaussian
+    integrate_box_1d
+    integrate_box
+    integrate_kde
+    pdf
+    logpdf
+    resample
+    set_bandwidth
+    covariance_factor
+
+    Notes
+    -----
+    Bandwidth selection strongly influences the estimate obtained from the KDE
+    (much more so than the actual shape of the kernel).  Bandwidth selection
+    can be done by a "rule of thumb", by cross-validation, by "plug-in
+    methods" or by other means; see [3]_, [4]_ for reviews.  `gaussian_kde`
+    uses a rule of thumb, the default is Scott's Rule.
+
+    Scott's Rule [1]_, implemented as `scotts_factor`, is::
+
+        n**(-1./(d+4)),
+
+    with ``n`` the number of data points and ``d`` the number of dimensions.
+    In the case of unequally weighted points, `scotts_factor` becomes::
+
+        neff**(-1./(d+4)),
+
+    with ``neff`` the effective number of datapoints.
+    Silverman's Rule [2]_, implemented as `silverman_factor`, is::
+
+        (n * (d + 2) / 4.)**(-1. / (d + 4)).
+
+    or in the case of unequally weighted points::
+
+        (neff * (d + 2) / 4.)**(-1. / (d + 4)).
+
+    Good general descriptions of kernel density estimation can be found in [1]_
+    and [2]_, the mathematics for this multi-dimensional implementation can be
+    found in [1]_.
+
+    With a set of weighted samples, the effective number of datapoints ``neff``
+    is defined by::
+
+        neff = sum(weights)^2 / sum(weights^2)
+
+    as detailed in [5]_.
+
+    References
+    ----------
+    .. [1] D.W. Scott, "Multivariate Density Estimation: Theory, Practice, and
+           Visualization", John Wiley & Sons, New York, Chicester, 1992.
+    .. [2] B.W. Silverman, "Density Estimation for Statistics and Data
+           Analysis", Vol. 26, Monographs on Statistics and Applied Probability,
+           Chapman and Hall, London, 1986.
+    .. [3] B.A. Turlach, "Bandwidth Selection in Kernel Density Estimation: A
+           Review", CORE and Institut de Statistique, Vol. 19, pp. 1-33, 1993.
+    .. [4] D.M. Bashtannyk and R.J. Hyndman, "Bandwidth selection for kernel
+           conditional density estimation", Computational Statistics & Data
+           Analysis, Vol. 36, pp. 279-298, 2001.
+    .. [5] Gray P. G., 1969, Journal of the Royal Statistical Society.
+           Series A (General), 132, 272
+
+    """
+    def __init__(self, dataset, bw_method=None, weights=None):
+        self.dataset = atleast_2d(asarray(dataset))
+        if not self.dataset.size > 1:
+            raise ValueError("`dataset` input should have multiple elements.")
+
+        self.d, self.n = self.dataset.shape
+
+        if weights is not None:
+            self._weights = atleast_1d(weights).astype(float)
+            self._weights /= sum(self._weights)
+            if self.weights.ndim != 1:
+                raise ValueError("`weights` input should be one-dimensional.")
+            if len(self._weights) != self.n:
+                raise ValueError("`weights` input should be of length n")
+            self._neff = 1/sum(self._weights**2)
+
+        self.set_bandwidth(bw_method=bw_method)
+
+    def evaluate(self, points):
+        """Evaluate the estimated pdf on a set of points.
+
+        Parameters
+        ----------
+        points : (# of dimensions, # of points)-array
+            Alternatively, a (# of dimensions,) vector can be passed in and
+            treated as a single point.
+
+        Returns
+        -------
+        values : (# of points,)-array
+            The values at each point.
+
+        Raises
+        ------
+        ValueError : if the dimensionality of the input points is different than
+                     the dimensionality of the KDE.
+
+        """
+        points = atleast_2d(asarray(points))
+
+        d, m = points.shape
+        if d != self.d:
+            if d == 1 and m == self.d:
+                # points was passed in as a row vector
+                points = reshape(points, (self.d, 1))
+                m = 1
+            else:
+                msg = "points have dimension %s, dataset has dimension %s" % (d,
+                    self.d)
+                raise ValueError(msg)
+
+        output_dtype = np.common_type(self.covariance, points)
+        result = zeros((m,), dtype=output_dtype)
+
+        whitening = linalg.cholesky(self.inv_cov)
+        scaled_dataset = dot(whitening, self.dataset)
+        scaled_points = dot(whitening, points)
+
+        if m >= self.n:
+            # there are more points than data, so loop over data
+            for i in range(self.n):
+                diff = scaled_dataset[:, i, newaxis] - scaled_points
+                energy = sum(diff * diff, axis=0) / 2.0
+                result += self.weights[i]*exp(-energy)
+        else:
+            # loop over points
+            for i in range(m):
+                diff = scaled_dataset - scaled_points[:, i, newaxis]
+                energy = sum(diff * diff, axis=0) / 2.0
+                result[i] = sum(exp(-energy)*self.weights, axis=0)
+
+        result = result / self._norm_factor
+
+        return result
+
+    __call__ = evaluate
+
+    def scotts_factor(self):
+        """Compute Scott's factor.
+
+        Returns
+        -------
+        s : float
+            Scott's factor.
+        """
+        return power(self.neff, -1./(self.d+4))
+
+    def silverman_factor(self):
+        """Compute the Silverman factor.
+
+        Returns
+        -------
+        s : float
+            The silverman factor.
+        """
+        return power(self.neff*(self.d+2.0)/4.0, -1./(self.d+4))
+
+    #  Default method to calculate bandwidth, can be overwritten by subclass
+    covariance_factor = scotts_factor
+    covariance_factor.__doc__ = """Computes the coefficient (`kde.factor`) that
+        multiplies the data covariance matrix to obtain the kernel covariance
+        matrix. The default is `scotts_factor`.  A subclass can overwrite this
+        method to provide a different method, or set it through a call to
+        `kde.set_bandwidth`."""
+
+    def set_bandwidth(self, bw_method=None):
+        """Compute the estimator bandwidth with given method.
+
+        The new bandwidth calculated after a call to `set_bandwidth` is used
+        for subsequent evaluations of the estimated density.
+
+        Parameters
+        ----------
+        bw_method : str, scalar or callable, optional
+            The method used to calculate the estimator bandwidth.  This can be
+            'scott', 'silverman', a scalar constant or a callable.  If a
+            scalar, this will be used directly as `kde.factor`.  If a callable,
+            it should take a `gaussian_kde` instance as only parameter and
+            return a scalar.  If None (default), nothing happens; the current
+            `kde.covariance_factor` method is kept.
+
+        Notes
+        -----
+        .. versionadded:: 0.11
+
+        """
+        if bw_method is None:
+            pass
+        elif bw_method == 'scott':
+            self.covariance_factor = self.scotts_factor
+        elif bw_method == 'silverman':
+            self.covariance_factor = self.silverman_factor
+        elif np.isscalar(bw_method) and not isinstance(bw_method, str):
+            self._bw_method = 'use constant'
+            self.covariance_factor = lambda: bw_method
+        elif callable(bw_method):
+            self._bw_method = bw_method
+            self.covariance_factor = lambda: self._bw_method(self)
+        else:
+            msg = "`bw_method` should be 'scott', 'silverman', a scalar " \
+                  "or a callable."
+            raise ValueError(msg)
+
+        self._compute_covariance()
+
+    def _compute_covariance(self):
+        """Computes the covariance matrix for each Gaussian kernel using
+        covariance_factor().
+        """
+        self.factor = self.covariance_factor()
+        # Cache covariance and inverse covariance of the data
+        if not hasattr(self, '_data_inv_cov'):
+            self._data_covariance = atleast_2d(cov(self.dataset, rowvar=1,
+                                               bias=False,
+                                               aweights=self.weights))
+            self._data_inv_cov = linalg.inv(self._data_covariance)
+
+        self.covariance = self._data_covariance * self.factor**2
+        self.inv_cov = self._data_inv_cov / self.factor**2
+        self._norm_factor = sqrt(linalg.det(2*pi*self.covariance))
+
+    def pdf(self, x):
+        """
+        Evaluate the estimated pdf on a provided set of points.
+
+        Notes
+        -----
+        This is an alias for `gaussian_kde.evaluate`.  See the ``evaluate``
+        docstring for more details.
+
+        """
+        return self.evaluate(x)
+
+    @property
+    def weights(self):
+        try:
+            return self._weights
+        except AttributeError:
+            self._weights = ones(self.n)/self.n
+            return self._weights
+
+    @property
+    def neff(self):
+        try:
+            return self._neff
+        except AttributeError:
+            self._neff = 1/sum(self.weights**2)
+            return self._neff
diff --git a/seaborn/matrix.py b/seaborn/matrix.py
index 00398038..e3cd6119 100644
--- a/seaborn/matrix.py
+++ b/seaborn/matrix.py
@@ -8,7 +8,11 @@ import matplotlib.pyplot as plt
 from matplotlib import gridspec
 import numpy as np
 import pandas as pd
-from scipy.cluster import hierarchy
+try:
+    from scipy.cluster import hierarchy
+    _no_scipy = False
+except ImportError:
+    _no_scipy = True
 
 from . import cm
 from .axisgrid import Grid
@@ -771,6 +775,9 @@ def dendrogram(
     dendrogramplotter.reordered_ind
 
     """
+    if _no_scipy:
+        raise RuntimeError("dendrogram requires scipy to be installed")
+
     plotter = _DendrogramPlotter(data, linkage=linkage, axis=axis,
                                  metric=metric, method=method,
                                  label=label, rotate=rotate)
@@ -786,6 +793,8 @@ class ClusterGrid(Grid):
                  figsize=None, row_colors=None, col_colors=None, mask=None,
                  dendrogram_ratio=None, colors_ratio=None, cbar_pos=None):
         """Grid object for organizing clustered heatmap input on to axes"""
+        if _no_scipy:
+            raise RuntimeError("ClusterGrid requires scipy to be available")
 
         if isinstance(data, pd.DataFrame):
             self.data = data
@@ -1250,6 +1259,8 @@ def clustermap(
     """
     Plot a matrix dataset as a hierarchically-clustered heatmap.
 
+    This function requires scipy to be available.
+
     Parameters
     ----------
     data : 2D array-like
@@ -1399,6 +1410,9 @@ def clustermap(
 
         >>> g = sns.clustermap(iris, z_score=0, cmap="vlag")
     """
+    if _no_scipy:
+        raise RuntimeError("clustermap requires scipy to be available")
+
     plotter = ClusterGrid(data, pivot_kws=pivot_kws, figsize=figsize,
                           row_colors=row_colors, col_colors=col_colors,
                           z_score=z_score, standard_scale=standard_scale,
diff --git a/seaborn/regression.py b/seaborn/regression.py
index 9f9b35c2..e302149b 100644
--- a/seaborn/regression.py
+++ b/seaborn/regression.py
@@ -4,7 +4,6 @@ from textwrap import dedent
 import warnings
 import numpy as np
 import pandas as pd
-from scipy.spatial import distance
 import matplotlib as mpl
 import matplotlib.pyplot as plt
 
@@ -318,17 +317,17 @@ class _RegressionPlotter(_LinearPlotter):
 
     def bin_predictor(self, bins):
         """Discretize a predictor by assigning value to closest bin."""
-        x = self.x
+        x = np.asarray(self.x)
         if np.isscalar(bins):
             percentiles = np.linspace(0, 100, bins + 2)[1:-1]
-            bins = np.c_[np.percentile(x, percentiles)]
+            bins = np.percentile(x, percentiles)
         else:
-            bins = np.c_[np.ravel(bins)]
+            bins = np.ravel(bins)
 
-        dist = distance.cdist(np.c_[x], bins)
+        dist = np.abs(np.subtract.outer(x, bins))
         x_binned = bins[np.argmin(dist, axis=1)].ravel()
 
-        return x_binned, bins.ravel()
+        return x_binned, bins
 
     def regress_out(self, a, b):
         """Regress b from a keeping a's original mean."""
diff --git a/seaborn/utils.py b/seaborn/utils.py
index 67cd4703..2b5cce5a 100644
--- a/seaborn/utils.py
+++ b/seaborn/utils.py
@@ -7,7 +7,6 @@ import colorsys
 from urllib.request import urlopen, urlretrieve
 
 import numpy as np
-from scipy import stats
 import pandas as pd
 import matplotlib as mpl
 import matplotlib.colors as mplcol
@@ -19,16 +18,6 @@ __all__ = ["desaturate", "saturate", "set_hls_values",
            "despine", "get_dataset_names", "get_data_home", "load_dataset"]
 
 
-def sort_df(df, *args, **kwargs):
-    """Wrapper to handle different pandas sorting API pre/post 0.17."""
-    msg = "This function is deprecated and will be removed in a future version"
-    warnings.warn(msg)
-    try:
-        return df.sort_values(*args, **kwargs)
-    except AttributeError:
-        return df.sort(*args, **kwargs)
-
-
 def ci_to_errsize(cis, heights):
     """Convert intervals to error arguments relative to plot heights.
 
@@ -59,34 +48,27 @@ def ci_to_errsize(cis, heights):
     return errsize
 
 
-def pmf_hist(a, bins=10):
-    """Return arguments to plt.bar for pmf-like histogram of an array.
-
-    DEPRECATED: will be removed in a future version.
-
-    Parameters
-    ----------
-    a: array-like
-        array to make histogram of
-    bins: int
-        number of bins
+def _normal_quantile_func(q):
+    """
+    Compute the quantile function of the standard normal distribution.
 
-    Returns
-    -------
-    x: array
-        left x position of bars
-    h: array
-        height of bars
-    w: float
-        width of bars
+    This wrapper exists because we are dropping scipy as a mandatory dependency
+    but statistics.NormalDist was added to the standard library in 3.8.
 
     """
-    msg = "This function is deprecated and will be removed in a future version"
-    warnings.warn(msg, FutureWarning)
-    n, x = np.histogram(a, bins)
-    h = n / n.sum()
-    w = x[1] - x[0]
-    return x[:-1], h, w
+    try:
+        from statistics import NormalDist
+        qf = np.vectorize(NormalDist().inv_cdf)
+    except ImportError:
+        try:
+            from scipy.stats import norm
+            qf = norm.ppf
+        except ImportError:
+            msg = (
+                "Standard normal quantile functions require either Python>=3.8 or scipy"
+            )
+            raise RuntimeError(msg)
+    return qf(q)
 
 
 def desaturate(color, prop):
@@ -326,89 +308,12 @@ def _kde_support(data, bw, gridsize, cut, clip):
     return support
 
 
-def percentiles(a, pcts, axis=None):
-    """Like scoreatpercentile but can take and return array of percentiles.
-
-    DEPRECATED: will be removed in a future version.
-
-    Parameters
-    ----------
-    a : array
-        data
-    pcts : sequence of percentile values
-        percentile or percentiles to find score at
-    axis : int or None
-        if not None, computes scores over this axis
-
-    Returns
-    -------
-    scores: array
-        array of scores at requested percentiles
-        first dimension is length of object passed to ``pcts``
-
-    """
-    msg = "This function is deprecated and will be removed in a future version"
-    warnings.warn(msg, FutureWarning)
-
-    scores = []
-    try:
-        n = len(pcts)
-    except TypeError:
-        pcts = [pcts]
-        n = 0
-    for i, p in enumerate(pcts):
-        if axis is None:
-            score = stats.scoreatpercentile(a.ravel(), p)
-        else:
-            score = np.apply_along_axis(stats.scoreatpercentile, axis, a, p)
-        scores.append(score)
-    scores = np.asarray(scores)
-    if not n:
-        scores = scores.squeeze()
-    return scores
-
-
 def ci(a, which=95, axis=None):
     """Return a percentile range from an array of values."""
     p = 50 - which / 2, 50 + which / 2
     return np.nanpercentile(a, p, axis)
 
 
-def sig_stars(p):
-    """Return a R-style significance string corresponding to p values.
-
-    DEPRECATED: will be removed in a future version.
-
-    """
-    msg = "This function is deprecated and will be removed in a future version"
-    warnings.warn(msg, FutureWarning)
-
-    if p < 0.001:
-        return "***"
-    elif p < 0.01:
-        return "**"
-    elif p < 0.05:
-        return "*"
-    elif p < 0.1:
-        return "."
-    return ""
-
-
-def iqr(a):
-    """Calculate the IQR for an array of numbers.
-
-    DEPRECATED: will be removed in a future version.
-
-    """
-    msg = "This function is deprecated and will be removed in a future version"
-    warnings.warn(msg, FutureWarning)
-
-    a = np.asarray(a)
-    q1 = stats.scoreatpercentile(a, 25)
-    q3 = stats.scoreatpercentile(a, 75)
-    return q3 - q1
-
-
 def get_dataset_names():
     """Report available example datasets, useful for reporting issues.
 
diff --git a/setup.py b/setup.py
index ca024b45..352163a1 100644
--- a/setup.py
+++ b/setup.py
@@ -31,11 +31,17 @@ PYTHON_REQUIRES = ">=3.7"
 
 INSTALL_REQUIRES = [
     'numpy>=1.16',
-    'scipy>=1.2',
     'pandas>=0.23',
     'matplotlib>=3.0',
 ]
 
+EXTRAS_REQUIRE = {
+    'all': [
+        'scipy>=1.2',
+        'statsmodels>=0.9',
+    ]
+}
+
 
 PACKAGES = [
     'seaborn',
@@ -79,6 +85,7 @@ if __name__ == "__main__":
         download_url=DOWNLOAD_URL,
         python_requires=PYTHON_REQUIRES,
         install_requires=INSTALL_REQUIRES,
+        extras_require=EXTRAS_REQUIRE,
         packages=PACKAGES,
         classifiers=CLASSIFIERS
     )


Test patch:
diff --git a/ci/deps_latest.txt b/ci/deps_latest.txt
deleted file mode 100644
index a8d32f99..00000000
--- a/ci/deps_latest.txt
+++ /dev/null
@@ -1,5 +0,0 @@
-numpy
-scipy
-matplotlib
-pandas
-statsmodels
diff --git a/seaborn/tests/test_categorical.py b/seaborn/tests/test_categorical.py
index b2fe279e..0d514c92 100644
--- a/seaborn/tests/test_categorical.py
+++ b/seaborn/tests/test_categorical.py
@@ -1,6 +1,5 @@
 import numpy as np
 import pandas as pd
-from scipy import stats, spatial
 import matplotlib as mpl
 import matplotlib.pyplot as plt
 from matplotlib.colors import rgb2hex
@@ -11,6 +10,7 @@ import numpy.testing as npt
 
 from .. import categorical as cat
 from .. import palettes
+from ..utils import _normal_quantile_func
 
 
 class CategoricalFixture:
@@ -450,10 +450,9 @@ class TestCategoricalStatPlotter(CategoricalFixture):
                                       y.groupby(g).mean())
 
         for ci, (_, grp_y) in zip(p.confint, y.groupby(g)):
-            sem = stats.sem(grp_y)
+            sem = grp_y.std() / np.sqrt(len(grp_y))
             mean = grp_y.mean()
-            stats.norm.ppf(.975)
-            half_ci = stats.norm.ppf(.975) * sem
+            half_ci = _normal_quantile_func(.975) * sem
             ci_want = mean - half_ci, mean + half_ci
             npt.assert_array_almost_equal(ci_want, ci, 2)
 
@@ -493,9 +492,10 @@ class TestCategoricalStatPlotter(CategoricalFixture):
         assert p.statistic.shape == (4,)
         assert p.confint.shape == (4, 2)
 
-        mean = y[g == "b"].mean()
-        sem = stats.sem(y[g == "b"])
-        half_ci = stats.norm.ppf(.975) * sem
+        rows = g == "b"
+        mean = y[rows].mean()
+        sem = y[rows].std() / np.sqrt(rows.sum())
+        half_ci = _normal_quantile_func(.975) * sem
         ci = mean - half_ci, mean + half_ci
         npt.assert_almost_equal(p.statistic[1], mean)
         npt.assert_array_almost_equal(p.confint[1], ci, 2)
@@ -522,9 +522,9 @@ class TestCategoricalStatPlotter(CategoricalFixture):
 
         for ci_g, (_, grp_y) in zip(p.confint, y.groupby(g)):
             for ci, hue_y in zip(ci_g, [grp_y[::2], grp_y[1::2]]):
-                sem = stats.sem(hue_y)
+                sem = hue_y.std() / np.sqrt(len(hue_y))
                 mean = hue_y.mean()
-                half_ci = stats.norm.ppf(.975) * sem
+                half_ci = _normal_quantile_func(.975) * sem
                 ci_want = mean - half_ci, mean + half_ci
                 npt.assert_array_almost_equal(ci_want, ci, 2)
 
@@ -584,9 +584,10 @@ class TestCategoricalStatPlotter(CategoricalFixture):
         assert p.statistic.shape == (4, 3)
         assert p.confint.shape == (4, 3, 2)
 
-        mean = y[(g == "b") & (h == "x")].mean()
-        sem = stats.sem(y[(g == "b") & (h == "x")])
-        half_ci = stats.norm.ppf(.975) * sem
+        rows = (g == "b") & (h == "x")
+        mean = y[rows].mean()
+        sem = y[rows].std() / np.sqrt(rows.sum())
+        half_ci = _normal_quantile_func(.975) * sem
         ci = mean - half_ci, mean + half_ci
         npt.assert_almost_equal(p.statistic[1, 2], mean)
         npt.assert_array_almost_equal(p.confint[1, 2], ci, 2)
@@ -1131,13 +1132,11 @@ class TestViolinPlotter(CategoricalFixture):
 
         # Test reference rule bandwidth
         kde, bw = p.fit_kde(data, "scott")
-        assert isinstance(kde, stats.gaussian_kde)
         assert kde.factor == kde.scotts_factor()
         assert bw == kde.scotts_factor() * data_std
 
         # Test numeric scale factor
         kde, bw = p.fit_kde(self.y, .2)
-        assert isinstance(kde, stats.gaussian_kde)
         assert kde.factor == .2
         assert bw == .2 * data_std
 
@@ -1749,7 +1748,7 @@ class TestSwarmPlotter(CategoricalFixture):
         y = np.sort(self.y)
         orig_xy = np.c_[x, y]
         swarm = p.beeswarm(orig_xy, d)
-        dmat = spatial.distance.cdist(swarm, swarm)
+        dmat = np.sqrt(np.sum(np.square(swarm[:, np.newaxis] - swarm), axis=-1))
         triu = dmat[np.triu_indices_from(dmat, 1)]
         npt.assert_array_less(d, triu)
         npt.assert_array_equal(y, swarm[:, 1])
diff --git a/seaborn/tests/test_distributions.py b/seaborn/tests/test_distributions.py
index bf9c1a7f..e1155078 100644
--- a/seaborn/tests/test_distributions.py
+++ b/seaborn/tests/test_distributions.py
@@ -5,7 +5,6 @@ import numpy as np
 import matplotlib as mpl
 import matplotlib.pyplot as plt
 from matplotlib.colors import to_rgb, to_rgba
-from scipy import stats, integrate
 
 import pytest
 from numpy.testing import assert_array_equal, assert_array_almost_equal
@@ -21,6 +20,7 @@ from .._core import (
 from .._statistics import (
     KDE,
     Histogram,
+    _no_scipy,
 )
 from ..distributions import (
     _DistributionPlotter,
@@ -84,9 +84,17 @@ class TestDistPlot(object):
             assert len(ax.lines) == 0
             assert len(ax.collections) == 1
 
+            class Norm:
+                """Dummy object that looks like a scipy RV"""
+                def fit(self, x):
+                    return ()
+
+                def pdf(self, x, *params):
+                    return np.zeros_like(x)
+
             plt.close(ax.figure)
-            ax = distplot(self.x,
-                          hist=False, kde=False, rug=False, fit=stats.norm)
+            ax = distplot(
+                self.x, hist=False, kde=False, rug=False, fit=Norm())
             assert len(ax.patches) == 0
             assert len(ax.lines) == 1
             assert len(ax.collections) == 0
@@ -555,8 +563,9 @@ class TestKDEPlotUnivariate:
 
         ax = kdeplot(data=long_df, x="x", cut=5)
         x, y = ax.lines[0].get_xydata().T
-        assert integrate.trapz(y, x) == pytest.approx(1)
+        assert integrate(y, x) == pytest.approx(1)
 
+    @pytest.mark.skipif(_no_scipy, reason="Test requires scipy")
     def test_cumulative(self, long_df):
 
         ax = kdeplot(data=long_df, x="x", cut=5, cumulative=True)
@@ -564,6 +573,12 @@ class TestKDEPlotUnivariate:
         assert y[0] == pytest.approx(0)
         assert y[-1] == pytest.approx(1)
 
+    @pytest.mark.skipif(not _no_scipy, reason="Test requires scipy's absence")
+    def test_cumulative_requires_scipy(self, long_df):
+
+        with pytest.raises(RuntimeError):
+            kdeplot(data=long_df, x="x", cut=5, cumulative=True)
+
     def test_common_norm(self, long_df):
 
         f, (ax1, ax2) = plt.subplots(ncols=2)
@@ -578,12 +593,12 @@ class TestKDEPlotUnivariate:
         total_area = 0
         for line in ax1.lines:
             xdata, ydata = line.get_xydata().T
-            total_area += integrate.trapz(ydata, xdata)
+            total_area += integrate(ydata, xdata)
         assert total_area == pytest.approx(1)
 
         for line in ax2.lines:
             xdata, ydata = line.get_xydata().T
-            assert integrate.trapz(ydata, xdata) == pytest.approx(1)
+            assert integrate(ydata, xdata) == pytest.approx(1)
 
     def test_common_grid(self, long_df):
 
@@ -705,7 +720,7 @@ class TestKDEPlotUnivariate:
         x = rng.lognormal(0, 1, 100)
         ax = kdeplot(x=x, log_scale=True, cut=10)
         xdata, ydata = ax.lines[0].get_xydata().T
-        integral = integrate.trapz(ydata, np.log10(xdata))
+        integral = integrate(ydata, np.log10(xdata))
         assert integral == pytest.approx(1)
 
     def test_weights(self):
@@ -1294,7 +1309,7 @@ class TestHistPlotUnivariate:
         hist_area = np.multiply(bar_widths, bar_heights).sum()
 
         density, = ax.lines
-        kde_area = integrate.trapz(density.get_ydata(), density.get_xdata())
+        kde_area = integrate(density.get_ydata(), density.get_xdata())
 
         assert kde_area == pytest.approx(hist_area)
 
@@ -1317,7 +1332,7 @@ class TestHistPlotUnivariate:
             hist_area = np.multiply(bar_widths, bar_heights).sum()
 
             x, y = ax.lines[i].get_xydata().T
-            kde_area = integrate.trapz(y, x)
+            kde_area = integrate(y, x)
 
             if multiple == "layer":
                 assert kde_area == pytest.approx(hist_area)
@@ -2118,3 +2133,11 @@ class TestDisPlot:
         l1 = sum(bool(c.get_segments()) for c in g.axes.flat[0].collections)
         l2 = sum(bool(c.get_segments()) for c in g.axes.flat[1].collections)
         assert l1 == l2
+
+
+def integrate(y, x):
+    """"Simple numerical integration for testing KDE code."""
+    y = np.asarray(y)
+    x = np.asarray(x)
+    dx = np.diff(x)
+    return (dx * y[:-1] + dx * y[1:]).sum() / 2
diff --git a/seaborn/tests/test_matrix.py b/seaborn/tests/test_matrix.py
index a4adf9f5..94a7d47f 100644
--- a/seaborn/tests/test_matrix.py
+++ b/seaborn/tests/test_matrix.py
@@ -6,27 +6,31 @@ import numpy as np
 import matplotlib as mpl
 import matplotlib.pyplot as plt
 import pandas as pd
-from scipy.spatial import distance
-from scipy.cluster import hierarchy
 
-import numpy.testing as npt
 try:
-    import pandas.testing as pdt
+    from scipy.spatial import distance
+    from scipy.cluster import hierarchy
+    _no_scipy = False
 except ImportError:
-    import pandas.util.testing as pdt
-import pytest
-
-from .. import matrix as mat
-from .. import color_palette
+    _no_scipy = True
 
 try:
     import fastcluster
-
     assert fastcluster
     _no_fastcluster = False
 except ImportError:
     _no_fastcluster = True
 
+import numpy.testing as npt
+try:
+    import pandas.testing as pdt
+except ImportError:
+    import pandas.util.testing as pdt
+import pytest
+
+from .. import matrix as mat
+from .. import color_palette
+
 
 class TestHeatmap:
     rs = np.random.RandomState(sum(map(ord, "heatmap")))
@@ -459,31 +463,34 @@ class TestHeatmap:
         assert len(ax2.collections) == 2
 
 
+@pytest.mark.skipif(_no_scipy, reason="Test requires scipy")
 class TestDendrogram:
+
     rs = np.random.RandomState(sum(map(ord, "dendrogram")))
 
+    default_kws = dict(linkage=None, metric='euclidean', method='single',
+                       axis=1, label=True, rotate=False)
+
     x_norm = rs.randn(4, 8) + np.arange(8)
     x_norm = (x_norm.T + np.arange(4)).T
     letters = pd.Series(["A", "B", "C", "D", "E", "F", "G", "H"],
                         name="letters")
 
     df_norm = pd.DataFrame(x_norm, columns=letters)
-    try:
-        import fastcluster
 
-        x_norm_linkage = fastcluster.linkage_vector(x_norm.T,
-                                                    metric='euclidean',
-                                                    method='single')
-    except ImportError:
-        x_norm_distances = distance.pdist(x_norm.T, metric='euclidean')
-        x_norm_linkage = hierarchy.linkage(x_norm_distances, method='single')
-    x_norm_dendrogram = hierarchy.dendrogram(x_norm_linkage, no_plot=True,
-                                             color_threshold=-np.inf)
-    x_norm_leaves = x_norm_dendrogram['leaves']
-    df_norm_leaves = np.asarray(df_norm.columns[x_norm_leaves])
+    if not _no_scipy:
+        if _no_fastcluster:
+            x_norm_distances = distance.pdist(x_norm.T, metric='euclidean')
+            x_norm_linkage = hierarchy.linkage(x_norm_distances, method='single')
+        else:
+            x_norm_linkage = fastcluster.linkage_vector(x_norm.T,
+                                                        metric='euclidean',
+                                                        method='single')
 
-    default_kws = dict(linkage=None, metric='euclidean', method='single',
-                       axis=1, label=True, rotate=False)
+        x_norm_dendrogram = hierarchy.dendrogram(x_norm_linkage, no_plot=True,
+                                                 color_threshold=-np.inf)
+        x_norm_leaves = x_norm_dendrogram['leaves']
+        df_norm_leaves = np.asarray(df_norm.columns[x_norm_leaves])
 
     def test_ndarray_input(self):
         p = mat._DendrogramPlotter(self.x_norm, **self.default_kws)
@@ -701,7 +708,9 @@ class TestDendrogram:
         plt.close(f)
 
 
+@pytest.mark.skipif(_no_scipy, reason="Test requires scipy")
 class TestClustermap:
+
     rs = np.random.RandomState(sum(map(ord, "clustermap")))
 
     x_norm = rs.randn(4, 8) + np.arange(8)
@@ -710,19 +719,6 @@ class TestClustermap:
                         name="letters")
 
     df_norm = pd.DataFrame(x_norm, columns=letters)
-    try:
-        import fastcluster
-
-        x_norm_linkage = fastcluster.linkage_vector(x_norm.T,
-                                                    metric='euclidean',
-                                                    method='single')
-    except ImportError:
-        x_norm_distances = distance.pdist(x_norm.T, metric='euclidean')
-        x_norm_linkage = hierarchy.linkage(x_norm_distances, method='single')
-    x_norm_dendrogram = hierarchy.dendrogram(x_norm_linkage, no_plot=True,
-                                             color_threshold=-np.inf)
-    x_norm_leaves = x_norm_dendrogram['leaves']
-    df_norm_leaves = np.asarray(df_norm.columns[x_norm_leaves])
 
     default_kws = dict(pivot_kws=None, z_score=None, standard_scale=None,
                        figsize=(10, 10), row_colors=None, col_colors=None,
@@ -738,6 +734,20 @@ class TestClustermap:
     row_colors = color_palette('Set2', df_norm.shape[0])
     col_colors = color_palette('Dark2', df_norm.shape[1])
 
+    if not _no_scipy:
+        if _no_fastcluster:
+            x_norm_distances = distance.pdist(x_norm.T, metric='euclidean')
+            x_norm_linkage = hierarchy.linkage(x_norm_distances, method='single')
+        else:
+            x_norm_linkage = fastcluster.linkage_vector(x_norm.T,
+                                                        metric='euclidean',
+                                                        method='single')
+
+        x_norm_dendrogram = hierarchy.dendrogram(x_norm_linkage, no_plot=True,
+                                                 color_threshold=-np.inf)
+        x_norm_leaves = x_norm_dendrogram['leaves']
+        df_norm_leaves = np.asarray(df_norm.columns[x_norm_leaves])
+
     def test_ndarray_input(self):
         cg = mat.ClusterGrid(self.x_norm, **self.default_kws)
         pdt.assert_frame_equal(cg.data, pd.DataFrame(self.x_norm))
@@ -1320,3 +1330,19 @@ class TestClustermap:
         for ax in [g.ax_col_dendrogram, g.ax_row_dendrogram]:
             tree, = ax.collections
             assert tuple(tree.get_color().squeeze())[:3] == rgb
+
+
+if _no_scipy:
+
+    def test_required_scipy_errors():
+
+        x = np.random.normal(0, 1, (10, 10))
+
+        with pytest.raises(RuntimeError):
+            mat.clustermap(x)
+
+        with pytest.raises(RuntimeError):
+            mat.ClusterGrid(x)
+
+        with pytest.raises(RuntimeError):
+            mat.dendrogram(x)
diff --git a/seaborn/tests/test_statistics.py b/seaborn/tests/test_statistics.py
index 4590849e..380f53e5 100644
--- a/seaborn/tests/test_statistics.py
+++ b/seaborn/tests/test_statistics.py
@@ -1,5 +1,4 @@
 import numpy as np
-from scipy import integrate
 
 try:
     import statsmodels.distributions as smdist
@@ -13,6 +12,7 @@ from .._statistics import (
     KDE,
     Histogram,
     ECDF,
+    _no_scipy,
 )
 
 
@@ -33,6 +33,12 @@ class DistributionFixtures:
 
 class TestKDE:
 
+    def integrate(self, y, x):
+        y = np.asarray(y)
+        x = np.asarray(x)
+        dx = np.diff(x)
+        return (dx * y[:-1] + dx * y[1:]).sum() / 2
+
     def test_gridsize(self, rng):
 
         x = rng.normal(0, 3, 1000)
@@ -75,8 +81,9 @@ class TestKDE:
         x = rng.normal(0, 3, 1000)
         kde = KDE()
         density, support = kde(x)
-        assert integrate.trapz(density, support) == pytest.approx(1, abs=1e-5)
+        assert self.integrate(density, support) == pytest.approx(1, abs=1e-5)
 
+    @pytest.mark.skipif(_no_scipy, reason="Test requires scipy")
     def test_cumulative(self, rng):
 
         x = rng.normal(0, 3, 1000)
@@ -138,6 +145,7 @@ class TestKDE:
         total = density.sum() * (dx * dy)
         assert total == pytest.approx(1, abs=1e-2)
 
+    @pytest.mark.skipif(_no_scipy, reason="Test requires scipy")
     def test_bivariate_cumulative(self, rng):
 
         x, y = rng.normal(0, 3, (2, 50))
diff --git a/seaborn/tests/test_utils.py b/seaborn/tests/test_utils.py
index e01bc46a..8e14d029 100644
--- a/seaborn/tests/test_utils.py
+++ b/seaborn/tests/test_utils.py
@@ -58,37 +58,6 @@ def _network(t=None, url="https://github.com"):
     return wrapper
 
 
-def test_pmf_hist_basics():
-    """Test the function to return barplot args for pmf hist."""
-    with pytest.warns(FutureWarning):
-        out = utils.pmf_hist(a_norm)
-    assert len(out) == 3
-    x, h, w = out
-    assert len(x) == len(h)
-
-    # Test simple case
-    a = np.arange(10)
-    with pytest.warns(FutureWarning):
-        x, h, w = utils.pmf_hist(a, 10)
-    assert np.all(h == h[0])
-
-    # Test width
-    with pytest.warns(FutureWarning):
-        x, h, w = utils.pmf_hist(a_norm)
-    assert x[1] - x[0] == w
-
-    # Test normalization
-    with pytest.warns(FutureWarning):
-        x, h, w = utils.pmf_hist(a_norm)
-    assert sum(h) == pytest.approx(1)
-    assert h.max() <= 1
-
-    # Test bins
-    with pytest.warns(FutureWarning):
-        x, h, w = utils.pmf_hist(a_norm, 20)
-    assert len(x) == 20
-
-
 def test_ci_to_errsize():
     """Test behavior of ci_to_errsize."""
     cis = [[.5, .5],
@@ -130,24 +99,6 @@ def test_saturate():
     assert out == (1, 0, 0)
 
 
-@pytest.mark.parametrize(
-    "p,annot", [(.0001, "***"), (.001, "**"), (.01, "*"), (.09, "."), (1, "")]
-)
-def test_sig_stars(p, annot):
-    """Test the sig stars function."""
-    with pytest.warns(FutureWarning):
-        stars = utils.sig_stars(p)
-        assert stars == annot
-
-
-def test_iqr():
-    """Test the IQR function."""
-    a = np.arange(5)
-    with pytest.warns(FutureWarning):
-        iqr = utils.iqr(a)
-    assert iqr == 2
-
-
 @pytest.mark.parametrize(
     "s,exp",
     [