code: 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
  
imports: import numpy
  from scipy import stats
label_desc: Use image node_id 4 for calling a global function with example usage: fit_bivariate_normal(x, y, robust) and returns: list,  sigma1,  sigma2,  alpha