Repository: scipy/scipy
Pull request: Make median_absolute_deviation scale argument aligned w/iqr
#### Reference issue
Closes gh-11090 -> #11090 

#### What does this implement/fix?
In gh-11090 there was some confusion around the default value of scale in `median_absolute_deviation` function and there is also an inconsistency between how 
`median_absolute_deviation` and `iqr` accept `scale=` argument. 

In `iqr` the scale can be either a string or a float. If string it checks if lower cased is 'raw' or 'normal' with 'raw' as default which corresponds to 1.0 float value. 

This PR adds similar functionality to `median_absolute_deviation` and makes the bevavior aligned with `iqr`.

#### Additional information
Note the scaling used for `iqr` when `scale='normal'` is used is ***different*** than what is needed in `median_absolute_deviation` is a normal random variable is assumed as input (the reason for the choice of keyword) and a normalization is required to provide a consistent estimator of the standard deviation. 

Rather than duplicate the scale dictionary for `median_absolute_deviation` from `iqr`:
https://github.com/scipy/scipy/blob/a903652679ecbe787449fb99b8f3a6e105ad39b0/scipy/stats/stats.py#L2602
I chose to reuse the existing one created for `iqr` and I did a rescale for `median_absolute_deviation` if `scale='normal'` to match the loss function/use case. 

Problem statement:
Misconception of the median absolute deviation in stats?
As far as I know, the median absolute deviation (MAD) is simply defined to be 
```MAD = median(|x - median(x)|)```
However, `stats.median_absolute_deviation` is misleading in that it returns `1.4826 * MAD`, which is really an approximation for the standard deviation of a gaussian distributed data set, not the MAD itself. 
I would suggest setting the default scalar to 1 (currently 1.4826) to avoid further confusion..!

Hints from issue discussion:
hi @tcjansen, in statistics, for better or worse, many things are compared against the Gaussian distribution. Here is a reference:
https://en.wikipedia.org/wiki/Median_absolute_deviation#Relation_to_standard_deviation

I agree that this should perhaps be better documented to make it more clear to the user. 
You can override the default value by changing the `scale=` argument  from the default. 
Here is a link to where that exists in the code: https://github.com/scipy/scipy/blob/f3ef128eaa8455e65d80b00c21422938f4935b13/scipy/stats/stats.py#L2767

Are you interested in opening a pull request to clarify the documentation? 

@rlucas7 In my opinion, the documentation is not the place to correct for a misleading function name. With the default `scale=1.4826` value, this is clearly not the MAD! I shouldn't have to read the documentation to figure out if a function named `median_absolute_deviation` actually returns the median absolute deviation. Either the default value should be changed to `scale=1.0`, or the name of the function should be changed to `median_absolute_deviation_approximation_to_gaussian_std`, which seems silly to me.
@u55 Yes, these are my thoughts exactly.
@u55 and @tcjansen thanks for your comments. The original proposal for the function was for consistency with an existing function in the `statsmodels` package which is used to provide a robust standard error estimate. One of the `statsmodels` maintainers noted in the issue (linked beneath) that there were other packages using `statsmodels` as a dependency solely for the `mad` function, the abbreviation used in that package. 

The name is overloaded in statistics (as are many names) to describe the ***loss function*** as well as the ***robust measure of scale*** which explains the confusion. The implementation here is for using a scaling factor for the robust measure of scale for a Gaussian as the default and using the value 1 as you both propose would be for the loss function calculation, or perhaps a robust measure of scale for some other probability distribution where the constant 1 provides a consistent estimator. 

You can read the original proposal:
https://github.com/scipy/scipy/issues/9635

And you can read the pull request here:
https://github.com/scipy/scipy/pull/9637

And the corresponding function in the `statsmodels` package here: 
https://github.com/statsmodels/statsmodels/blob/fa16d00b8e66ee75a8045856ece7e7ac449631db/statsmodels/robust/scale.py#L18-L49
which describes the similar functionality.

A couple of concerns here would be: 
1. Making this change would have different defaults from `statsmodels` and might confuse other users who are using the method as a robust scale estimate. Or who are familiar with the `statsmodels` implementation. 
2. The function is already out in 1.3.X release so systems which took a dependency on that default would receive a breaking change if there was a switch. 
3. I think the usual approach is to issue a warning for 2 or more releases and then make the change. At least that is what we usually do for deprecations.

Given the consideration of 1 and 2, perhaps an alternative is to make the scale argument a ***required*** rather than an ***optional*** value and not have a default? 

This would remove the ambiguity and confusion on the part of the end user.
This would also give us a chance to change the name, I regret the full spelling as it is a lot of keystrokes. I'd rather have a name like `med_abs_dev()` instead.
iqr has options 'raw' as default and 'normal'

It really depends on the context and usage whether one or the other scaling is more convenient.
As dispersion measures or robust scale measures it is convenient to have var, iqr and mad on the same scale. That's the main usage in statsmodels robust, and the defaults are for that use case.

For scipy it's ambiguous, for scipy.stats normal distribution scaled mad is more appropriate if they are considered as dispersion measures.
If it were in numpy, then mad would just be raw MAD,

However, I think it would be better if iqr and mad follow the same pattern for the default.
@josef-pkt thanks for your comments, they are very helpful. I agree it is ambiguous w/scipy relative to other related packages.

>iqr has options 'raw' as default and 'normal'

Good call out, I hadn't considered the defaults of that method. I see what you mean:
https://github.com/scipy/scipy/blob/master/scipy/stats/stats.py#L2601-L2602

>However, I think it would be better if iqr and mad follow the same pattern for the default.

I agree that the interface here should be consistent with existing methods inside the package that are also robust measures of scale. So the change would be to make the `scale` argument in `median_absolute_deviation()` be consistent with `iqr()`, maintain the name and switch the default to 'raw' but also allow values to be specified by user, similar to what `iqr` does in scipy here:

https://github.com/scipy/scipy/blob/master/scipy/stats/stats.py#L2736-L2740

@tcjansen and @u55 would either of you be interested in opening a pull request to implement the change?
@rlucas7 

> 1. Making this change would have different defaults from statsmodels and might confuse other users who are using the method as a robust scale estimate. Or who are familiar with the statsmodels implementation.

I disagree with this justification for the current scipy API for the following reasons:
1. Scipy has a much larger user base than statsmodels, and users unfamiliar with statsmodels are _not_ going to expect the MAD to be multiplied by a magic number by default.
2. A poor API choice in statsmodels is not a good reason for scipy to replicate that bad design.
3. The current scipy implementation of `median_absolute_deviation` is _multiplied_ by its `scale` parameter, making it inconsistent with `statsmodels.robust.scale.mad`, which is _divided_ by its `c` parameter.
4. Scipy's default `scale=1.4826` value is a poor approximation to the exact value `1/c = 1/scipy.stats.norm.ppf(3/4.)` that `statsmodels.robust.scale.mad` uses, again making it inconsistent with statsmodels.

> So the change would be to make the scale argument in median_absolute_deviation() be consistent with iqr()

Apart from changing the default value of `scale` to `'raw'`, which I am in favor of, making `scale` in `median_absolute_deviation` consistent with `iqr` is even more of a breaking change since, currently, `iqr` is divided by its `scale` value and `median_absolute_deviation` is multiplied by its `scale` value. Making them consistent at this point, even with a warning cycle, seems just too confusing. If you want to make them consistent, I wonder if it wouldn't be better to just deprecate the entire function in favor of a differently-named function with a better API.

mad is in `scipy.stats`, and not in numpy or scipy.special, so it's perfectly valid to have mainly a statistical interpretation for it, or defaults for the main usage in statistics.

In statistic, I find unscaled iqr and mad pretty meaningless or useless, because the raw values don't have a scale that can be compared to anything else.
(e.g. use iqr as robust scale estimator in KDE)



@josef-pkt I'm not arguing that the MAD doesn't have a statistical interpretation. I am simply arguing that multiplying it by a magic number means that it is not the MAD anymore, so calling it the `median_absolute_deviation` is misleading! "Main usage" must always take a back seat to the principle of least surprise. Statistics users who want to use the MAD to estimate the Gaussian standard deviation are more likely to see a function named `scipy.stats.median_absolute_deviation` and multiply it by 1.4826 themselves (erroneously applying the scaling twice) not realizing that scipy is trying to be too cute. An optional scaling parameter that changes the meaning of the function should be a _convenience_ that one has to opt-in to, not a _default_ that one has to opt-out of.
in R https://www.rdocumentation.org/packages/stats/versions/3.6.1/topics/mad is scaled
NIST/Dataplot has mad and scaled mad
I cannot figure out from the SAS documentation what they are doing, the basic mad function is "raw"
**edit** new SAS docs: `MAD` computes both and raw by default
https://documentation.sas.com/?docsetId=imlug&docsetTarget=imlug_langref_sect259.htm&docsetVersion=14.2&locale=en
Stata doesn't seem to have it, user provided mad use scaled

"Although practicality beats purity."  Zen 9

I like the iqr way of specifying "normal"
In statsmodels and R you have to specify a number. It is easy to specify 1 for unscaled, but I wouldn't want to have to find the number for "normal" each time I use it.

It's still a `mad` Python function that allows us to compute scaled and unscaled mad. The default is just defined for the most common use cases in statistics.
i.e. I wouldn't change it in statsmodels to something that satisfies a "pure" definition by default.
But I'm indifferent to the default in scipy.stats. IMO both defaults can be justified in `stats`.

We have a lot of "magic" numbers in `statsmodels.robust` to calibrate it to the normal case. I only looked up a few of those as developer or user. 

(Aside: This reminds me of some function in scipy.special, where it also wasn't obvious whether the functions were scaled or not and if scaled by what. This was a long time ago when the scipy.special docs where half a sentence per function. Is there a `2 pi` in hermite polynomials? i.e. are they scaled for normal distribution or "raw". AFAIR)

Commit message: MAINT: make 'raw' default scale behavior and clarify docstring and test

TST: change default scaling behaviour in median_absolute_deviation

  * This assumes the change of scale in tests would be from
    `\sqrt(2)\times erfinv(0.5)` to 1. See gh-issue 11090
    for conversation on the change.

MAINT: make signature of median_absolute_deviation aligned w/iqr

DOC: make 'raw' docstring for median_absolute_deviation

DOC: fix doctests with proposed new scale default 1

TST: group in test class and add test for scale arg changes

MAINT: make 'raw' default scale behaviorand clarify docstring

DEP: deprecate median_absolute_deviation

  * This changeset deprecates median_absolute_deviation in favor
    of median_abs_devation. The latter will have a scale by 1 default
    whereas the former has a scale by 1.4876 default which was an
    imprecise approximation to the normal distribution scale factor.

DOC: add median_abs_deviation listing to refguide

MAINT: change how deprecation is called backwards compatible defaults

MAINT: clean up E501 and E128 errors and reorder key check-throw

MAINT: deprecate scale='raw' in iqr and median_abs_deviation

STY: fix pep8 issues from linter no whitespace after comma

DOC: add newlines and format docstring

DOC: add format errors from doc build

FIX: Work around numpy/gh-16202

FIX: Better stats page (single line)

MAINT: remove depr warning from new function and clean up

MAINT: fix scale and clean up pep8 and docstrings

clarify deprecation scale and language

TST: fix test name collision

Patch to review:
diff --git a/scipy/stats/__init__.py b/scipy/stats/__init__.py
index cd4849b00201..10ce7a356e7c 100644
--- a/scipy/stats/__init__.py
+++ b/scipy/stats/__init__.py
@@ -201,6 +201,7 @@
    mvsdist
    entropy
    median_absolute_deviation
+   median_abs_deviation
 
 Frequency statistics
 ====================
diff --git a/scipy/stats/stats.py b/scipy/stats/stats.py
index 62224225796b..8d36f487a3df 100644
--- a/scipy/stats/stats.py
+++ b/scipy/stats/stats.py
@@ -84,7 +84,7 @@
     zscore
     gstd
     iqr
-    median_absolute_deviation
+    median_abs_deviation
 
 Trimming Functions
 ------------------
@@ -193,6 +193,7 @@
            'scoreatpercentile', 'percentileofscore',
            'cumfreq', 'relfreq', 'obrientransform',
            'sem', 'zmap', 'zscore', 'iqr', 'gstd', 'median_absolute_deviation',
+           'median_abs_deviation',
            'sigmaclip', 'trimboth', 'trim1', 'trim_mean',
            'f_oneway', 'F_onewayConstantInputWarning',
            'F_onewayBadInputSizesWarning',
@@ -2650,7 +2651,7 @@ def gstd(a, axis=0, ddof=1):
                       'normal': special.erfinv(0.5) * 2.0 * math.sqrt(2.0)}
 
 
-def iqr(x, axis=None, rng=(25, 75), scale='raw', nan_policy='propagate',
+def iqr(x, axis=None, rng=(25, 75), scale=1.0, nan_policy='propagate',
         interpolation='linear', keepdims=False):
     r"""
     Compute the interquartile range of the data along the specified axis.
@@ -2686,7 +2687,8 @@ def iqr(x, axis=None, rng=(25, 75), scale='raw', nan_policy='propagate',
           'raw' : No scaling, just return the raw IQR.
           'normal' : Scale by :math:`2 \sqrt{2} erf^{-1}(\frac{1}{2}) \approx 1.349`.
 
-        The default is 'raw'. Array-like scale is also allowed, as long
+        The default is 1.0. The use of scale='raw' is deprecated.
+        Array-like scale is also allowed, as long
         as it broadcasts correctly to the output such that
         ``out / scale`` is a valid operation. The output dimensions
         depend on the input array, `x`, the `axis` argument, and the
@@ -2785,6 +2787,11 @@ def iqr(x, axis=None, rng=(25, 75), scale='raw', nan_policy='propagate',
         scale_key = scale.lower()
         if scale_key not in _scale_conversions:
             raise ValueError("{0} not a valid scale for `iqr`".format(scale))
+        if scale_key == 'raw':
+            warnings.warn(
+                "use of scale='raw' is deprecated, use scale=1.0 instead",
+                np.VisibleDeprecationWarning
+                )
         scale = _scale_conversions[scale_key]
 
     # Select the percentile function to use based on nans and policy
@@ -2812,9 +2819,9 @@ def iqr(x, axis=None, rng=(25, 75), scale='raw', nan_policy='propagate',
     return out
 
 
-def median_absolute_deviation(x, axis=0, center=np.median, scale=1.4826,
-                              nan_policy='propagate'):
-    """
+def median_abs_deviation(x, axis=0, center=np.median, scale=1.0,
+                         nan_policy='propagate'):
+    r"""
     Compute the median absolute deviation of the data along the given axis.
 
     The median absolute deviation (MAD, [1]_) computes the median over the
@@ -2823,7 +2830,7 @@ def median_absolute_deviation(x, axis=0, center=np.median, scale=1.4826,
 
     The MAD of an empty array is ``np.nan``.
 
-    .. versionadded:: 1.3.0
+    .. versionadded:: 1.5.0
 
     Parameters
     ----------
@@ -2834,19 +2841,25 @@ def median_absolute_deviation(x, axis=0, center=np.median, scale=1.4826,
         the MAD over the entire array.
     center : callable, optional
         A function that will return the central value. The default is to use
-        np.median. Any user defined function used will need to have the function
-        signature ``func(arr, axis)``.
-    scale : int, optional
-        The scaling factor applied to the MAD. The default scale (1.4826)
-        ensures consistency with the standard deviation for normally distributed
-        data.
+        np.median. Any user defined function used will need to have the
+        function signature ``func(arr, axis)``.
+    scale : scalar or str, optional
+        The numerical value of scale will be divided out of the final
+        result. The following string value is also recognized (default is 1.0):
+
+        * 'normal' : Scale by :math:`\frac{1}{\sqrt{2}erf^{-1}(1/2)}`.
+
+        Array-like scale is also allowed, as long
+        as it broadcasts correctly to the output such that
+        ``out / scale`` is a valid operation. The output dimensions
+        depend on the input array, `x`, and the `axis` argument.
     nan_policy : {'propagate', 'raise', 'omit'}, optional
         Defines how to handle when input contains nan.
         The following options are available (default is 'propagate'):
 
-          * 'propagate': returns nan
-          * 'raise': throws an error
-          * 'omit': performs the calculations ignoring nan values
+        * 'propagate': returns nan
+        * 'raise': throws an error
+        * 'omit': performs the calculations ignoring nan values
 
     Returns
     -------
@@ -2870,12 +2883,14 @@ def median_absolute_deviation(x, axis=0, center=np.median, scale=1.4826,
 
     References
     ----------
-    .. [1] "Median absolute deviation" https://en.wikipedia.org/wiki/Median_absolute_deviation
-    .. [2] "Robust measures of scale" https://en.wikipedia.org/wiki/Robust_measures_of_scale
+    .. [1] "Median absolute deviation",
+           https://en.wikipedia.org/wiki/Median_absolute_deviation
+    .. [2] "Robust measures of scale",
+           https://en.wikipedia.org/wiki/Robust_measures_of_scale
 
     Examples
     --------
-    When comparing the behavior of `median_absolute_deviation` with ``np.std``,
+    When comparing the behavior of `median_abs_deviation` with ``np.std``,
     the latter is affected when we change a single value of an array to have an
     outlier value while the MAD hardly changes:
 
@@ -2883,13 +2898,13 @@ def median_absolute_deviation(x, axis=0, center=np.median, scale=1.4826,
     >>> x = stats.norm.rvs(size=100, scale=1, random_state=123456)
     >>> x.std()
     0.9973906394005013
-    >>> stats.median_absolute_deviation(x)
-    1.2280762773108278
+    >>> stats.median_abs_deviation(x)
+    0.82832610097857
     >>> x[0] = 345.6
     >>> x.std()
     34.42304872314415
-    >>> stats.median_absolute_deviation(x)
-    1.2340335571164334
+    >>> stats.median_abs_deviation(x)
+    0.8323442311590675
 
     Axis handling example:
 
@@ -2897,14 +2912,30 @@ def median_absolute_deviation(x, axis=0, center=np.median, scale=1.4826,
     >>> x
     array([[10,  7,  4],
            [ 3,  2,  1]])
-    >>> stats.median_absolute_deviation(x)
-    array([5.1891, 3.7065, 2.2239])
-    >>> stats.median_absolute_deviation(x, axis=None)
-    2.9652
+    >>> stats.median_abs_deviation(x)
+    array([3.5, 2.5, 1.5])
+    >>> stats.median_abs_deviation(x, axis=None)
+    2.0
+
+    Scale normal example:
+
+    >>> x = stats.norm.rvs(size=1000000, scale=2, random_state=123456)
+    >>> stats.median_abs_deviation(x)
+    1.3487398527041636
+    >>> stats.median_abs_deviation(x, scale='normal')
+    1.9996446978061115
 
     """
     x = asarray(x)
 
+    # An error may be raised here, so fail-fast, before doing lengthy
+    # computations, even though `scale` is not used until later
+    if isinstance(scale, str):
+        if scale.lower() == 'normal':
+            scale = 0.6744897501960817  # special.ndtri(0.75)
+        else:
+            raise ValueError(f"{scale} is not a valid scale value.")
+
     # Consistent with `np.var` and `np.std`.
     if not x.size:
         nan_shape = [item for i, item in enumerate(x.shape) if i != axis]
@@ -2932,13 +2963,134 @@ def median_absolute_deviation(x, axis=0, center=np.median, scale=1.4826,
         med = np.apply_over_axes(center, arr, axis)
         mad = np.median(np.abs(arr - med), axis=axis)
 
-    return scale * mad
+    return mad / scale
+
+
+# Keep the top newline so that the message does not show up on the stats page
+_median_absolute_deviation_deprec_msg = """
+To preserve the existing default behavior, use
+`scipy.stats.median_abs_deviation(..., scale=1/1.4826)`.
+The value 1.4826 is not numerically precise for scaling
+with a normal distribution. For a numerically precise value, use
+`scipy.stats.median_abs_deviation(..., scale='normal')`.
+"""
+
+
+# Due to numpy/gh-16349 we need to unindent the entire docstring
+@np.deprecate(old_name='median_absolute_deviation',
+              new_name='median_abs_deviation',
+              message=_median_absolute_deviation_deprec_msg)
+def median_absolute_deviation(x, axis=0, center=np.median, scale=1.4826,
+                              nan_policy='propagate'):
+    r"""
+Compute the median absolute deviation of the data along the given axis.
 
+The median absolute deviation (MAD, [1]_) computes the median over the
+absolute deviations from the median. It is a measure of dispersion
+similar to the standard deviation but more robust to outliers [2]_.
+
+The MAD of an empty array is ``np.nan``.
+
+.. versionadded:: 1.3.0
+
+Parameters
+----------
+x : array_like
+    Input array or object that can be converted to an array.
+axis : int or None, optional
+    Axis along which the range is computed. Default is 0. If None, compute
+    the MAD over the entire array.
+center : callable, optional
+    A function that will return the central value. The default is to use
+    np.median. Any user defined function used will need to have the function
+    signature ``func(arr, axis)``.
+scale : int, optional
+    The scaling factor applied to the MAD. The default scale (1.4826)
+    ensures consistency with the standard deviation for normally distributed
+    data.
+nan_policy : {'propagate', 'raise', 'omit'}, optional
+    Defines how to handle when input contains nan.
+    The following options are available (default is 'propagate'):
+
+    * 'propagate': returns nan
+    * 'raise': throws an error
+    * 'omit': performs the calculations ignoring nan values
+
+Returns
+-------
+mad : scalar or ndarray
+    If ``axis=None``, a scalar is returned. If the input contains
+    integers or floats of smaller precision than ``np.float64``, then the
+    output data-type is ``np.float64``. Otherwise, the output data-type is
+    the same as that of the input.
+
+See Also
+--------
+numpy.std, numpy.var, numpy.median, scipy.stats.iqr, scipy.stats.tmean,
+scipy.stats.tstd, scipy.stats.tvar
+
+Notes
+-----
+The `center` argument only affects the calculation of the central value
+around which the MAD is calculated. That is, passing in ``center=np.mean``
+will calculate the MAD around the mean - it will not calculate the *mean*
+absolute deviation.
+
+References
+----------
+.. [1] "Median absolute deviation",
+       https://en.wikipedia.org/wiki/Median_absolute_deviation
+.. [2] "Robust measures of scale",
+       https://en.wikipedia.org/wiki/Robust_measures_of_scale
+
+Examples
+--------
+When comparing the behavior of `median_absolute_deviation` with ``np.std``,
+the latter is affected when we change a single value of an array to have an
+outlier value while the MAD hardly changes:
+
+>>> from scipy import stats
+>>> x = stats.norm.rvs(size=100, scale=1, random_state=123456)
+>>> x.std()
+0.9973906394005013
+>>> stats.median_absolute_deviation(x)
+1.2280762773108278
+>>> x[0] = 345.6
+>>> x.std()
+34.42304872314415
+>>> stats.median_absolute_deviation(x)
+1.2340335571164334
+
+Axis handling example:
+
+>>> x = np.array([[10, 7, 4], [3, 2, 1]])
+>>> x
+array([[10,  7,  4],
+       [ 3,  2,  1]])
+>>> stats.median_absolute_deviation(x)
+array([5.1891, 3.7065, 2.2239])
+>>> stats.median_absolute_deviation(x, axis=None)
+2.9652
+"""
+    if isinstance(scale, str):
+        if scale.lower() == 'raw':
+            warnings.warn(
+                "use of scale='raw' is deprecated, use scale=1.0 instead",
+                np.VisibleDeprecationWarning
+                )
+            scale = 1.0
+
+    if not isinstance(scale, str):
+        scale = 1 / scale
+
+    return median_abs_deviation(x, axis=axis, center=center, scale=scale,
+                                nan_policy=nan_policy)
 
 #####################################
 #         TRIMMING FUNCTIONS        #
 #####################################
 
+
 SigmaclipResult = namedtuple('SigmaclipResult', ('clipped', 'lower', 'upper'))
 
 
diff --git a/scipy/stats/tests/test_stats.py b/scipy/stats/tests/test_stats.py
index 270e03102431..6329d6c8d2ea 100644
--- a/scipy/stats/tests/test_stats.py
+++ b/scipy/stats/tests/test_stats.py
@@ -1876,65 +1876,99 @@ def test_zscore_nan_raise(self):
 
         assert_raises(ValueError, stats.zscore, x, nan_policy='raise')
 
-    def test_mad(self):
-        dat = np.array([2.20, 2.20, 2.4, 2.4, 2.5, 2.7, 2.8, 2.9, 3.03,
-                3.03, 3.10, 3.37, 3.4, 3.4, 3.4, 3.5, 3.6, 3.7, 3.7,
-                3.7, 3.7,3.77, 5.28, 28.95])
-        assert_almost_equal(stats.median_absolute_deviation(dat, axis=None), 0.526323)
 
-        dat = dat.reshape(6, 4)
-        mad = stats.median_absolute_deviation(dat, axis=0)
-        mad_expected = np.asarray([0.644931, 0.7413, 0.66717, 0.59304])
+class TestMedianAbsDeviation(object):
+    def setup_class(self):
+        self.dat_nan = np.array([2.20, 2.20, 2.4, 2.4, 2.5, 2.7, 2.8, 2.9,
+                3.03, 3.03, 3.10, 3.37, 3.4, 3.4, 3.4, 3.5, 3.6, 3.7, 3.7,
+                3.7, 3.7, 3.77, 5.28, np.nan])
+        self.dat = np.array([2.20, 2.20, 2.4, 2.4, 2.5, 2.7, 2.8, 2.9, 3.03,
+                3.03, 3.10, 3.37, 3.4, 3.4, 3.4, 3.5, 3.6, 3.7, 3.7,
+                3.7, 3.7, 3.77, 5.28, 28.95])
+
+    def test_median_abs_deviation(self):
+        assert_almost_equal(stats.median_abs_deviation(self.dat, axis=None),
+                            0.355)
+        dat = self.dat.reshape(6, 4)
+        mad = stats.median_abs_deviation(dat, axis=0)
+        mad_expected = np.asarray([0.435, 0.5, 0.45, 0.4])
         assert_array_almost_equal(mad, mad_expected)
 
+    def test_mad_nan_omit(self):
+        mad = stats.median_abs_deviation(self.dat_nan, nan_policy='omit')
+        assert_almost_equal(mad, 0.34)
+
+
+class TestMedianAbsoluteDeviation(object):
+    def setup_class(self):
+        self.dat_nan = np.array([2.20, 2.20, 2.4, 2.4, 2.5, 2.7, 2.8, 2.9, 3.03,
+                3.03, 3.10, 3.37, 3.4, 3.4, 3.4, 3.5, 3.6, 3.7, 3.7,
+                3.7, 3.7, 3.77, 5.28, np.nan])
+        self.dat = np.array([2.20, 2.20, 2.4, 2.4, 2.5, 2.7, 2.8, 2.9, 3.03,
+                3.03, 3.10, 3.37, 3.4, 3.4, 3.4, 3.5, 3.6, 3.7, 3.7,
+                3.7, 3.7, 3.77, 5.28, 28.95])
+
     def test_mad_empty(self):
         dat = []
-        mad = stats.median_absolute_deviation(dat)
+        with suppress_warnings() as sup:
+            sup.filter(DeprecationWarning)
+            mad = stats.median_absolute_deviation(dat)
         assert_equal(mad, np.nan)
 
     def test_mad_nan_shape1(self):
         z = np.ones((3, 0))
-        mad_axis0 = stats.median_absolute_deviation(z, axis=0)
+        with suppress_warnings() as sup:
+            sup.filter(DeprecationWarning)
+            mad_axis0 = stats.median_absolute_deviation(z, axis=0)
+            mad_axis1 = stats.median_absolute_deviation(z, axis=1)
         assert_equal(mad_axis0, np.nan)
-        mad_axis1 = stats.median_absolute_deviation(z, axis=1)
         assert_equal(mad_axis1, np.array([np.nan, np.nan, np.nan]))
         assert_equal(mad_axis1.shape, (3,))
 
     def test_mad_nan_shape2(self):
         z = np.ones((3, 0, 2))
-        mad_axis0 = stats.median_absolute_deviation(z, axis=0)
+        with suppress_warnings() as sup:
+            sup.filter(DeprecationWarning)
+            mad_axis0 = stats.median_absolute_deviation(z, axis=0)
+            mad_axis1 = stats.median_absolute_deviation(z, axis=1)
+            mad_axis2 = stats.median_absolute_deviation(z, axis=2)
         assert_equal(mad_axis0, np.nan)
-        mad_axis1 = stats.median_absolute_deviation(z, axis=1)
         assert_equal(mad_axis1, np.array([[np.nan, np.nan],
                                           [np.nan, np.nan],
                                           [np.nan, np.nan]]))
         assert_equal(mad_axis1.shape, (3, 2))
-        mad_axis2 = stats.median_absolute_deviation(z, axis=2)
         assert_equal(mad_axis2, np.nan)
 
     def test_mad_nan_propagate(self):
-        dat = np.array([2.20, 2.20, 2.4, 2.4, 2.5, 2.7, 2.8, 2.9, 3.03,
-                3.03, 3.10, 3.37, 3.4, 3.4, 3.4, 3.5, 3.6, 3.7, 3.7,
-                3.7, 3.7,3.77, 5.28, np.nan])
-
-        mad = stats.median_absolute_deviation(dat, nan_policy='propagate')
+        with suppress_warnings() as sup:
+            sup.filter(DeprecationWarning)
+            mad = stats.median_absolute_deviation(self.dat_nan,
+                                                  nan_policy='propagate')
         assert_equal(mad, np.nan)
 
     def test_mad_nan_raise(self):
-        dat = np.array([2.20, 2.20, 2.4, 2.4, 2.5, 2.7, 2.8, 2.9, 3.03,
-                3.03, 3.10, 3.37, 3.4, 3.4, 3.4, 3.5, 3.6, 3.7, 3.7,
-                3.7, 3.7,3.77, 5.28, np.nan])
-
         with assert_raises(ValueError):
-            stats.median_absolute_deviation(dat, nan_policy='raise')
+            with suppress_warnings() as sup:
+                sup.filter(DeprecationWarning)
+                stats.median_absolute_deviation(self.dat_nan,
+                                                nan_policy='raise')
 
-    def test_mad_nan_omit(self):
-        dat = np.array([2.20, 2.20, 2.4, 2.4, 2.5, 2.7, 2.8, 2.9, 3.03,
-                3.03, 3.10, 3.37, 3.4, 3.4, 3.4, 3.5, 3.6, 3.7, 3.7,
-                3.7, 3.7,3.77, 5.28, np.nan])
+    def test_mad_scale_default(self):
+        with suppress_warnings() as sup:
+            sup.filter(DeprecationWarning)
+            mad = stats.median_absolute_deviation(self.dat, scale=1.0)
+            mad_float = stats.median_absolute_deviation(self.dat, scale=1.0)
+        assert_almost_equal(mad, 0.355)
+        assert_almost_equal(mad, mad_float)
 
-        mad = stats.median_absolute_deviation(dat, nan_policy='omit')
-        assert_almost_equal(mad, 0.504084)
+    def test_mad_scale_normal(self):
+        with suppress_warnings() as sup:
+            sup.filter(DeprecationWarning)
+            mad = stats.median_absolute_deviation(self.dat, scale="normal")
+            scale = 1.4826022185056018
+            mad_float = stats.median_absolute_deviation(self.dat, scale=scale)
+        assert_almost_equal(mad, 0.526323787)
+        assert_almost_equal(mad, mad_float)
 
 
 def _check_warnings(warn_list, expected_type, expected_len):
@@ -1962,7 +1996,7 @@ def test_api(self):
         stats.iqr(d, 1)
         stats.iqr(d, (0, 1))
         stats.iqr(d, None, (10, 90))
-        stats.iqr(d, None, (30, 20), 'raw')
+        stats.iqr(d, None, (30, 20), 1.0)
         stats.iqr(d, None, (25, 75), 1.5, 'propagate')
         stats.iqr(d, None, (50, 50), 'normal', 'raise', 'linear')
         stats.iqr(d, None, (25, 75), -0.4, 'omit', 'lower', True)
@@ -2138,7 +2172,7 @@ def test_scale(self):
         x = np.arange(15.0).reshape((3, 5))
 
         # No NaNs
-        assert_equal(stats.iqr(x, scale='raw'), 7)
+        assert_equal(stats.iqr(x, scale=1.0), 7)
         assert_almost_equal(stats.iqr(x, scale='normal'), 7 / 1.3489795)
         assert_equal(stats.iqr(x, scale=2.0), 3.5)
 
@@ -2146,11 +2180,11 @@ def test_scale(self):
         x[1, 2] = np.nan
         with warnings.catch_warnings(record=True):
             warnings.simplefilter("always")
-            assert_equal(stats.iqr(x, scale='raw', nan_policy='propagate'), np.nan)
+            assert_equal(stats.iqr(x, scale=1.0, nan_policy='propagate'), np.nan)
             assert_equal(stats.iqr(x, scale='normal', nan_policy='propagate'), np.nan)
             assert_equal(stats.iqr(x, scale=2.0, nan_policy='propagate'), np.nan)
             # axis=1 chosen to show behavior with both nans and without
-            assert_equal(stats.iqr(x, axis=1, scale='raw',
+            assert_equal(stats.iqr(x, axis=1, scale=1.0,
                                    nan_policy='propagate'), [2, np.nan, 2])
             assert_almost_equal(stats.iqr(x, axis=1, scale='normal',
                                           nan_policy='propagate'),
@@ -2161,7 +2195,7 @@ def test_scale(self):
             # np.percentile with nans, so we don't check the number of
             # warnings here. See https://github.com/numpy/numpy/pull/12679.
 
-        assert_equal(stats.iqr(x, scale='raw', nan_policy='omit'), 7.5)
+        assert_equal(stats.iqr(x, scale=1.0, nan_policy='omit'), 7.5)
         assert_almost_equal(stats.iqr(x, scale='normal', nan_policy='omit'),
                             7.5 / 1.3489795)
         assert_equal(stats.iqr(x, scale=2.0, nan_policy='omit'), 3.75)


Code review comments:
diff_hunk: @@ -2686,7 +2687,8 @@ def iqr(x, axis=None, rng=(25, 75), scale='raw', nan_policy='propagate',
           'raw' : No scaling, just return the raw IQR.
           'normal' : Scale by :math:`2 \sqrt{2} erf^{-1}(\frac{1}{2}) \approx 1.349`.
path: scipy/stats/stats.py
text: @user1:
```suggestion
          * 'raw' : No scaling, just return the raw IQR.
            **Deprecated!**  Use `scale=1` instead.
          * 'normal' : Scale by
            :math:`2 \sqrt{2} erf^{-1}(\frac{1}{2}) \approx 1.349`.
```

diff_hunk: @@ -2834,19 +2841,25 @@ def median_absolute_deviation(x, axis=0, center=np.median, scale=1.4826,
         the MAD over the entire array.
     center : callable, optional
         A function that will return the central value. The default is to use
-        np.median. Any user defined function used will need to have the function
-        signature ``func(arr, axis)``.
-    scale : int, optional
-        The scaling factor applied to the MAD. The default scale (1.4826)
-        ensures consistency with the standard deviation for normally distributed
-        data.
+        np.median. Any user defined function used will need to have the
+        function signature ``func(arr, axis)``.
+    scale : scalar or str, optional
+        The numerical value of scale will be divided out of the final
+        result. The following string value is also recognized (default is 1.0):
+
+        * 'normal' : Scale by :math:`\frac{1}{\sqrt{2}erf^{-1}(1/2)}`.
+
+        Array-like scale is also allowed, as long
+        as it broadcasts correctly to the output such that
+        ``out / scale`` is a valid operation. The output dimensions
+        depend on the input array, `x`, and the `axis` argument.
path: scipy/stats/stats.py
text: @user1:
```suggestion
        The numerical value of scale will be divided out of the final
        result. The default is 1.0. The string "normal" is also accepted,
        and results in `scale` being the inverse of the standard normal
        quantile function at 0.75, which is approximately 0.67449.
        Array-like scale is also allowed, as long as it broadcasts correctly
        to the output such that ``out / scale`` is a valid operation. The
        output dimensions depend on the input array, `x`, and the `axis`
        argument.
```