code: """
  Notes:
      These functions were used with modification from the astroML python
      function `fit_bivariate_normal` which can be found [here]
      (http://www.astroml.org/book_figures/chapter3/fig_robust_pca.html). In
      this modified version, the percentiles are computed using the numpy
      function `nanpercentile`.
  License:
      These functions were borrowed from the AstroML library in agreement with
      their 3-clause BSD license. These functions retain this 3-clause BSD
      licensing and are the copyright of the owners.
  """
  import numpy as np
  from scipy import stats
  
  #from scipy.special import erfinv
  #sigmaG_factor = 1. / (2 * np.sqrt(2) * erfinv(0.5))
  sigmaG_factor = 0.74130110925280102
  
  
  def mean_sigma(a, axis=None, dtype=None, ddof=0, keepdims=False):
      """
      Compute mean and standard deviation for an array
  
      Parameters
      ----------
      a : array_like
          Array containing numbers whose mean is desired. If `a` is not an
          array, a conversion is attempted.
      axis : int, optional
          Axis along which the means are computed. The default is to compute
          the mean of the flattened array.
      dtype : dtype, optional
          Type to use in computing the standard deviation. For arrays of
          integer type the default is float64, for arrays of float types it is
          the same as the array type.
      keepdims : bool, optional
          If this is set to True, the axes which are reduced are left
          in the result as dimensions with size one. With this option,
          the result will broadcast correctly against the original `arr`.
  
      Returns
      -------
      mu : ndarray, see dtype parameter above
          array containing the mean values
      sigma : ndarray, see dtype parameter above.
          array containing the standard deviation
  
      See Also
      --------
      median_sigmaG : robust rank-based version of this calculation.
      Notes
  
      -----
      This routine simply calls ``np.mean`` and ``np.std``, passing the
      keyword arguments to them.  It is provided for ease of comparison
      with the function median_sigmaG()
      """
      mu = np.mean(a, axis=axis, dtype=dtype)
      sigma = np.std(a, axis=axis, dtype=dtype, ddof=ddof)
  
      if keepdims:
          if axis is None:
              newshape = a.ndim * (1,)
          else:
              newshape = np.asarray(a.shape)
              newshape[axis] = 1
  
          mu = mu.reshape(newshape)
          sigma = sigma.reshape(newshape)
  
      return mu, sigma
  
  
  def median_sigmaG(a, axis=None, overwrite_input=False, keepdims=False):
      """
      Compute median and rank-based estimate of the standard deviation
  
      Parameters
      ----------
      a : array_like
          Array containing numbers whose mean is desired. If `a` is not an
          array, a conversion is attempted.
      axis : int, optional
          Axis along which the means are computed. The default is to compute
          the mean of the flattened array.
      overwrite_input : bool, optional
         If True, then allow use of memory of input array `a` for
         calculations. The input array will be modified by the call to
         median. This will save memory when you do not need to preserve
         the contents of the input array. Treat the input as undefined,
         but it will probably be fully or partially sorted.
         Default is False. Note that, if `overwrite_input` is True and the
         input is not already an array, an error will be raised.
      keepdims : bool, optional
          If this is set to True, the axes which are reduced are left
          in the result as dimensions with size one. With this option,
          the result will broadcast correctly against the original `arr`.
  
      Returns
      -------
      median : ndarray, see dtype parameter above
          array containing the median values
      sigmaG : ndarray, see dtype parameter above.
          array containing the robust estimator of the standard deviation
  
      See Also
      --------
      mean_sigma : non-robust version of this calculation
      sigmaG : robust rank-based estimate of standard deviation
  
      Notes
      -----
      This routine uses a single call to ``np.nanpercentile`` to find the
      quartiles along the given axis, and uses these to compute the
      median and sigmaG:
      median = q50
      sigmaG = (q75 - q25) * 0.7413
      where 0.7413 ~ 1 / (2 sqrt(2) erf^-1(0.5))
      """
      q25, median, q75 = np.nanpercentile(a, [25, 50, 75],
                                          axis=axis,
                                          overwrite_input=overwrite_input)
      sigmaG = sigmaG_factor * (q75 - q25)
  
      if keepdims:
          if axis is None:
              newshape = a.ndim * (1,)
          else:
              newshape = np.asarray(a.shape)
              newshape[axis] = 1
  
          median = median.reshape(newshape)
          sigmaG = sigmaG.reshape(newshape)
  
      return median, sigmaG
  
  
  def sigmaG(a, axis=None, overwrite_input=False, keepdims=False):
      """
      Compute the rank-based estimate of the standard deviation
  
      Parameters
      ----------
      a : array_like
          Array containing numbers whose mean is desired. If `a` is not an
          array, a conversion is attempted.
      axis : int, optional
          Axis along which the means are computed. The default is to compute
          the mean of the flattened array.
      overwrite_input : bool, optional
         If True, then allow use of memory of input array `a` for
         calculations. The input array will be modified by the call to
         median. This will save memory when you do not need to preserve
         the contents of the input array. Treat the input as undefined,
         but it will probably be fully or partially sorted.
         Default is False. Note that, if `overwrite_input` is True and the
         input is not already an array, an error will be raised.
      keepdims : bool, optional
          If this is set to True, the axes which are reduced are left
          in the result as dimensions with size one. With this option,
          the result will broadcast correctly against the original `arr`.
  
      Returns
      -------
      median : ndarray, see dtype parameter above
          array containing the median values
      sigmaG : ndarray, see dtype parameter above.
          array containing the robust estimator of the standard deviation
  
      See Also
      --------
      median_sigmaG : robust rank-based estimate of mean and standard deviation
  
      Notes
      -----
      This routine uses a single call to ``np.nanpercentile`` to find the
      quartiles along the given axis, and uses these to compute the
      sigmaG, a robust estimate of the standard deviation sigma:
      sigmaG = 0.7413 * (q75 - q25)
      where 0.7413 ~ 1 / (2 sqrt(2) erf^-1(0.5))
      """
      q25, q75 = np.nanpercentile(a, [25, 75],
                                  axis=axis,
                                  overwrite_input=overwrite_input)
      sigmaG = sigmaG_factor * (q75 - q25)
  
      if keepdims:
          if axis is None:
              newshape = a.ndim * (1,)
          else:
              newshape = np.asarray(a.shape)
              newshape[axis] = 1
  
          sigmaG = sigmaG.reshape(newshape)
  
      return sigmaG
  
  
  def fit_bivariate_normal(x, y, robust=False):
      """
      Fit bivariate normal parameters to a 2D distribution of points
  
      Parameters
      ----------
      x, y : array_like
          The x, y coordinates of the points
      robust : boolean (optional, default=False)
          If True, then use rank-based statistics which are robust to outliers
          Otherwise, use mean/std statistics which are not robust
  
      Returns
      -------
      mu : tuple
          (x, y) location of the best-fit bivariate normal
      sigma_1, sigma_2 : float
          The best-fit gaussian widths in the uncorrelated frame
      alpha : float
          The rotation angle in radians of the uncorrelated frame
      """
      x = np.asarray(x)
      y = np.asarray(y)
  
      assert x.shape == y.shape
  
      if robust:
          # use quartiles to compute center and spread
          med_x, sigmaG_x = median_sigmaG(x)
          med_y, sigmaG_y = median_sigmaG(y)
  
          # define the principal variables from Shevlyakov & Smirnov (2011)
          sx = 2 * sigmaG_x
          sy = 2 * sigmaG_y
  
          u = (x / sx + y / sy) / np.sqrt(2)
          v = (x / sx - y / sy) / np.sqrt(2)
  
          med_u, sigmaG_u = median_sigmaG(u)
          med_v, sigmaG_v = median_sigmaG(v)
  
          r_xy = ((sigmaG_u ** 2 - sigmaG_v ** 2) /
                  (sigmaG_u ** 2 + sigmaG_v ** 2))
  
          # rename estimators
          mu_x, mu_y = med_x, med_y
          sigma_x, sigma_y = sigmaG_x, sigmaG_y
      else:
          mu_x = np.mean(x)
          sigma_x = np.std(x)
  
          mu_y = np.mean(y)
          sigma_y = np.std(y)
  
          r_xy = stats.pearsonr(x, y)[0]
  
      # We need to use the full (-180, 180) version of arctan: this is
      # np.arctan2(x, y) = np.arctan(x / y), modulo 180 degrees
      sigma_xy = r_xy * sigma_x * sigma_y
      alpha = 0.5 * np.arctan2(2 * sigma_xy, sigma_x ** 2 - sigma_y ** 2)
  
      sigma1 = np.sqrt((0.5 * (sigma_x ** 2 + sigma_y ** 2) +
                        np.sqrt(0.25 * (sigma_x ** 2 - sigma_y ** 2) ** 2 +
                                sigma_xy ** 2)))
      sigma2 = np.sqrt((0.5 * (sigma_x ** 2 + sigma_y ** 2) -
                        np.sqrt(0.25 * (sigma_x ** 2 - sigma_y ** 2) ** 2 +
                                sigma_xy ** 2)))
  
      return [mu_x, mu_y], sigma1, sigma2, alpha
  
apis: numpy.arctan2 scipy.stats.pearsonr numpy.asarray numpy.nanpercentile numpy.sqrt numpy.std numpy.mean