returns confidence interval (dev_min, dev_max) 
        for a given deviation dev, equivalent degrees of freedom edf,
        and degree of confidence ci.
        
    Parameters
    ----------
    dev: float
        Mean value (e.g. adev) around which we produce the confidence interval
    edf: float
        Equivalent degrees of freedon
    ci: float, defaults to scipy.special.erf(1/math.sqrt(2))
        for 1-sigma standard error set
        ci = scipy.special.erf(1/math.sqrt(2))
            = 0.68268949213708585
        
    Returns
    -------
    (dev_min, dev_max): (float, float)
        Confidence interval

def confidence_interval(dev, edf, ci=ONE_SIGMA_CI):
    """ returns confidence interval (dev_min, dev_max) 
        for a given deviation dev, equivalent degrees of freedom edf,
        and degree of confidence ci.
        
    Parameters
    ----------
    dev: float
        Mean value (e.g. adev) around which we produce the confidence interval
    edf: float
        Equivalent degrees of freedon
    ci: float, defaults to scipy.special.erf(1/math.sqrt(2))
        for 1-sigma standard error set
        ci = scipy.special.erf(1/math.sqrt(2))
            = 0.68268949213708585
        
    Returns
    -------
    (dev_min, dev_max): (float, float)
        Confidence interval
    """
    ci_l = min(np.abs(ci), np.abs((ci-1))) / 2
    ci_h = 1 - ci_l

    # function from scipy, works OK, but scipy is large and slow to build
    chi2_l = scipy.stats.chi2.ppf(ci_l, edf)
    chi2_h = scipy.stats.chi2.ppf(ci_h, edf)

    variance = dev*dev
    var_l = float(edf) * variance / chi2_h  # NIST SP1065 eqn (45)
    var_h = float(edf) * variance / chi2_l
    return (np.sqrt(var_l), np.sqrt(var_h))