Repository: alteryx/woodwork
Pull request #1547 (2022-10-24 06:19:41)

Added medcouple as a method for outlier detection

codecov[bot]: # [Codecov](https://codecov.io/gh/alteryx/woodwork/pull/1547?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alteryx) Report
> Merging [#1547](https://codecov.io/gh/alteryx/woodwork/pull/1547?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alteryx) (3f09b39) into [main](https://codecov.io/gh/alteryx/woodwork/commit/2e8e5e9905db4f4433344242d67da935c0af3b00?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alteryx) (2e8e5e9) will **decrease** coverage by `81.63%`.
> The diff coverage is `17.74%`.

```diff
@@             Coverage Diff             @@
##             main    #1547       +/-   ##
===========================================
- Coverage   98.91%   17.28%   -81.64%     
===========================================
  Files          91       92        +1     
  Lines       10813    10933      +120     
===========================================
- Hits        10696     1890     -8806     
- Misses        117     9043     +8926     
```


| [Impacted Files](https://codecov.io/gh/alteryx/woodwork/pull/1547?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alteryx) | Coverage Δ | |
|---|---|---|
| [woodwork/tests/accessor/test\_statistics.py](https://codecov.io/gh/alteryx/woodwork/pull/1547/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alteryx#diff-d29vZHdvcmsvdGVzdHMvYWNjZXNzb3IvdGVzdF9zdGF0aXN0aWNzLnB5) | `0.00% <0.00%> (-100.00%)` | :arrow_down: |
| [.../statistics\_utils/\_get\_box\_plot\_info\_for\_column.py](https://codecov.io/gh/alteryx/woodwork/pull/1547/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alteryx#diff-d29vZHdvcmsvc3RhdGlzdGljc191dGlscy9fZ2V0X2JveF9wbG90X2luZm9fZm9yX2NvbHVtbi5weQ==) | `15.51% <13.04%> (-84.49%)` | :arrow_down: |
| [...statistics\_utils/\_determine\_best\_outlier\_method.py](https://codecov.io/gh/alteryx/woodwork/pull/1547/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alteryx#diff-d29vZHdvcmsvc3RhdGlzdGljc191dGlscy9fZGV0ZXJtaW5lX2Jlc3Rfb3V0bGllcl9tZXRob2QucHk=) | `33.33% <33.33%> (ø)` | |
| [woodwork/tests/conftest.py](https://codecov.io/gh/alteryx/woodwork/pull/1547/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alteryx#diff-d29vZHdvcmsvdGVzdHMvY29uZnRlc3QucHk=) | `57.07% <50.00%> (-42.93%)` | :arrow_down: |
| [woodwork/column\_accessor.py](https://codecov.io/gh/alteryx/woodwork/pull/1547/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alteryx#diff-d29vZHdvcmsvY29sdW1uX2FjY2Vzc29yLnB5) | `51.63% <71.42%> (-48.37%)` | :arrow_down: |
| [woodwork/tests/test\_config.py](https://codecov.io/gh/alteryx/woodwork/pull/1547/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alteryx#diff-d29vZHdvcmsvdGVzdHMvdGVzdF9jb25maWcucHk=) | `0.00% <0.00%> (-100.00%)` | :arrow_down: |
| [woodwork/tests/test\_version.py](https://codecov.io/gh/alteryx/woodwork/pull/1547/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alteryx#diff-d29vZHdvcmsvdGVzdHMvdGVzdF92ZXJzaW9uLnB5) | `0.00% <0.00%> (-100.00%)` | :arrow_down: |
| [woodwork/tests/utils/test\_utils.py](https://codecov.io/gh/alteryx/woodwork/pull/1547/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alteryx#diff-d29vZHdvcmsvdGVzdHMvdXRpbHMvdGVzdF91dGlscy5weQ==) | `0.00% <0.00%> (-100.00%)` | :arrow_down: |
| [woodwork/tests/utils/test\_concat.py](https://codecov.io/gh/alteryx/woodwork/pull/1547/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alteryx#diff-d29vZHdvcmsvdGVzdHMvdXRpbHMvdGVzdF9jb25jYXQucHk=) | `0.00% <0.00%> (-100.00%)` | :arrow_down: |
| [woodwork/tests/type\_system/conftest.py](https://codecov.io/gh/alteryx/woodwork/pull/1547/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alteryx#diff-d29vZHdvcmsvdGVzdHMvdHlwZV9zeXN0ZW0vY29uZnRlc3QucHk=) | `0.00% <0.00%> (-100.00%)` | :arrow_down: |
| ... and [60 more](https://codecov.io/gh/alteryx/woodwork/pull/1547/diff?src=pr&el=tree-more&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alteryx) | |

Help us with your feedback. Take ten seconds to tell us [how you rate us](https://about.codecov.io/nps?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alteryx). Have a feature suggestion? [Share it here.](https://app.codecov.io/gh/feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=alteryx)

ParthivNaresh: Here are the perf [tests](https://alteryx.atlassian.net/wiki/spaces/PS/pages/1576730738/Box+Plot+vs+Medcouple)
Cmancuso: Some small comments, but looking good to me!

Patch:
diff --git a/docs/source/guides/statistical_insights.ipynb b/docs/source/guides/statistical_insights.ipynb
index 528a87f..42dc87d 100644
--- a/docs/source/guides/statistical_insights.ipynb
+++ b/docs/source/guides/statistical_insights.ipynb
@@ -178,7 +178,7 @@
    "source": [
     "## Outlier Detection with Series.ww.box_plot_dict\n",
     "\n",
-    "Woodwork allows for univariate outlier detection using the IQR, or interquartile range, method. This can be done on a by-column basis using the `series.ww.box_plot_dict` method that identifies outliers and includes the statistical data necessary for building a box and whisker plot."
+    "One of the ways that Woodwork allows for univariate outlier detection is the IQR, or interquartile range, method. This can be done on a by-column basis using the `series.ww.box_plot_dict` method that identifies outliers and includes the statistical data necessary for building a box and whisker plot."
    ]
   },
   {
@@ -214,6 +214,8 @@
     "s = pd.Series(rnd.normal(50, 10, 401604))\n",
     "s.ww.init()\n",
     "box_plot_dict = s.ww.box_plot_dict()\n",
+    "\n",
+    "print(\"high bound: \", box_plot_dict[\"method\"])\n",
     "print(\"high bound: \", box_plot_dict[\"high_bound\"])\n",
     "print(\"low_bound: \", box_plot_dict[\"low_bound\"])\n",
     "print(\"quantiles: \", box_plot_dict[\"quantiles\"])\n",
@@ -228,6 +230,79 @@
     "With the normally distributed set of data, the likelyhood of outliers is closer to what we'd expect, around .7%. "
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Outlier Detection with Series.ww.medcouple_dict\n",
+    "\n",
+    "Another way that Woodwork allows for univariate outlier detection is through the medcouple statistic, implemented via `series.ww.medcouple_dict`. Similar to `series.ww.box_plot_dict`, this method also returns the information necessary to build a box and whiskers plot"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "total = df.ww[\"total\"]\n",
+    "medcouple_dict = total.ww.medcouple_dict()\n",
+    "\n",
+    "print(\"medcouple: \", medcouple_dict[\"method\"])\n",
+    "print(\"medcouple: \", medcouple_dict[\"medcouple_stat\"])\n",
+    "print(\"high bound: \", medcouple_dict[\"high_bound\"])\n",
+    "print(\"low_bound: \", medcouple_dict[\"low_bound\"])\n",
+    "print(\"quantiles: \", medcouple_dict[\"quantiles\"])\n",
+    "print(\"number of low outliers: \", len(medcouple_dict[\"low_values\"]))\n",
+    "print(\"number of high outliers: \", len(medcouple_dict[\"high_values\"]))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The number of outliers identified has decreased to around 7% of the rows. Additionally the direction of these outliers has changed, with more lower outliers and fewer higher outliers being identified. This is because the medcouple statistic takes into account the skewness of the distribution when determining which points are most likely to be outliers. For distributions that have a long tail, the Box Plot method might not be as appropriate.\n",
+    "The medcouple statistic is calculated based on at most 10,000 randomly sampled observations in the column. This value can be changed within `Config` settings via `ww.config.set_option(\"medcouple_sample_size\", N)`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "rnd = np.random.RandomState(33)\n",
+    "s = pd.Series(rnd.normal(50, 10, 401604))\n",
+    "s.ww.init()\n",
+    "medcouple_dict = s.ww.medcouple_dict()\n",
+    "print(\"medcouple: \", medcouple_dict[\"medcouple_stat\"])\n",
+    "print(\"high bound: \", medcouple_dict[\"high_bound\"])\n",
+    "print(\"low_bound: \", medcouple_dict[\"low_bound\"])\n",
+    "print(\"quantiles: \", medcouple_dict[\"quantiles\"])\n",
+    "print(\"number of low outliers: \", len(medcouple_dict[\"low_values\"]))\n",
+    "print(\"number of high outliers: \", len(medcouple_dict[\"high_values\"]))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "As you can see, the medcouple statistic is close to 0, which indicates a data with little to no skew. The outliers identified via the Medcouple statistic are around the same quantity as those via the Box Plot method when the data is normal."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Outlier Detection with Series.ww.get_outliers\n",
+    "\n",
+    "Knowing whether to use the box plot method or the medcouple approach can be confusing. You might not know whether a distrbiution is skewed or not, and if it is, whether its skewed enough to warrant using `medcouple_dict`. The method `get_outliers` can be used to address these issues. Running `series.ww.get_outliers` with the default method of `best` will return the outlier information based on the best approach.\n",
+    "\n",
+    "This is determined by taking the medcouple statistic of a random sampling of the series. If the absolute value of the medcouple statistic of that sampling is under the default value of 0.3, then the method chosen will be box plot. If it is greater than or equal to 0.3, it will be the medcouple approach that is used (as this indicates that there is at least an average amount of skewness in the distribution).\n",
+    "\n",
+    "To change this default medcouple threshold of 0.3, feel free to change the value in `Config` settings via `ww.config.set_option(\"medcouple_threshold\", 0.5)`."
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
diff --git a/docs/source/release_notes.rst b/docs/source/release_notes.rst
index b8d7fcc..99a1e5d 100644
--- a/docs/source/release_notes.rst
+++ b/docs/source/release_notes.rst
@@ -7,6 +7,7 @@ Future Release
 ==============
     * Enhancements
         * Improved ``Boolean`` and ``BooleanNullable`` inference to detect common string representations of boolean values (:pr:`1549`)
+        * Added the ``get_outliers`` and ``medcouple_dict`` functions to ``WoodworkColumnAccessor`` so that the medcouple statistic can be used for outlier detection (:pr:`1547`)
     * Fixes
         * Resolved FutureWarning in ``_get_box_plot_info_for_column`` (:pr:`1563`)
         * Fixed error message in validate method in logical_types.py (:pr:`1565`)
@@ -18,10 +19,11 @@ Future Release
         * Unpin dask dependency (:pr:`1561`)
         * Changed the sampling strategy for type inference from ``head`` to random (:pr:`1566`)
     * Documentation Changes
+        * Updated documentation to include the ``get_outliers`` and ``medcouple_dict`` (:pr:`1547`)
     * Testing Changes
 
     Thanks to the following people for contributing to this release:
-    :user:`gsheni`, :user:`sbadithe`, :user:`bchen1116`
+    :user:`bchen1116`, :user:`gsheni`, :user:`ParthivNaresh`, :user:`sbadithe`
 
 Breaking Changes
 ++++++++++++++++
diff --git a/woodwork/column_accessor.py b/woodwork/column_accessor.py
index 027f58e..7084b54 100644
--- a/woodwork/column_accessor.py
+++ b/woodwork/column_accessor.py
@@ -1,6 +1,7 @@
 import copy
 import warnings
 import weakref
+from typing import Dict, Optional, Union
 
 import pandas as pd
 from pandas.api import types as pdtypes
@@ -388,12 +389,64 @@ class WoodworkColumnAccessor:
         """
         self._schema._set_semantic_tags(semantic_tags)
 
+    @_check_column_schema
+    def get_outliers(
+        self,
+        method="best",
+        quantiles: Optional[Dict[float, Union[int, float]]] = None,
+        include_indices_and_values: bool = True,
+        ignore_zeros: bool = False,
+    ):
+        """Gets the information necessary to create a box and whisker plot with outliers for a numeric column
+        using the selected method.
+
+        Args:
+            method (str): The method to use when calculating the box and whiskers plot. Options are 'best', 'box_plot' and 'medcouple'.
+                Defaults to 'best' at which point a heuristic will determine the appropriate method to use.
+            quantiles (dict[float -> int or float], optional): A dictionary containing the quantiles for the data
+                where the key indicates the quantile, and the value is the quantile's value for the data. If
+                no quantiles are provided, they will be computed from the data.
+            include_indices_and_values (bool, optional): Whether or not the lists containing individual
+                outlier values and their indices will be included in the returned dictionary.
+                Defaults to True.
+            ignore_zeros (bool): Whether to ignore 0 values (not NaN values) when calculating the box plot and outliers.
+                Defaults to False.
+
+        Returns:
+            (dict[str -> float,list[number]]): Returns a dictionary containing box plot information for the Series.
+                The following elements will be found in the dictionary:
+
+                - low_bound (float): the lowest data point in the dataset excluding any outliers - to be used as a whisker
+                - high_bound (float): the highest point in the dataset excluding any outliers - to be used as a whisker
+                - quantiles (list[float]): the quantiles used to determine the bounds.
+                    If quantiles were passed in, will contain all quantiles passed in. Otherwise, contains the five
+                    quantiles {0.0, 0.25, 0.5, 0.75, 1.0}.
+                - low_values (list[float, int], optional): the values of the lower outliers.
+                    Will not be included if ``include_indices_and_values`` is False.
+                - high_values (list[float, int], optional): the values of the upper outliers
+                    Will not be included if ``include_indices_and_values`` is False.
+                - low_indices (list[int], optional): the corresponding index values for each of the lower outliers
+                    Will not be included if ``include_indices_and_values`` is False.
+                - high_indices (list[int], optional): the corresponding index values for each of the upper outliers
+                    Will not be included if ``include_indices_and_values`` is False.
+                - method (str): the method used to identify outliers
+                - medcouple_stat (float): the medcouple statistic will be added in the response if the method selected
+                    is medcouple
+        """
+        return _get_box_plot_info_for_column(
+            self._series,
+            method=method,
+            quantiles=quantiles,
+            include_indices_and_values=include_indices_and_values,
+            ignore_zeros=ignore_zeros,
+        )
+
     @_check_column_schema
     def box_plot_dict(
         self,
-        quantiles=None,
-        include_indices_and_values=True,
-        ignore_zeros=False,
+        quantiles: Optional[Dict[int, int]] = None,
+        include_indices_and_values: bool = True,
+        ignore_zeros: bool = False,
     ):
         """Gets the information necessary to create a box and whisker plot with outliers for a numeric column
         using the IQR method.
@@ -401,7 +454,7 @@ class WoodworkColumnAccessor:
         Args:
             quantiles (dict[float -> float], optional): A dictionary containing the quantiles for the data
                 where the key indicates the quantile, and the value is the quantile's value for the data. If
-                no qantiles are provided, they will be computed from the data.
+                no quantiles are provided, they will be computed from the data.
             include_indices_and_values (bool, optional): Whether or not the lists containing individual
                 outlier values and their indices will be included in the returned dictionary.
                 Defaults to True.
@@ -431,9 +484,65 @@ class WoodworkColumnAccessor:
                     Will not be included if ``include_indices_and_values`` is False.
                 - high_indices (list[int], optional): the corresponding index values for each of the upper outliers
                     Will not be included if ``include_indices_and_values`` is False.
+                - method (str): the method used to identify outliers, in this case box_plot
+        """
+        return _get_box_plot_info_for_column(
+            self._series,
+            method="box_plot",
+            quantiles=quantiles,
+            include_indices_and_values=include_indices_and_values,
+            ignore_zeros=ignore_zeros,
+        )
+
+    @_check_column_schema
+    def medcouple_dict(
+        self,
+        quantiles: Optional[Dict[int, int]] = None,
+        include_indices_and_values: bool = True,
+        ignore_zeros: bool = False,
+    ):
+        """Gets the information necessary to create a box and whisker plot with outliers for a numeric column
+        using the Medcouple statistic method.
+
+        Args:
+            quantiles (dict[float -> float], optional): A dictionary containing the quantiles for the data
+                where the key indicates the quantile, and the value is the quantile's value for the data. If
+                no quantiles are provided, they will be computed from the data.
+            include_indices_and_values (bool, optional): Whether or not the lists containing individual
+                outlier values and their indices will be included in the returned dictionary.
+                Defaults to True.
+            ignore_zeros (bool): Whether to ignore 0 values (not NaN values) when calculating the box plot and outliers.
+                Defaults to False.
+
+        Note:
+            The minimum quantiles necessary for building a box plot using the Medcouple statistic method are the
+            minimum value (0.0 in the quantiles dict), first quartile (0.25), third quartile (0.75), and maximum value (1.0).
+            If no quantiles are provided, the following quantiles will be calculated:
+            {0.0, 0.25, 0.5, 0.75, 1.0}, which correspond to {min, first quantile, median, third quantile, max}.
+
+        Returns:
+            (dict[str -> float,list[number]]): Returns a dictionary containing box plot information for the Series.
+                The following elements will be found in the dictionary:
+
+                - low_bound (float): the lowest data point in the dataset excluding any outliers - to be used as a whisker
+                - high_bound (float): the highest point in the dataset excluding any outliers - to be used as a whisker
+                - quantiles (list[float]): the quantiles used to determine the bounds.
+                    If quantiles were passed in, will contain all quantiles passed in. Otherwise, contains the five
+                    quantiles {0.0, 0.25, 0.5, 0.75, 1.0}.
+                - low_values (list[float, int], optional): the values of the lower outliers.
+                    Will not be included if ``include_indices_and_values`` is False.
+                - high_values (list[float, int], optional): the values of the upper outliers
+                    Will not be included if ``include_indices_and_values`` is False.
+                - low_indices (list[int], optional): the corresponding index values for each of the lower outliers
+                    Will not be included if ``include_indices_and_values`` is False.
+                - high_indices (list[int], optional): the corresponding index values for each of the upper outliers
+                    Will not be included if ``include_indices_and_values`` is False.
+                - method (str): the method used to identify outliers, in this case medcouple.
+                - medcouple_stat (float): the medcouple statistic
         """
         return _get_box_plot_info_for_column(
             self._series,
+            method="medcouple",
             quantiles=quantiles,
             include_indices_and_values=include_indices_and_values,
             ignore_zeros=ignore_zeros,
diff --git a/woodwork/config.py b/woodwork/config.py
index 236709d..3583c8c 100644
--- a/woodwork/config.py
+++ b/woodwork/config.py
@@ -51,6 +51,8 @@ CONFIG_DEFAULTS = {
     "frequence_inference_window_length": 15,
     "frequence_inference_threshold": 0.9,
     "correlation_metrics": ["mutual_info", "pearson", "spearman", "max", "all"],
+    "medcouple_threshold": 0.3,  # Must be between 0.0 and 1.0
+    "medcouple_sample_size": 10000,
     "boolean_inference_strings": [
         ["yes", "no"],
         ["y", "n"],
diff --git a/woodwork/statistics_utils/__init__.py b/woodwork/statistics_utils/__init__.py
index 23e7f56..db46389 100644
--- a/woodwork/statistics_utils/__init__.py
+++ b/woodwork/statistics_utils/__init__.py
@@ -4,12 +4,20 @@ from woodwork.statistics_utils._bin_numeric_cols_into_categories import (
 from woodwork.statistics_utils._convert_ordinal_to_numeric import (
     _convert_ordinal_to_numeric,
 )
+from woodwork.statistics_utils._get_medcouple_statistic import (
+    _get_medcouple_statistic,
+    _sample_for_medcouple,
+)
 from woodwork.statistics_utils._get_box_plot_info_for_column import (
     _get_box_plot_info_for_column,
+    _get_low_high_bound,
 )
 from woodwork.statistics_utils._get_dependence_dict import _get_dependence_dict
 from woodwork.statistics_utils._get_describe_dict import _get_describe_dict
 from woodwork.statistics_utils._get_histogram_values import _get_histogram_values
+from woodwork.statistics_utils._medcouple_implementation import (
+    _calculate_medcouple_statistic,
+)
 from woodwork.statistics_utils._get_mode import _get_mode
 from woodwork.statistics_utils._get_numeric_value_counts_in_range import (
     _get_numeric_value_counts_in_range,
diff --git a/woodwork/statistics_utils/_get_box_plot_info_for_column.py b/woodwork/statistics_utils/_get_box_plot_info_for_column.py
index a5883a4..21292f1 100644
--- a/woodwork/statistics_utils/_get_box_plot_info_for_column.py
+++ b/woodwork/statistics_utils/_get_box_plot_info_for_column.py
@@ -1,24 +1,110 @@
+from collections import namedtuple
+
 import numpy as np
 import pandas as pd
+from scipy.stats import skew
 
+import woodwork as ww
+from woodwork.statistics_utils._get_medcouple_statistic import _sample_for_medcouple
 from woodwork.utils import import_or_none
 
 dd = import_or_none("dask.dataframe")
 ps = import_or_none("pyspark.pandas")
 
 
+method_result = namedtuple(
+    "MedcoupleHeuristicResult",
+    (
+        "method",
+        "mc",
+    ),
+)
+
+
+def _determine_best_outlier_method(series):
+    """Determines the best outlier method to use based on the distribution of the series and outcome of the medcouple statistic.
+
+    Args:
+        series (Series): Data on which the medcouple statistic will be run in order to determine skewness.
+
+    Note:
+        The calculation of the medcouple statistic has a large memory requirement of O(N**2), therefore larger series will
+        have a random subset selected in order to determine skewness and the best outlier method.
+
+    Returns:
+        MedcoupleHeuristicResult - Named tuple with 2 fields
+            method (str): Name of the outlier method to use.
+            mc (float): The medcouple statistic (if the method chosen is medcouple, otherwise None).
+    """
+    mc = _sample_for_medcouple(series)
+    method = "medcouple"
+    if np.abs(mc) < ww.config.get_option("medcouple_threshold"):
+        method = "box_plot"
+        mc = None
+    return method_result(method, mc)
+
+
+def _determine_coefficients(series, mc):
+    """Determines the coefficients to use in the medcouple bounds equation using the moment of skewness.
+
+    Args:
+        series (Series): Data on which the medcouple statistic will be run in order to determine coefficients.
+
+    Note:
+        Based off of the paper I. H. Adil and A. R. Irshad, "A Modified Approach for Detection of Outliers",
+        Pak.j.stat.oper.res. Vol.XI No.1 2015 pp. 91-102, 2015.
+
+    Returns:
+        MedcoupleHeuristicResult - Named tuple with 2 fields
+            method (str): Name of the outlier method to use.
+            mc (float): The medcouple statistic (if the method chosen is medcouple, otherwise None).
+    """
+    coeff = np.abs(skew(series))
+    coeff = min(coeff, 3.5)
+    if mc >= 0:
+        return -coeff, coeff
+    else:
+        return coeff, -coeff
+
+
+def _get_low_high_bound(series, method, q1, q3, min_value, max_value, mc=None):
+    iqr = q3 - q1
+    if method == "medcouple":
+        if mc is None:
+            raise ValueError(
+                "If the method selected is medcouple, then mc cannot be None."
+            )
+        # Medcouple coefficients change based on the skew direction
+        lower_bound_coef, higher_bound_coef = _determine_coefficients(series, mc)
+        low_bound = q1 - 1.5 * np.exp(lower_bound_coef * mc) * iqr
+        high_bound = q3 + 1.5 * np.exp(higher_bound_coef * mc) * iqr
+    elif method == "box_plot":
+        low_bound = q1 - (iqr * 1.5)
+        high_bound = q3 + (iqr * 1.5)
+    else:
+        raise ValueError(
+            f"Acceptable methods are 'box_plot' and 'medcouple'. The value passed was '{method}'."
+        )
+    # Box plot bounds calculation - the bounds should never be beyond the min and max values
+    low_bound = max(low_bound, min_value)
+    high_bound = min(high_bound, max_value)
+    return low_bound, high_bound
+
+
 def _get_box_plot_info_for_column(
     series,
+    method="best",
     quantiles=None,
     include_indices_and_values=True,
     ignore_zeros=False,
 ):
-    """Gets the information necessary to create a box and whisker plot with outliers for a numeric column
-        using the IQR method.
+    """Gets the information necessary to create a box and whisker plot with outliers for a numeric column.
 
     Args:
         series (Series): Data for which the box plot and outlier information will be gathered.
             Will be used to calculate quantiles if none are provided.
+        method (str): The method to use when calculating the box and whiskers plot. Options are 'box_plot' and 'medcouple'.
+        Defaults to 'best' at which point a heuristic will determine the appropriate method to use.
         quantiles (dict[float -> float], optional): A dictionary containing the quantiles for the data
             where the key indicates the quantile, and the value is the quantile's value for the data.
         include_indices_and_values (bool, optional): Whether or not the lists containing individual
@@ -104,10 +190,20 @@ def _get_box_plot_info_for_column(
     q3 = quantiles[0.75]
     max_value = quantiles[1.0]
 
-    # Box plot bounds calculation - the bounds should never be beyond the min and max values
-    iqr = q3 - q1
-    low_bound = max(q1 - (iqr * 1.5), min_value)
-    high_bound = min(q3 + (iqr * 1.5), max_value)
+    mc = None
+    if method == "best":
+        method, mc = _determine_best_outlier_method(series)
+    if method == "medcouple" and mc is None:
+        mc = _sample_for_medcouple(series)
+
+    if isinstance(mc, np.ndarray):
+        mc = mc.tolist()
+
+    low_bound, high_bound = _get_low_high_bound(
+        series, method, q1, q3, min_value, max_value, mc
+    )
+    low_bound = round(low_bound, 5) if isinstance(low_bound, float) else low_bound
+    high_bound = round(high_bound, 5) if isinstance(high_bound, float) else high_bound
 
     if include_indices_and_values:
         # identify outliers in the series
@@ -129,6 +225,10 @@ def _get_box_plot_info_for_column(
             "high_indices": high_series.index.tolist(),
         }
 
+    outliers_dict["method"] = method
+    if method == "medcouple":
+        outliers_dict["medcouple_stat"] = mc
+
     return {
         "low_bound": low_bound,
         "high_bound": high_bound,
diff --git a/woodwork/statistics_utils/_get_medcouple_statistic.py b/woodwork/statistics_utils/_get_medcouple_statistic.py
new file mode 100644
index 0000000..3da602c
--- /dev/null
+++ b/woodwork/statistics_utils/_get_medcouple_statistic.py
@@ -0,0 +1,43 @@
+import numpy as np
+import pandas as pd
+
+import woodwork as ww
+from woodwork.statistics_utils._medcouple_implementation import (
+    _calculate_medcouple_statistic,
+)
+
+
+def _get_medcouple_statistic(array_, axis=0):
+    array_ = np.asarray(array_, dtype=np.double)
+    mc = np.apply_along_axis(_calculate_medcouple_statistic, axis, array_)
+    if isinstance(mc, np.ndarray) and isinstance(mc.tolist(), float):
+        mc = mc.tolist()
+        mc = round(mc, 3)
+    return mc
+
+
+def _sample_for_medcouple(series, seed=42):
+    """Samples the data in order to calculate the medcouple statistic.
+
+    Args:
+        series (Series): Data on which sampling will occur.
+        seed (int): Seed for random sampling.
+
+    Note:
+        The calculation of the medcouple statistic has a large memory requirement of O(N**2), therefore a series over 10,000
+        observations will be randomly sampled up to 10,000.
+
+    Returns:
+        float: Medcouple statistic from the sampled data
+    """
+    np.random.seed(seed)
+    series_size = len(series)
+    if series_size < ww.config.get_option("medcouple_sample_size"):
+        series_sample = series.copy()
+    else:
+        series_sample = np.random.choice(
+            series, ww.config.get_option("medcouple_sample_size"), replace=False
+        )
+    col = pd.Series(series_sample)
+    mc = _get_medcouple_statistic(col)
+    return mc
diff --git a/woodwork/statistics_utils/_medcouple_implementation.py b/woodwork/statistics_utils/_medcouple_implementation.py
new file mode 100644
index 0000000..d025ea5
--- /dev/null
+++ b/woodwork/statistics_utils/_medcouple_implementation.py
@@ -0,0 +1,39 @@
+import numpy as np
+
+
+def _calculate_medcouple_statistic(y):
+    """Calculates the medcouple statistic. Based on the implementation by statsmodels and off of the paper
+    M. Hubert and E. Vandervieren, "An adjusted boxplot for skewed distributions" Computational Statistics
+    & Data Analysis, vol. 52, pp. 5186-5201, August 2008.
+
+    Args:
+        y (np.ndarray, pd.Series): Data on which sampling will occur.
+
+    Returns:
+        float: Medcouple statistic
+    """
+    y = np.squeeze(np.asarray(y))
+    y = np.sort(y)
+
+    y_length = y.shape[0]
+    if y_length % 2 == 0:
+        median = (y[y_length // 2 - 1] + y[y_length // 2]) / 2
+    else:
+        median = y[(y_length - 1) // 2]
+
+    deviations = y - median
+    lower = deviations[deviations <= 0.0]
+    upper = deviations[deviations >= 0.0]
+    upper = upper[:, None]
+    standardization = upper - lower
+    is_zero = np.logical_and(lower == 0.0, upper == 0.0)
+    standardization[is_zero] = np.inf
+    spread = upper + lower
+    h = spread / standardization
+    num_ties = np.sum(lower == 0.0)
+    if num_ties:
+        replacements = np.ones((num_ties, num_ties)) - np.eye(num_ties)
+        replacements -= 2 * np.triu(replacements)
+        replacements = np.fliplr(replacements)
+        h[:num_ties, -num_ties:] = replacements
+    return np.median(h)


Test patch:
diff --git a/woodwork/tests/accessor/test_statistics.py b/woodwork/tests/accessor/test_statistics.py
index bef7c3f..9856df7 100644
--- a/woodwork/tests/accessor/test_statistics.py
+++ b/woodwork/tests/accessor/test_statistics.py
@@ -1,3 +1,4 @@
+import math
 import re
 import sys
 from datetime import datetime
@@ -7,6 +8,7 @@ from unittest.mock import patch
 import numpy as np
 import pandas as pd
 import pytest
+from scipy.stats import skew
 
 from woodwork.accessor_utils import _is_spark_dataframe, init_series
 from woodwork.config import CONFIG_DEFAULTS, config
@@ -41,11 +43,17 @@ from woodwork.statistics_utils import (
     _convert_ordinal_to_numeric,
     _get_describe_dict,
     _get_histogram_values,
+    _get_low_high_bound,
+    _get_medcouple_statistic,
     _get_mode,
     _get_numeric_value_counts_in_range,
     _get_recent_value_counts,
     _get_top_values_categorical,
 )
+from woodwork.statistics_utils._get_box_plot_info_for_column import (
+    _determine_best_outlier_method,
+    _determine_coefficients,
+)
 from woodwork.statistics_utils._parse_measures import _parse_measures
 from woodwork.tests.testing_utils import (
     _check_close,
@@ -1646,6 +1654,7 @@ def test_box_plot_outliers_with_quantiles(outliers_df):
     has_outliers_dict = outliers_series.ww.box_plot_dict(
         quantiles={0.0: -16.0, 0.25: 36.25, 0.5: 42.0, 0.75: 55.0, 1.0: 93.0},
     )
+    assert has_outliers_dict["method"] == "box_plot"
     assert has_outliers_dict["low_bound"] == 8.125
     assert has_outliers_dict["high_bound"] == 83.125
     assert has_outliers_dict["quantiles"] == {
@@ -1660,10 +1669,10 @@ def test_box_plot_outliers_with_quantiles(outliers_df):
     assert has_outliers_dict["low_indices"] == [3]
     assert has_outliers_dict["high_indices"] == [0]
 
-    no_outliers_dict = no_outliers_series.ww.box_plot_dict(
+    no_outliers_dict = no_outliers_series.ww.medcouple_dict(
         quantiles={0.0: 23.0, 0.25: 36.25, 0.5: 42.0, 0.75: 55.0, 1.0: 60.0},
     )
-
+    assert no_outliers_dict["method"] == "medcouple"
     assert no_outliers_dict["low_bound"] == 23.0
     assert no_outliers_dict["high_bound"] == 60.0
     assert no_outliers_dict["quantiles"] == {
@@ -1679,12 +1688,12 @@ def test_box_plot_outliers_with_quantiles(outliers_df):
     assert no_outliers_dict["high_indices"] == []
 
 
-def test_get_outliers_for_column_with_nans(outliers_df):
+def test_get_outliers_for_column_with_nans_box_plot(outliers_df):
     contains_nans_series = outliers_df["has_outliers_with_nans"]
     contains_nans_series.ww.init()
 
     box_plot_dict = contains_nans_series.ww.box_plot_dict()
-
+    assert box_plot_dict["method"] == "box_plot"
     assert box_plot_dict["low_bound"] == 4.5
     assert box_plot_dict["high_bound"] == 88.5
     assert box_plot_dict["quantiles"] == {
@@ -1776,6 +1785,7 @@ def test_box_plot_different_quantiles(outliers_df):
     box_plot_info = has_outliers_series.ww.box_plot_dict()
 
     assert set(box_plot_info.keys()) == {
+        "method",
         "low_bound",
         "high_bound",
         "quantiles",
@@ -1784,6 +1794,7 @@ def test_box_plot_different_quantiles(outliers_df):
         "low_indices",
         "high_indices",
     }
+    assert box_plot_info["method"] == "box_plot"
     assert box_plot_info["low_bound"] == 8.125
     assert box_plot_info["high_bound"] == 83.125
     assert len(box_plot_info["quantiles"]) == 5
@@ -1795,6 +1806,7 @@ def test_box_plot_different_quantiles(outliers_df):
     box_plot_info = has_outliers_series.ww.box_plot_dict(quantiles=partial_quantiles)
 
     assert set(box_plot_info.keys()) == {
+        "method",
         "low_bound",
         "high_bound",
         "quantiles",
@@ -1816,6 +1828,7 @@ def test_box_plot_different_quantiles(outliers_df):
     box_plot_info = no_outliers_series.ww.box_plot_dict(quantiles=partial_quantiles)
 
     assert set(box_plot_info.keys()) == {
+        "method",
         "low_bound",
         "high_bound",
         "quantiles",
@@ -1864,10 +1877,11 @@ def test_box_plot_optional_return_values(outliers_df):
         include_indices_and_values=True,
     )
 
-    assert {"low_bound", "high_bound", "quantiles"} == set(
+    assert {"low_bound", "high_bound", "quantiles", "method"} == set(
         has_outliers_box_plot_info_without_optional.keys(),
     )
     assert {
+        "method",
         "low_bound",
         "high_bound",
         "quantiles",
@@ -1887,10 +1901,11 @@ def test_box_plot_optional_return_values(outliers_df):
         include_indices_and_values=True,
     )
 
-    assert {"low_bound", "high_bound", "quantiles"} == set(
+    assert {"low_bound", "high_bound", "quantiles", "method"} == set(
         no_outliers_box_plot_info_without_optional.keys(),
     )
     assert {
+        "method",
         "low_bound",
         "high_bound",
         "quantiles",
@@ -1901,6 +1916,236 @@ def test_box_plot_optional_return_values(outliers_df):
     } == set(no_outliers_box_plot_info_with_optional.keys())
 
 
+def test_medcouple_outliers(skewed_outliers_df):
+    outliers_series_skewed_right = skewed_outliers_df["right_skewed_outliers"]
+    outliers_series_skewed_right.ww.init()
+
+    outliers_series_skewed_left = skewed_outliers_df["left_skewed_outliers"]
+    outliers_series_skewed_left.ww.init()
+
+    right_skewed_dict = outliers_series_skewed_right.ww.medcouple_dict()
+    left_skewed_dict = outliers_series_skewed_left.ww.medcouple_dict()
+
+    assert set(right_skewed_dict.keys()) == {
+        "method",
+        "low_bound",
+        "high_bound",
+        "quantiles",
+        "low_values",
+        "high_values",
+        "low_indices",
+        "high_indices",
+        "medcouple_stat",
+    }
+
+    expected_right_skewed_dict = {
+        "method": "medcouple",
+        "low_bound": 1.58676,
+        "high_bound": 20.32873,
+        "quantiles": {
+            0.0: 1.0,
+            0.25: 3.0,
+            0.5: 4.0,
+            0.75: 6.0,
+            1.0: 30.0,
+        },
+        "low_values": [1, 1],
+        "high_values": [30],
+        "low_indices": [0, 1],
+        "high_indices": [65],
+        "medcouple_stat": 0.333,
+    }
+
+    expected_left_skewed_dict = {
+        "method": "medcouple",
+        "low_bound": 23.58676,
+        "high_bound": 30.0,
+        "quantiles": {
+            0.0: 1.0,
+            0.25: 25.0,
+            0.5: 27.0,
+            0.75: 28.0,
+            1.0: 30.0,
+        },
+        "low_values": [23, 23, 22, 22, 21, 20, 18, 17, 15, 1],
+        "high_values": [],
+        "low_indices": [56, 57, 58, 59, 60, 61, 62, 63, 64, 65],
+        "high_indices": [],
+        "medcouple_stat": -0.333,
+    }
+
+    assert right_skewed_dict == expected_right_skewed_dict
+    assert left_skewed_dict == expected_left_skewed_dict
+
+
+def test_medcouple_outliers_with_quantiles(skewed_outliers_df):
+    outliers_series_skewed_right = skewed_outliers_df["right_skewed_outliers"]
+    outliers_series_skewed_right.ww.init()
+
+    outliers_series_skewed_left = skewed_outliers_df["left_skewed_outliers"]
+    outliers_series_skewed_left.ww.init()
+
+    override_quantiles_right = {
+        0.0: 1.0,
+        0.25: 8.0,
+        0.5: 11.0,
+        0.75: 25.0,
+        1.0: 30.0,
+    }
+
+    override_quantiles_left = {
+        0.0: 1.0,
+        0.25: 3.0,
+        0.5: 4.0,
+        0.75: 26.0,
+        1.0: 30.0,
+    }
+
+    expected_skewed_dict = {
+        "method": "medcouple",
+        "low_bound": 1.0,
+        "high_bound": 30.0,
+        "quantiles": None,
+        "low_values": [],
+        "high_values": [],
+        "low_indices": [],
+        "high_indices": [],
+        "medcouple_stat": 0.333,
+    }
+
+    right_skewed_dict = outliers_series_skewed_right.ww.medcouple_dict(
+        quantiles=override_quantiles_right,
+    )
+    left_skewed_dict = outliers_series_skewed_left.ww.medcouple_dict(
+        quantiles=override_quantiles_left,
+    )
+
+    expected_skewed_dict["quantiles"] = override_quantiles_right
+    assert right_skewed_dict == expected_skewed_dict
+    expected_skewed_dict["quantiles"] = override_quantiles_left
+    expected_skewed_dict["medcouple_stat"] = -0.333
+    assert left_skewed_dict == expected_skewed_dict
+
+
+def test_get_outliers_for_column_with_nans_medcouple(skewed_outliers_df):
+    contains_nans_series = skewed_outliers_df["has_outliers_with_nans"]
+    contains_nans_series.ww.init()
+
+    medcouple_dict = contains_nans_series.ww.medcouple_dict()
+
+    expected_skewed_dict = {
+        "method": "medcouple",
+        "low_bound": 1.94779,
+        "high_bound": 16.0754,
+        "quantiles": {
+            0.0: 1.0,
+            0.25: 3.0,
+            0.5: 4.0,
+            0.75: 5.25,
+            1.0: 30.0,
+        },
+        "low_values": [1.0, 1.0],
+        "high_values": [30.0],
+        "low_indices": [0, 1],
+        "high_indices": [65],
+        "medcouple_stat": 0.333,
+    }
+
+    assert medcouple_dict == expected_skewed_dict
+
+
+@pytest.mark.parametrize("mc", [-1.0, -0.5, -0.1, 0, 0.3333333, 1.0])
+def test_determine_coefficients(mc, skewed_outliers_df):
+    if _is_spark_dataframe(skewed_outliers_df):
+        pytest.xfail("spark hasn't implemented __iter__() for series")
+    right_skewed = skewed_outliers_df["right_skewed_outliers"]
+    left_skewed = skewed_outliers_df["left_skewed_outliers"]
+
+    right_coeff = np.abs(skew(right_skewed))
+    right_coeff = min(right_coeff, 3.5)
+
+    left_coeff = np.abs(skew(left_skewed))
+    left_coeff = min(left_coeff, 3.5)
+
+    if mc >= 0:
+        assert _determine_coefficients(right_skewed, mc) == (-right_coeff, right_coeff)
+        assert _determine_coefficients(left_skewed, mc) == (-left_coeff, left_coeff)
+    else:
+        assert _determine_coefficients(right_skewed, mc) == (right_coeff, -right_coeff)
+        assert _determine_coefficients(left_skewed, mc) == (left_coeff, -left_coeff)
+
+
+def test_get_low_high_bound_warnings():
+    error = "If the method selected is medcouple, then mc cannot be None."
+    with pytest.raises(ValueError, match=error):
+        _get_low_high_bound(None, "medcouple", 1, 1, None, None, None)
+
+    error = "Acceptable methods are 'box_plot' and 'medcouple'. The value passed was 'something_else'"
+    with pytest.raises(ValueError, match=error):
+        _get_low_high_bound(None, "something_else", 1, 1, None, None, None)
+
+
+def test_get_medcouple(outliers_df_pandas, skewed_outliers_df_pandas):
+    has_outliers_series = outliers_df_pandas["has_outliers"]
+    has_outliers_series = has_outliers_series.append(pd.Series([39]), ignore_index=True)
+    has_outliers_series.ww.init()
+    mc = _get_medcouple_statistic(has_outliers_series)
+    assert mc == 0.122
+
+    outliers_series_skewed_right = skewed_outliers_df_pandas["right_skewed_outliers"]
+    outliers_series_skewed_right.ww.init()
+    mc = _get_medcouple_statistic(outliers_series_skewed_right)
+    assert mc == 0.333
+
+    outliers_series_skewed = skewed_outliers_df_pandas[
+        ["right_skewed_outliers", "left_skewed_outliers"]
+    ]
+    outliers_series_skewed.ww.init()
+    mc = _get_medcouple_statistic(outliers_series_skewed)
+    assert isinstance(mc, np.ndarray)
+    np.testing.assert_almost_equal(mc, np.array([0.33333333, -0.33333333]))
+
+
+def test_determine_best_outlier_method_sampling_outcome(skewed_outliers_df_pandas):
+    # Column of 66,000, far above the 10,000 limit
+    contains_nans_series_skewed = (
+        skewed_outliers_df_pandas["right_skewed_outliers"]
+        .repeat(1000)
+        .reset_index(drop=True)
+    )
+    contains_nans_series_skewed.ww.init()
+
+    mc_result = _determine_best_outlier_method(contains_nans_series_skewed)
+
+    assert mc_result.method == "medcouple"
+    assert math.isclose(mc_result.mc, 0.33, rel_tol=0.01)
+
+
+def test_determine_best_outlier_method_equivalent_outcome(
+    outliers_df_pandas, skewed_outliers_df_pandas
+):
+    contains_nans_series_skewed = skewed_outliers_df_pandas["right_skewed_outliers"]
+    contains_nans_series_skewed.ww.init()
+
+    contains_nans_series = outliers_df_pandas["has_outliers"]
+    contains_nans_series.ww.init()
+
+    outliers_mc_skewed = contains_nans_series_skewed.ww.get_outliers(method="medcouple")
+    outliers_best_skewed = contains_nans_series_skewed.ww.get_outliers(method="best")
+
+    outliers_bp = contains_nans_series.ww.get_outliers(method="box_plot")
+    outliers_best = contains_nans_series.ww.get_outliers(method="best")
+
+    assert "medcouple_stat" not in outliers_bp.keys()
+    assert "medcouple_stat" in outliers_mc_skewed.keys()
+
+    assert outliers_bp == outliers_best
+    assert _get_medcouple_statistic(contains_nans_series) < 0.3
+
+    assert outliers_mc_skewed == outliers_best_skewed
+    assert _get_medcouple_statistic(contains_nans_series_skewed) >= 0.3
+
+
 @patch.object(
     sys.modules["woodwork.statistics_utils._infer_temporal_frequencies"],
     "infer_frequency",
@@ -2165,12 +2410,22 @@ def test_box_plot_ignore_zeros_null(dtype):
     zeros_df.ww.init(logical_type=dtype)
     no_zeros_df.ww.init(logical_type=dtype)
 
+    zeros_box_ignored_skewed = zeros_df.ww.medcouple_dict(ignore_zeros=True)
+    zeros_box_not_ignored_skewed = zeros_df.ww.medcouple_dict()
+
     zeros_box_ignored = zeros_df.ww.box_plot_dict(ignore_zeros=True)
     zeros_box_not_ignored = zeros_df.ww.box_plot_dict()
 
+    no_zeros_box_ignored_skewed = no_zeros_df.ww.medcouple_dict(ignore_zeros=True)
+    no_zeros_box_not_ignored_skewed = no_zeros_df.ww.medcouple_dict()
+
     no_zeros_box_ignored = no_zeros_df.ww.box_plot_dict(ignore_zeros=True)
     no_zeros_box_not_ignored = no_zeros_df.ww.box_plot_dict()
 
+    assert zeros_box_ignored_skewed == no_zeros_box_ignored_skewed
+    assert zeros_box_ignored_skewed == no_zeros_box_not_ignored_skewed
+    assert zeros_box_not_ignored_skewed != no_zeros_box_ignored_skewed
+
     assert zeros_box_ignored == no_zeros_box_ignored
     assert zeros_box_ignored == no_zeros_box_not_ignored
     assert zeros_box_not_ignored != no_zeros_box_ignored
diff --git a/woodwork/tests/conftest.py b/woodwork/tests/conftest.py
index 98565af..095327b 100644
--- a/woodwork/tests/conftest.py
+++ b/woodwork/tests/conftest.py
@@ -971,6 +971,59 @@ def outliers_df(request):
     return request.getfixturevalue(request.param)
 
 
+@pytest.fixture()
+def skewed_outliers_df_pandas():
+    outliers_df = pd.DataFrame(
+        {
+            "right_skewed_outliers": [1] * 2
+            + [2] * 6
+            + [3] * 20
+            + [4] * 12
+            + [5] * 8
+            + [6] * 5
+            + [7] * 3
+            + [8, 8, 9, 9, 10, 11, 13, 14, 16, 30],
+            "no_outliers": [60, 42, 37, 23, 49, 42, 36, 57, 60, 23.0] * 6
+            + [35, 54, 43, 47, 41, 39],
+            "non_numeric": ["a"] * 66,
+            "has_outliers_with_nans": [1] * 2
+            + [2] * 6
+            + [3] * 20
+            + [4] * 12
+            + [5] * 8
+            + [6] * 5
+            + [7] * 3
+            + [8, None, 9, 9, 10, 11, None, 14, 16, 30],
+            "nans": pd.Series([None] * 66, dtype="float64"),
+        },
+    )
+    outliers_df["left_skewed_outliers"] = 31 - outliers_df["right_skewed_outliers"]
+    return outliers_df
+
+
+@pytest.fixture()
+def skewed_outliers_df_dask(skewed_outliers_df_pandas):
+    dd = pytest.importorskip("dask.dataframe", reason="Dask not installed, skipping")
+    return dd.from_pandas(skewed_outliers_df_pandas, npartitions=2)
+
+
+@pytest.fixture()
+def skewed_outliers_df_spark(skewed_outliers_df_pandas):
+    ps = pytest.importorskip("pyspark.pandas", reason="Spark not installed, skipping")
+    return ps.from_pandas(skewed_outliers_df_pandas)
+
+
+@pytest.fixture(
+    params=[
+        "skewed_outliers_df_pandas",
+        "skewed_outliers_df_dask",
+        "skewed_outliers_df_spark",
+    ]
+)
+def skewed_outliers_df(request):
+    return request.getfixturevalue(request.param)
+
+
 class MockCallback:
     def __init__(self):
         self.progress_history = []