code: def statistical_axes(fit, **kw):
      
      method = kw.pop('method', 'noise')
      confidence_level = kw.pop('confidence_level', 0.95)
      dof = kw.pop('dof',2)
  
      nominal = fit.eigenvalues
  
      if method == 'sampling':
          cov = sampling_covariance(fit,**kw)
      elif method == 'noise':
          cov = noise_covariance(fit,**kw)
  
      if kw.pop('chisq', False):
          # Model the incorrect behaviour of using the
          # Chi2 distribution instead of the Fisher
          # distribution (which is a measure of the
          # ratio between the two).
          z = chi2.ppf(confidence_level,dof)
      else:
          z = fisher_statistic(fit.n,confidence_level,dof=dof)
  
      # Apply two fisher F parameters (one along each axis)
      # Since we apply to each axis without division,
      # it is as if we are applying N.sqrt(2*F) to the entire
      # distribution, aligning us with (Francq, 2014)
      err = z*N.sqrt(cov)
      return apply_error_scaling(nominal, err, n=fit.n, **kw)
docs: Hyperbolic error using a statistical process (either sampling
      or noise errors)
  
      Integrates covariance with error level
      and degrees of freedom for plotting
      confidence intervals.
  
      Degrees of freedom is set to 2, which is the
      relevant number of independent dimensions
      to planar fitting of *a priori* centered data.