Repository: scikit-learn/scikit-learn
Pull request: ENH Add option to scale output to unit var in RobustScaler
<!--
Thanks for contributing a pull request! Please ensure you have taken a look at
the contribution guidelines: https://github.com/scikit-learn/scikit-learn/blob/master/CONTRIBUTING.md#pull-request-checklist
-->

#### Reference Issues/PRs
Follows from PR #10140
closes #10139 

#### What does this implement/fix? Explain your changes.
Adds gauss_adjust parameter to `RobustScaler` so transformed data is standard Gaussian when the input data was a Gaussian distribution.

Adds test to check that this option makes StandardScaler and RobustScaler equivalent on normal data as suggested by @amueller. Note these are only equivalent with a large dataset. Not sure if the parameters I used are optimal, open to suggestions.

#### Any other comments?


<!--
Please be aware that we are a loose team of volunteers so patience is
necessary; assistance handling other issues is very welcome. We value
all user contributions, no matter how minor they are. If we are slow to
review, either the pull request needs some benchmarking, tinkering,
convincing, etc. or more likely the reviewers are simply busy. In either
case, we ask for your understanding during the review process.
For more information, see our FAQ on this topic:
http://scikit-learn.org/dev/faq.html#why-is-my-pull-request-not-getting-any-attention.

Thanks for contributing!
-->


Problem statement:
Feature Request: Gaussian adjust option in RobustScaler
It would be nice if sklearn.preprocessing.RobustScaler transformed data to a standard Gaussian when the input data was a Gaussian distribution.  Currently, with the default arguments of quantile_range=(25.0, 75.0), if Gaussian data is passed to RobustScaler, then the transformed data has:
mean ~= 0 and std ~= 1/(norm.ppf(.75) - norm.ppf(.25)) = 0.741
It would be nice if this were 0 and 1 to be more consistent with the resulting scales in StandardScaler.

I propose adding an extra option to the constructor gauss_adjust, which if True adjusts the scales accordingly:
```
self.scale_ = (q[1] - q[0])
```
on line 1057 of data.py needs to be
```
self.scale_ = (q[1] - q[0]) / self.adjust
```
where in the constructor
```
self.adjust = norm.ppf(q_max / 100.0) - norm.ppf(q_min / 100.0) if gauss_adjust else 1.0
```

Hints from issue discussion:
Can't you just specify a different quantile_range?

On 15 November 2017 at 04:08, rdturnermtl <notifications@github.com> wrote:

> It would be nice if sklearn.preprocessing.RobustScaler transformed data
> to a standard Gaussian when the input data was a Gaussian distribution.
> Currently, with the default arguments of quantile_range=(25.0, 75.0), if
> Gaussian data is passed to RobustScaler, then the transformed data has:
> mean ~= 0 and std ~= 1/(norm.ppf(.75) - norm.ppf(.25)) = 0.741
> It would be nice if this were 0 and 1 to be more consistent with the
> resulting scales in StandardScaler.
>
> I propose adding an extra option to the constructor gauss_adjust, which if
> True adjusts the scales accordingly:
>
> self.scale_ = (q[1] - q[0])
>
> on line 1057 of data.py needs to be
>
> self.scale_ = (q[1] - q[0]) / self.adjust
>
> where in the constructor
>
> self.adjust = norm.ppf(q_max / 100.0) - norm.ppf(q_min / 100.0) if gauss_adjust else 1.0
>
> —
> You are receiving this because you are subscribed to this thread.
> Reply to this email directly, view it on GitHub
> <https://github.com/scikit-learn/scikit-learn/issues/10139>, or mute the
> thread
> <https://github.com/notifications/unsubscribe-auth/AAEz65HZQODk6ibHsY3ikJndARhNirl7ks5s2cgXgaJpZM4QdrEm>
> .
>

mean=0 and std=0.741 is still a Gaussian :)
I am not sure what is the benefit in practice to have all features with an std of 0.741 or 1.0
If really this is an issue, one can always pipeline a robust scaler and standard scaler
> I am not sure what is the benefit in practice to have all features with an std of 0.741 or 1.0

Priors/regularization on regression weights depend on the scale of the input. If the scale of robust scale is smaller than standard scaler, then robust scaler implicitly results in stronger priors on the regression weights. This is resolved if the robust scaler gives around the same scales as standard scaler, but with outlier robustness.

> one can always pipeline a robust scaler and standard scaler

Pipelining with a standard scaler will undo the scaling of robust scaler and we lose the robustness to outliers.

> Can't you just specify a different quantile_range?

The desired standardization can be acheived with a smaller range:
    quantile_range=(30.9, 69.1)
However, different ranges have different statistical efficiencies:
    (30.9, 69.1) has 26% efficiency, (25.0, 75.0) has 37%.
However, a bigger range is what is needed if anything:
    (7.0, 93.0) has 66% efficiency.
Lower efficiency means there will be a bigger shift in scale between the training set and test set.
See:
https://en.wikipedia.org/wiki/Interdecile_range
https://en.wikipedia.org/wiki/Efficiency_(statistics)

This is why we would like the flexibility to specify a large inter-quartile range but also have around the same approximate scale to standard scaler.

Here is a nice little demo:
![eff](https://user-images.githubusercontent.com/28273671/32856144-f57935aa-ca11-11e7-8165-85d0d6078124.png)

```
import numpy as np
import scipy.stats as ss

np.random.seed(1234)

N = 1000
runs = 20000

w = np.linspace(1, 49, 500)
pp = np.vstack((50.0 - w, 50.0 + w))

normalizer = ss.norm.ppf(pp[1, :] / 100) - ss.norm.ppf(pp[0, :] / 100)

x = np.random.randn(N, runs)
v = np.percentile(x, pp[1, :], axis=0) - np.percentile(x, pp[0, :], axis=0)
v = v / normalizer[:, None]  # Now estimator of std_dev

base_var = np.var(np.std(x, axis=0))

mean_est = np.mean(v, axis=1)
var_est = np.var(v, axis=1)
eff = base_var / var_est

import matplotlib.pyplot as plt
plt.plot(2 * w, 100 * eff)
plt.xlabel('range')
plt.ylabel('efficiency')
plt.grid()
```
~but what your "gaussian adjust" is doing is equivalent to changing the quantiles, right?~ (no it's not)
One could make the point that you need to search for the correct regularization anyway, but having different default behavior because of robust scaling is indeed a bit odd. So I can see the point.
This is only relevant if the data is really Gaussian, though, right?
Yeah, this example uses a Gaussian as the default reference case, since that is the most obvious choice if we need to pick a default reference.  It seems a lot of the efficiency properties derived in robust statistics uses the same principle.

The distribution could be made adjustable too if we wanted to be fancy.  It looks like the implementation in #10140 already envisions that.
Another thought: Does the QuantileTransformer help your case?

> Does the QuantileTransformer help your case?

Well QuantileTransformer is a nonlinear preprocessor while RobustScaler is still just a linear transform so I think they are really different things.
They are, but you get more control of the output distribution while
minimising the impact of outliers in the input.

On 17 November 2017 at 03:45, rdturnermtl <notifications@github.com> wrote:

> Does the QuantileTransformer help your case?
>
> Well QuantileTransformer is a nonlinear preprocessor while RobustScaler is
> still just a linear transform so I think they are really different things.
>
> —
> You are receiving this because you commented.
> Reply to this email directly, view it on GitHub
> <https://github.com/scikit-learn/scikit-learn/issues/10139#issuecomment-344983888>,
> or mute the thread
> <https://github.com/notifications/unsubscribe-auth/AAEz63z82iG8hjl9V3fqOWAPup8WQuaVks5s3GbFgaJpZM4QdrEm>
> .
>

Commit message: add whats new

Patch to review:
diff --git a/doc/whats_new/v0.24.rst b/doc/whats_new/v0.24.rst
index 8661f8138f003..e949cff3d3316 100644
--- a/doc/whats_new/v0.24.rst
+++ b/doc/whats_new/v0.24.rst
@@ -65,6 +65,14 @@ Changelog
   samples between the train and test set on each fold.
   :pr:`13204` by :user:`Kyle Kosic <kykosic>`.
 
+:mod:`sklearn.preprocessing`
+..............................
+
+- |Enhancement| Add ``unit_variance`` to :class:`preprocessing.RobustScaler`,
+  which scales output data such that normally distributed features have a
+  variance of 1. :pr:`17193` by :user:`Lucy Liu <lucyleeow>` and
+  :user:`Mabel Villalba <mabelvj>`.
+
 
 Code and Documentation Contributors
 -----------------------------------
diff --git a/sklearn/preprocessing/_data.py b/sklearn/preprocessing/_data.py
index cc8776951f114..c68dff15e63d1 100644
--- a/sklearn/preprocessing/_data.py
+++ b/sklearn/preprocessing/_data.py
@@ -1130,6 +1130,10 @@ class RobustScaler(TransformerMixin, BaseEstimator):
 
         .. versionadded:: 0.18
 
+    unit_variance : boolean, False by default
+        If True, scale data so that normally distributed features have a
+        variance of 1.
+
     copy : boolean, optional, default is True
         If False, try to avoid a copy and do inplace scaling instead.
         This is not guaranteed to always work inplace; e.g. if the data is
@@ -1180,10 +1184,12 @@ class RobustScaler(TransformerMixin, BaseEstimator):
     """
     @_deprecate_positional_args
     def __init__(self, *, with_centering=True, with_scaling=True,
-                 quantile_range=(25.0, 75.0), copy=True):
+                 quantile_range=(25.0, 75.0), unit_variance=False,
+                 copy=True):
         self.with_centering = with_centering
         self.with_scaling = with_scaling
         self.quantile_range = quantile_range
+        self.unit_variance = unit_variance
         self.copy = copy
 
     def fit(self, X, y=None):
@@ -1233,6 +1239,12 @@ def fit(self, X, y=None):
 
             self.scale_ = quantiles[1] - quantiles[0]
             self.scale_ = _handle_zeros_in_scale(self.scale_, copy=False)
+            if self.unit_variance:
+                # Create scipy.stats.norm object
+                output_distribution = getattr(stats, 'norm')
+                self.adjust_ = output_distribution.ppf(q_max / 100.0) - \
+                    output_distribution.ppf(q_min / 100.0)
+                self.scale_ = self.scale_ / self.adjust_
         else:
             self.scale_ = None
 
@@ -1290,7 +1302,7 @@ def _more_tags(self):
 
 @_deprecate_positional_args
 def robust_scale(X, *, axis=0, with_centering=True, with_scaling=True,
-                 quantile_range=(25.0, 75.0), copy=True):
+                 quantile_range=(25.0, 75.0), unit_variance=False, copy=True):
     """Standardize a dataset along any axis
 
     Center to the median and component wise scale
@@ -1321,6 +1333,10 @@ def robust_scale(X, *, axis=0, with_centering=True, with_scaling=True,
 
         .. versionadded:: 0.18
 
+    unit_variance : boolean, False by default
+        If True, scale data so that normally distributed features have a
+        variance of 1.
+
     copy : boolean, optional, default is True
         set to False to perform inplace row normalization and avoid a
         copy (if the input is already a numpy array or a scipy.sparse
@@ -1357,7 +1373,8 @@ def robust_scale(X, *, axis=0, with_centering=True, with_scaling=True,
         X = X.reshape(X.shape[0], 1)
 
     s = RobustScaler(with_centering=with_centering, with_scaling=with_scaling,
-                     quantile_range=quantile_range, copy=copy)
+                     quantile_range=quantile_range,
+                     unit_variance=unit_variance, copy=copy)
     if axis == 0:
         X = s.fit_transform(X)
     else:
diff --git a/sklearn/preprocessing/tests/test_data.py b/sklearn/preprocessing/tests/test_data.py
index f79703610bee5..bc8bc46bf009e 100644
--- a/sklearn/preprocessing/tests/test_data.py
+++ b/sklearn/preprocessing/tests/test_data.py
@@ -1651,6 +1651,18 @@ def test_robust_scaler_zero_variance_features():
     assert_array_almost_equal(X_trans_new, X_expected_new, decimal=3)
 
 
+def test_robust_scaler_unit_variance():
+    # Check RobustScaler with unit_variance=True is equivalent to
+    # StandardScaler on large normal data
+    rng = np.random.RandomState(0)
+    X = rng.randn(100000, 1) * 5 + 2
+    robust_scaler = RobustScaler(unit_variance=True)
+    standard_scaler = StandardScaler()
+    assert_array_almost_equal(robust_scaler.fit_transform(X),
+                              standard_scaler.fit_transform(X),
+                              decimal=2)
+
+
 def test_maxabs_scaler_zero_variance_features():
     # Check MaxAbsScaler on toy data with zero variance features
     X = [[0., 1., +0.5],


Code review comments:
diff_hunk: @@ -1321,6 +1333,10 @@ def robust_scale(X, *, axis=0, with_centering=True, with_scaling=True,
 
         .. versionadded:: 0.18
 
+    unit_variance : boolean, False by default
+        If True, scale data so that normally distributed features have a
+        variance of 1.
+
path: sklearn/preprocessing/_data.py
text: @user1:
`.. versionadded` and move to be the last arg here as well.

diff_hunk: @@ -1233,6 +1239,12 @@ def fit(self, X, y=None):
 
             self.scale_ = quantiles[1] - quantiles[0]
             self.scale_ = _handle_zeros_in_scale(self.scale_, copy=False)
+            if self.unit_variance:
+                # Create scipy.stats.norm object
+                output_distribution = getattr(stats, 'norm')
path: sklearn/preprocessing/_data.py
text: @user1:
why this instead of `stats.norm`? Is there any difference?

@author:
No i don't think so, updated.

diff_hunk: @@ -1130,6 +1130,10 @@ class RobustScaler(TransformerMixin, BaseEstimator):
 
         .. versionadded:: 0.18
 
+    unit_variance : boolean, False by default
+        If True, scale data so that normally distributed features have a
+        variance of 1.
+
path: sklearn/preprocessing/_data.py
text: @user1:
this needs a `.. versionadded` and needs to be moved to be the last argument. We're still accepting positional args in the next release.

@user1:
Could we also expand on this a bit and explain what happens if the feature is not normally distributed? Since this is the only place we're having the explanation of this parameter (no user guide), it'd be nice to elaborate a bit.

@author:
Hmm in general if `norm.ppf(quantile_range[1]) - norm.ppf(quantile_range[0])` is greater than 1, the output data will be scaled down (decrease variance) and if less than 1 it will be scaled up. Would this be adequate..?