title: Python Library: Yeo-Johnson transformation
body: <p>I'm looking for an implementation of the <a href="https://www.stat.umn.edu/arc/yjpower.pdf" rel="nofollow">Yeo-Johnson power transformation</a> in Python. Unfortunately, SciPy library has the implementation of box-cox (applies only to positive numbers)!</p>
  
  <p>At the moment, I'm using an R package (<a href="http://www.inside-r.org/packages/cran/vgam/docs/yeo.johnson" rel="nofollow">caret::preProcess</a>) via RPY2; however, it has a big performance hit.</p>
  
  <p>Thanks</p>
  
  <p>P.S.: calling caret::preProcess from Python:</p>
  
  <pre><code>#!/usr/bin/env python
  from rpy2.robjects.functions import SignatureTranslatedFunction
  import rpy2.robjects as robjects
  from rpy2.robjects.packages import importr
  
  r_caret = importr("caret")
  r_stats = importr("stats")
  
  feature_in = list(range(-5, 5))
  
  r_caret.preProcess = SignatureTranslatedFunction(r_caret.preProcess, {'r_lambda': 'lambda'})
  r_scale = r_caret.preProcess(robjects.DataFrame({'a': robjects.FloatVector(feature_in)}),
                            method=robjects.StrVector(["center", "scale", "YeoJohnson"]),
                            r_lambda=robjects.IntVector([-10, 10, 0.1]),
                            na_remove=False)
  temp = r_stats.predict(r_scale, robjects.DataFrame({'a': robjects.FloatVector(feature_in)}))
  feature_out = np.array(temp.rx2(1), dtype=np.float64)
  
  print(feature_in)  # output [-5, -4, -3, -2, -1, 0, 1, 2, 3, 4]
  print(feature_out)  # output [-1.46625361 -1.14751589 -0.82725836 -0.50508336 -0.1803113  0.14850228  0.48337191  0.8224039  1.16416754  1.5079769]
  </code></pre>
  
  <p>The solution is as the following:</p>
  
  <pre><code>#!/usr/bin/env python
  # -*- coding: UTF-8 -*-
  import warnings
  import numpy as np
  import pandas as pd
  import sys
  
  __author__ = "Mohsen Mesgarpour"
  __copyright__ = "Copyright 2016, https://github.com/mesgarpour"
  __credits__ = ["Mohsen Mesgarpour"]
  __license__ = "GPL"
  __version__ = "1.0"
  __maintainer__ = "Mohsen Mesgarpour"
  __email__ = "mohsen.mesgarpour@gmail.com"
  
  
  class YeoJohnson:
      """
      Computing Yeo-Johnson transofrmation, which is an extension of Box-Cox transformation
      but can handle both positive and negative values.
  
      References:
      Weisberg, S. (2001). Yeo-Johnson Power Transformations.
      Department of Applied Statistics, University of Minnesota. Retrieved June, 1, 2003.
      https://www.stat.umn.edu/arc/yjpower.pdf
  
      Adapted from CRAN - Package VGAM
      """
      def fit(self, y, lmbda, derivative=0, epsilon=np.finfo(np.float).eps, inverse=False):
          """
          :param y: The variable to be transformed (numeric array).
          :param lmbda: The function's Lambda value (numeric value or array).
          :param derivative: The derivative with respect to lambda
          (non-negative integer; default: ordinary function evaluation).
          :param epsilon: The lambda's tolerance (positive value).
          :param inverse: The inverse transformation option (logical value).
          :return: The Yeo-Johnson transformation or its inverse, or its derivatives with respect to lambda, of y.
          """
          # Validate arguments
          self.__validate(y, lmbda, derivative, epsilon, inverse)
  
          # initialise
          y = np.array(y, dtype=float)
          result = y
          if not (isinstance(lmbda, list) or isinstance(lmbda, np.ndarray)):
              lmbda, y = np.broadcast_arrays(lmbda, y)
              lmbda = np.array(lmbda, dtype=float)
          l0 = np.abs(lmbda) &gt; epsilon
          l2 = np.abs(lmbda - 2) &gt; epsilon
  
          # Inverse
          with warnings.catch_warnings():  # suppress warnings
              warnings.simplefilter("ignore")
              if inverse is True:
                  mask = np.where(((y &gt;= 0) &amp; l0) == True)
                  result[mask] = np.power(np.multiply(y[mask], lmbda[mask]) + 1, 1 / lmbda[mask]) - 1
  
                  mask = np.where(((y &gt;= 0) &amp; ~l0) == True)
                  result[mask] = np.expm1(y[mask])
  
                  mask = np.where(((y &lt; 0) &amp; l2) == True)
                  result[mask] = 1 - np.power(np.multiply(-(2 - lmbda[mask]), y[mask]) + 1, 1 / (2 - lmbda[mask]))
  
                  mask = np.where(((y &lt; 0) &amp; ~l2) == True)
                  result[mask] = -np.expm1(-y[mask])
  
              # Derivative
              else:
                  if derivative == 0:
                      mask = np.where(((y &gt;= 0) &amp; l0) == True)
                      result[mask] = np.divide(np.power(y[mask] + 1, lmbda[mask]) - 1, lmbda[mask])
  
                      mask = np.where(((y &gt;= 0) &amp; ~l0) == True)
                      result[mask] = np.log1p(y[mask])
  
                      mask = np.where(((y &lt; 0) &amp; l2) == True)
                      result[mask] = np.divide(-(np.power(-y[mask] + 1, 2 - lmbda[mask]) - 1), 2 - lmbda[mask])
  
                      mask = np.where(((y &lt; 0) &amp; ~l2) == True)
                      result[mask] = -np.log1p(-y[mask])
  
                  # Not Derivative
                  else:
                      p = self.fit(y, lmbda, derivative=derivative - 1, epsilon=epsilon, inverse=inverse)
  
                      mask = np.where(((y &gt;= 0) &amp; l0) == True)
                      result[mask] = np.divide(np.multiply(np.power(y[mask] + 1, lmbda[mask]),
                                                           np.power(np.log1p(y[mask]), derivative)) -
                                               np.multiply(derivative, p[mask]), lmbda[mask])
  
                      mask = np.where(((y &gt;= 0) &amp; ~l0) == True)
                      result[mask] = np.divide(np.power(np.log1p(y[mask]), derivative + 1), derivative + 1)
  
                      mask = np.where(((y &lt; 0) &amp; l2) == True)
                      result[mask] = np.divide(-(np.multiply(np.power(-y[mask] + 1, 2 - lmbda[mask]),
                                                                      np.power(-np.log1p(-y[mask]), derivative)) -
                                                          np.multiply(derivative, p[mask])), 2 - lmbda[mask])
  
                      mask = np.where(((y &lt; 0) &amp; ~l2) == True)
                      result[mask] = np.divide(np.power(-np.log1p(-y[mask]), derivative + 1), derivative + 1)
  
          return result
  
      @staticmethod
      def __validate(y, lmbda, derivative, epsilon, inverse):
          try:
              if not isinstance(y, (list, np.ndarray, pd.Series)):
                  raise Exception("Argument 'y' must be a list!")
              if not isinstance(lmbda, (int, float, np.int, np.float)):
                  if not isinstance(lmbda, (list, np.ndarray, pd.Series)) or len(lmbda) != len(y):
                      raise Exception("Argument 'lmbda' must be a number "
                                      "or a list, which its length matches 'y' argument!")
              if not isinstance(derivative, (int, float, np.int, np.float)) or derivative &lt; 0:
                  raise Exception("Argument 'derivative' must be a non-negative integer!")
              if not isinstance(epsilon, (int, float, np.int, np.float)) or epsilon &lt;= 0:
                  raise Exception("Argument 'epsilon' must be a positive number!")
              if not isinstance(inverse, bool):
                  raise Exception("Argument 'inverse' must be boolean!")
              if inverse is True and derivative != 0:
                  raise Exception("Argument 'derivative' must be zero "
                                  "when argument 'inverse' is 'True'!")
          except ():
              sys.exit()
  </code></pre>