def statistical_axes(fit, **kw):
    """
    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.
    """
    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)