Repository: scipy/scipy
Pull request: BUG: Fix yeojohnson when transformed data has zero variance
<!-- 
Thanks for contributing a pull request! Please ensure that
your PR satisfies the checklist before submitting:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#checklist-before-submitting-a-pr

Also, please name and describe your PR as you would write a
commit message:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#writing-the-commit-message

Note that we are a team of volunteers; we appreciate your
patience during the review process.

Again, thanks for contributing!
-->

#### Reference issue
<!--Example: Closes gh-WXYZ.-->
Closes https://github.com/scipy/scipy/issues/10821

#### What does this implement/fix?
<!--Please explain your changes.-->
This PR updates `yeojohnson_llf` to handle cases where the variance is zero in the transformed space. Note that `yeojohnson_llf` returns `_XMAX` when the variance is too low, because the formula is `-N/2 * log(variance)`. `_neg_llf` will consider those to be uninteresting by replacing the `_XMAX` with `-_XMAX`.

#### Additional information
The `np.errstate(over='ignore')` is used because `_neg_llf` will return `_XMAX` which causes an overflow in `optimize.brent`.

Problem statement:
Errors with the Yeo-Johnson Transform that also Appear in Scikit-Learn
My issue is about issues with the Yeo-Johnson transform that also appear in scikit-learn

#### Reproducing code example:

```
import numpy as np
import scipy
import sklearn
sklearn.show_versions()
a = np.array([
 3251637.22,620695.44,11642969.00,2223468.22,85307500.00,16494389.89,
  917215.88,11642969.00,2145773.87,4962000.00,620695.44,651234.50,
 1907876.71,4053297.88,3251637.22,3259103.08,9547969.00,20631286.23,
12807072.08,2383819.84,90114500.00,17209575.46,12852969.00,2414609.99,
 2170368.23])
scipy.stats.yeojohnson(a)
```

#### Error message:


```
System:
    python: 3.6.4 |Anaconda, Inc.| (default, Jan 16 2018, 12:04:33)  [GCC 4.2.1 Compatible Clang 4.0.1 (tags/RELEASE_401/final)]
executable: /anaconda3/bin/python
   machine: Darwin-17.7.0-x86_64-i386-64bit

Python deps:
       pip: 19.2.3
setuptools: 38.4.0
   sklearn: 0.22.dev0
     numpy: 1.17.0
     scipy: 1.3.1
    Cython: 0.29.13
    pandas: 0.24.2
matplotlib: 2.1.2
    joblib: 0.13.2
/anaconda3/lib/python3.6/site-packages/scipy/stats/morestats.py:1465: RuntimeWarning: divide by zero encountered in log
  loglike = -n_samples / 2 * np.log(trans.var(axis=0))
/anaconda3/lib/python3.6/site-packages/scipy/optimize/optimize.py:1961: RuntimeWarning: invalid value encountered in double_scalars
  p = (x - v) * tmp2 - (x - w) * tmp1
/anaconda3/lib/python3.6/site-packages/scipy/optimize/optimize.py:1962: RuntimeWarning: invalid value encountered in double_scalars
  tmp2 = 2.0 * (tmp2 - tmp1)
/anaconda3/lib/python3.6/site-packages/scipy/optimize/optimize.py:1960: RuntimeWarning: invalid value encountered in double_scalars
  tmp2 = (x - v) * (fx - fw)
/anaconda3/lib/python3.6/site-packages/scipy/optimize/optimize.py:1959: RuntimeWarning: invalid value encountered in double_scalars
  tmp1 = (x - w) * (fx - fv)
Out[7]:
(array([0.23389584, 0.23389584, 0.23389584, 0.23389584, 0.23389584,
        0.23389584, 0.23389584, 0.23389584, 0.23389584, 0.23389584,
        0.23389584, 0.23389584, 0.23389584, 0.23389584, 0.23389584,
        0.23389584, 0.23389584, 0.23389584, 0.23389584, 0.23389584,
        0.23389584, 0.23389584, 0.23389584, 0.23389584, 0.23389584]),
 -4.275407328878737)
```

#### Scipy/Numpy/Python version information:

Please see above

Hints from issue discussion:
The issue is that I supply an array with large variance and am returned an array with 0 variance plus many error messages.

I have reported a very similar error to scikit-learn, from whence I believe the code originated:

https://github.com/scikit-learn/scikit-learn/issues/14959

There is also a comment in a scikit-learn pull request that bears on this:

https://github.com/scikit-learn/scikit-learn/pull/12522

    "There's a numerical instability issue in scipy's box-cox (scipy/scipy#6873) that is reproducible for yeo-johnson as well. I'm keeping an eye on this and will port the fix once it's merged on scipy."
The trick that is used in https://github.com/scipy/scipy/pull/10072 (fox BoxCox)  is that:

1. The loglikelihood only depends on the transformed values through their variance
2. Var(X) = Var(X + constant)
3. The numerical instability is caused by a constant offset = `1 / λ` in the transformed values. We can just remove it: the variance stays the same and the loglikelihood is still correct.

While that offset is constant for boxcox, in the yeo johnson transform the offset is different depending on the sign of the input (see [original paper](https://scholar.google.com/scholar?hl=en&as_sdt=0%2C44&q=https%3A%2F%2Fwww.jstor.org%2Fstable%2Fpdf%2F2673623.pdf&btnG=): it's either `1 / λ` or `1/(2-λ)`

Since we wouldn't be offsetting the data with a constant offset, I don't think we can use the same trick here :/

@WarrenWeckesser any ideas?
@NicolasHug, I tried the same thing when I saw first saw this issue, but discovered like you did that the trick doesn't work for Yeo-Johnson.  I don't know of any other neat ways to avoid the problem.

The standard advice for a dataset such as the one given above is to rescale the data to avoid those huge values.  Also note that for such a dataset of positive values `y`, the Yeo-Johnson transform is the same as the Box-Cox transform of `y + 1`.  So the practical advice I would give for handling such a dataset is to just use Box-Cox.

(Edited to give a more definitive statement about the relation between Yeo-Johnson and Box-Cox for a positive dataset.)
> ...just use Box-Cox

Here's a demonstration that shows Box-Cox works and is basically the same as Yeo-Johnson for this dataset:

```
In [51]: from scipy.stats import yeojohnson, boxcox                             

In [52]: a = np.array([ 
    ...:  3251637.22,620695.44,11642969.00,2223468.22,85307500.00,16494389.89, 
    ...:   917215.88,11642969.00,2145773.87,4962000.00,620695.44,651234.50, 
    ...:  1907876.71,4053297.88,3251637.22,3259103.08,9547969.00,20631286.23, 
    ...: 12807072.08,2383819.84,90114500.00,17209575.46,12852969.00,2414609.99, 
    ...:  2170368.23])                                                          
```
Compute the Box-Cox transform of `a + 1`:
```
In [106]: t, lam = boxcox(a + 1)

In [107]: t
Out[107]: 
array([6.65421632, 6.38491716, 6.82593329, 6.59736951, 7.04368037,
       6.86814054, 6.45366817, 6.82593329, 6.59190665, 6.71425592,
       6.38491716, 6.3935601 , 6.57368014, 6.68592749, 6.65421632,
       6.65455097, 6.80103851, 6.89427917, 6.83766875, 6.60799291,
       7.04892209, 6.87315608, 6.8381065 , 6.60994041, 6.59365976])

In [108]: lam
Out[108]: -0.1283482601400123

```
Use `lam` as the parameter in the Yeo-Johnson transform.  Note that the transformed values agree with the Box-Cox output.
```
In [109]: yeojohnson(a, lam)
Out[109]: 
array([6.65421632, 6.38491716, 6.82593329, 6.59736951, 7.04368037,
       6.86814054, 6.45366817, 6.82593329, 6.59190665, 6.71425592,
       6.38491716, 6.3935601 , 6.57368014, 6.68592749, 6.65421632,
       6.65455097, 6.80103851, 6.89427917, 6.83766875, 6.60799291,
       7.04892209, 6.87315608, 6.8381065 , 6.60994041, 6.59365976])

```
@thomasjpfan would you be willing to submit the same fix as in scikit-learn/scikit-learn#20653 over here? Thanks for considering it! (Curious if you know why scikit-learn and scipy don't just share an implementation for this one?)

Commit message: BUG: Fix yeojohnson when transformed data has zero variance

Patch to review:
diff --git a/scipy/stats/_morestats.py b/scipy/stats/_morestats.py
index 06978177e20d..dd2045ca9003 100644
--- a/scipy/stats/_morestats.py
+++ b/scipy/stats/_morestats.py
@@ -18,6 +18,7 @@
 from . import _stats_py
 from ._stats_py import find_repeats, _contains_nan, _normtest_finish
 from .contingency import chi2_contingency
+from ._constants import _XMIN, _XMAX
 from . import distributions
 from ._distn_infrastructure import rv_generic
 from ._hypotests import _get_wilcoxon_distr
@@ -1525,6 +1526,9 @@ def yeojohnson_llf(lmb, data):
     where :math:`\hat{\sigma}^2` is estimated variance of the the Yeo-Johnson
     transformed input data ``x``.
 
+    If the Yeo-Johnson transformed input data has zero variance, then the
+    log-likelihood is set to the largest representable float.
+
     .. versionadded:: 1.2.0
 
     Examples
@@ -1579,10 +1583,17 @@ def yeojohnson_llf(lmb, data):
         return np.nan
 
     trans = _yeojohnson_transform(data, lmb)
+    trans_var = trans.var(axis=0)
+    loglike = np.empty_like(trans_var)
 
-    loglike = -n_samples / 2 * np.log(trans.var(axis=0))
-    loglike += (lmb - 1) * (np.sign(data) * np.log(np.abs(data) + 1)).sum(axis=0)
+    # Avoid RuntimeWarning raised by np.log when the variance is too low
+    tiny_variance = trans_var < _XMIN
+    loglike[tiny_variance] = _XMAX
 
+    loglike[~tiny_variance] = (
+        -n_samples / 2 * np.log(trans_var[~tiny_variance]))
+    loglike[~tiny_variance] += (
+        (lmb - 1) * (np.sign(data) * np.log(np.abs(data) + 1)).sum(axis=0))
     return loglike
 
 
@@ -1634,9 +1645,14 @@ def yeojohnson_normmax(x, brack=(-2, 2)):
 
     """
     def _neg_llf(lmbda, data):
-        return -yeojohnson_llf(lmbda, data)
-
-    return optimize.brent(_neg_llf, brack=brack, args=(x,))
+        llf = yeojohnson_llf(lmbda, data)
+        # reject likelihoods that are _XMAX which are likely due to small
+        # variance in the transformed space
+        llf[llf == _XMAX] = -_XMAX
+        return -llf
+
+    with np.errstate(over='ignore'):
+        return optimize.brent(_neg_llf, brack=brack, args=(x,))
 
 
 def yeojohnson_normplot(x, la, lb, plot=None, N=80):
diff --git a/scipy/stats/tests/test_morestats.py b/scipy/stats/tests/test_morestats.py
index 903a1e6201fa..e09087002e01 100644
--- a/scipy/stats/tests/test_morestats.py
+++ b/scipy/stats/tests/test_morestats.py
@@ -2103,6 +2103,19 @@ def test_input_dtype_integer(self, dtype):
         assert_allclose(xt_int, xt_float, rtol=1e-7)
         assert_allclose(lmbda_int, lmbda_float, rtol=1e-7)
 
+    def test_input_high_variance(self):
+        # non-regression test for gh-10821
+        x = np.array([3251637.22, 620695.44, 11642969.00, 2223468.22,
+                      85307500.00, 16494389.89, 917215.88, 11642969.00,
+                      2145773.87, 4962000.00, 620695.44, 651234.50,
+                      1907876.71, 4053297.88, 3251637.22, 3259103.08,
+                      9547969.00, 20631286.23, 12807072.08, 2383819.84,
+                      90114500.00, 17209575.46, 12852969.00, 2414609.99,
+                      2170368.23])
+        xt_yeo, lam_yeo = stats.yeojohnson(x)
+        xt_box = stats.boxcox(x + 1, lmbda=lam_yeo)
+        assert_allclose(xt_yeo, xt_box)
+
 
 class TestYeojohnsonNormmax:
     def setup_method(self):


Code review comments:
diff_hunk: @@ -2103,6 +2103,19 @@ def test_input_dtype_integer(self, dtype):
         assert_allclose(xt_int, xt_float, rtol=1e-7)
         assert_allclose(lmbda_int, lmbda_float, rtol=1e-7)
 
+    def test_input_high_variance(self):
+        # non-regression test for gh-10821
+        x = np.array([3251637.22, 620695.44, 11642969.00, 2223468.22,
+                      85307500.00, 16494389.89, 917215.88, 11642969.00,
+                      2145773.87, 4962000.00, 620695.44, 651234.50,
+                      1907876.71, 4053297.88, 3251637.22, 3259103.08,
+                      9547969.00, 20631286.23, 12807072.08, 2383819.84,
+                      90114500.00, 17209575.46, 12852969.00, 2414609.99,
+                      2170368.23])
+        xt_yeo, lam_yeo = stats.yeojohnson(x)
+        xt_box = stats.boxcox(x + 1, lmbda=lam_yeo)
+        assert_allclose(xt_yeo, xt_box)
path: scipy/stats/tests/test_morestats.py
text: @user1:
Oops, one more thing. This test would pass even in `main` because the relationship is true in general. I think what we want is:
```suggestion
        xt_yeo, lam_yeo = stats.yeojohnson(x)
        xt_box, lam_box = stats.boxcox(x + 1)
        assert_allclose(lam_yeo, lam_box, rtol=1e-6)
```
The tolerance is because the relative error is `1.46440174e-07` on my machine, slightly over the default tolerance.