Outlier detection in small samples using pull-clipping

Author: Yannick Copin y.copin@ipnl.in2p3.fr

Abstract: We present and tests different algorithms to compute the mean of very small samples (\(n=4\)) in presence of a small fraction of outliers. The favorite choice appears to be the inverse-variance weighted mean of pull-clipped values, both more robust and statistically efficient (from 10% to 40%) than the median.

%matplotlib notebook

import numpy as N
import scipy.stats as SS

import astropy.stats as AS
import astropy.table as AT

import matplotlib.pyplot as P

from inspec.combine import *

N.random.seed(12345)

Introduction

Consider a sample of independent (possibly-fallacious) measurements \(\{x_i\}\) of the same quantity \(X\), with presumably normal-errors \(\{\sigma_i\}\). One may want to combine these measurements into a single sample mean estimate. Typical user case is the combination of multiple images or spectra in presence of cosmic rays.

In absence of outliers, the Maximum Likelihood Estimate of the mean of measurements \(\{x_i\}\) normally distributed with standard deviation \(\{\sigma_i\}\) is the inverse-variance weighted mean:

\begin{align} \mu_w &= \frac{1}{w}\sum_i w_i x_i \quad\text{with}\quad w_i = \sigma_i^{-2} \quad\text{and}\quad w = \sum_i w_i \\ \sigma_w &= w^{-1/2} \end{align}

Unfortunately, the (weighted) mean is linearly sensitive to outliers: if a value is off by a factor \(\eta\), the mean is off by a factor directly \(\propto \eta\).

In presence of outliers, it is traditional to apply the median statistic, known to be both robust (breakdown point of 50%: up to half of the points can be infinite outliers) and easy to compute:

\begin{align} m &= \text{med}(\{x_i\}) \\ \sigma_m &= \left(\eta_n \langle 1/\sigma_i \rangle_i^2\right)^{-1/2} \\ \eta_n &= \frac{2}{\pi} \times \begin{cases} (n + \pi - 2) & n\;\text{even} \\ (n + \pi/2 - 1) & n\;\text{odd} \end{cases} \end{align}

However, the median has major drawbacks, specially in the case of very small samples:

• The median is known to be significantly less efficient than standard mean, with a variance increased by a factor \(\pi/2 \times n/(n-1)\), ranging from \(\pi/2 \sim 1.6\) for large normal samples (\(n \gg 1\)) to a factor 2 for \(n=4\) samples (in the noiseless case). This typically corresponds to an equivalent waste of effective exposure time needed to reach a given signal-to-noise ratio on the combined measurement.
• The median does not make any standard use of measured errors, treating all measurements as equally accurate (even though a weighted percentile can be defined).
• Furthermore, it is generally the case that there are much less than 50% of outliers in the sample, and the usage of the median is then overdone.

Other robust location estimators (e.g. trimmed/winsorized means) could be of interest in the general case, but their usage is not adapted to very small samples.

Note: We explicitely consider here samples without intrinsic dispersion, i.e. where all the observed scatter of \(\{x_i\}\) can be statistically explained by measurement errors \(\{\sigma_i\}\). In case of non-null instrinsic dispersion \(\sigma\), the proper weighting is \(w_i = (\sigma^2 + \sigma_i^2)^{-1}\).

For very small sample sizes, it therefore appears the optimally-weighted average on an outlier-clipped sample is the best alternative for both a robust and accurate estimate of the sample mean.

Outlier detection

𝜎-clipping (Grubbs’ tests for outliers)

A standard way to detect a single outlier in a normally-distributed sample is the Grubbs’ test (see e.g. NIST/SEMATECH e-Handbook of Statistical Methods):

1. estimate the mean \(\bar{x}\) and standard deviation \(s\) for the sample;

2. compute the z-score for all measurement: \(z_i = (x_i - \bar{x})/s\);

3. the measurement with maximal (resp. minimal) z-score is declared an upper (resp. lower) outlier if

   \[z > z_{\max} = \frac{n-1}{\sqrt{n}}\sqrt{\frac{t_{\alpha/2n,n-2}^2}{n-2+t_{\alpha/2n,n-2}^2}}\]

where \(t_{\alpha/2n,n-2}\) is the critical value of the t-distribution with \((n-2)\) degrees of freedom and a significance level of \(\alpha/(2N)\) (for a two-sided test).

t = N.array([199.31, 199.53, 200.19, 200.82, 201.92, 201.95, 202.18, 245.57])
grubbs_test(t, side=+1, alpha=0.05, modified=False, verbose=True)

H0:  there are no outliers in the data
Ha:  the maximum value is an outlier

Test statistic:      G = [2.46876461121245]
Significance level:  α = 0.05
Critical value for an upper one-tailed test: G_max = 2.031652001549952
Critical region:     Reject H0 if G > G_max

masked_array(data=[False, False, False, False, False, False, False,  True],
             mask=False,
       fill_value=True)

Grubbs’ test for outliers actually closely relates to the intuitive \(\sigma\)-clipping but with a statistically-founded clip value \(z_{\max}\).

ns = N.unique(N.logspace(N.log10(3), 3, 50).astype(int))
fig, ax = P.subplots()
for alpha in (5e-2, 1e-2, 1e-3):
    ax.plot(ns, [ grubbs_gmax(_, alpha=alpha) for _ in ns ], label=u"α = {:.1%}".format(alpha))
ax.set(xlabel='Sample size', xscale='log', ylabel='zmax', title=u"Grubbs' test two-sided critical zmax")
ax.xaxis.set_major_locator(P.matplotlib.ticker.LogLocator(base=10, subs=range(1,10)))
ax.legend(fontsize='small', loc='upper left');

Grubbs’ test can be generalized in two different ways:

• by using robust sample estimates, median and (normalized) Median Absolute Deviation, to compute the modified z-score;
• by using (inverse-variance) weighted sample mean \(\bar{x} = \mu_w\) for location and error estimate \(s = \sigma_i\) for scale (since \(\sigma_i\) is by definition the estimate of the dispersion of \(x_i\)).

Note: Formally, because of potential screening effects (multiple outliers weaken the single-outlier test performance), it is supposedly not appropriate to apply the single-outlier Grubbs’ test sequentially in order to detect multiple outliers. In the following application, we will still apply it sequentially until no more outlier is detected.

Pull-clipping

As traditionnaly defined, Grubbs’ test has two drawbacks:

• all measurements, including a potential outlier, are equally considered in the estimate of the location and scale of the sample needed for the z-score,
• measurement errors \(\{\sigma_i\}\) are not (fully) used in the traditional Grubbs’ test for outliers.

The pull statistic is a seamingly more pertinent quantity than the z-score. For each measurement of the sample, the pull \(p_i\) is defined by:

\[p_i = \frac{x_i - \bar{x}_i}{\sqrt{\sigma_i^2 + \bar{\sigma}_i^2}}\]

where \(\bar{x}_i\) and \(\bar{\sigma}_i\) are the inverse-variance weighted average and its associated error of the sample without point i.

We therefore introduce the pull-clipping, a procedure similar to Grubbs’ test for outliers but using the pull in place of the z-score. We did not try to derive a statistically-founded clipping value \(p_{\max}\), which will have to be estimated from numerical experimentations.

If errors \(\{\sigma_i\}\) are correct and fully represent the observed dispersion of measurements \(\{x_i\}\), the pull have a normal distribution \(\mathcal{N}(0, 1)\).

dt = SS.lognorm(0.25).rvs(size=4)   # dt ~ 1
t = SS.norm(loc=0, scale=dt).rvs()  # t ~ 0 ± 1
t[0] += 10                          # 10-sigma outlier

table = AT.Table([t, dt, zscore(t), zscore(t, modified=True), zscore(t, dt), pull(t, dt)],
                  names=('t', 'dt', 'z-score', 'robust', 'weighted', 'pull'))
#table.pprint(max_width=-1)
table.show_in_notebook()

Table length=4

idx	t	dt	z-score	robust	weighted	pull
0	11.867708894287606	0.9501105669742944	1.4895453801415957	9.928816330449406	9.035352497246803	10.402203491172893
1	1.57064572445153	1.1271990451628713	-0.330674762767773	0.6003816023551078	-1.5192341523114672	-1.6720597058969242
2	0.08159343055294702	0.8782186550018675	-0.5938957350484306	-0.7485978980370558	-3.6454834585085023	-4.318429580845206
3	0.24519993088463823	0.8702867046153093	-0.5649748823253922	-0.6003816023551078	-3.4907175570078923	-4.150447126710755

Implementation and numerical experiments

We develop here a numerical experiment supposed to mimic the spectrum combination cases when some of the pixel are affected by cosmic rays residuals.

Input population

Define a population of m = 10000 samples of n = 4 measurements (e.g. n spectra of m pixels to be combined). Each point have a distribution \(\mathcal{N}(\mu=0, \sigma^2 \sim 1)\), the actual errors being distributed around 1 with a log-normal distribution \(\ln\mathcal{N}(\mu=0, \sigma^2=0.25^2)\).

n = 4                       # Nb of points per sample
m = 10000                   # Nb of samples
size = (m, n)
print("{} realizations of {}-samples".format(m, n))

10000 realizations of 4-samples

shape = 0.25
dx_dist = SS.lognorm(shape)  # Log-normal distribution

Each value has a uniform probability \(\epsilon = 5\%\) of being an outlier:

eps = 0.05
print("Probability of outliers: {:.1%}".format(eps))

outliers = N.random.rand(m, n) <= eps  # Outlier mask

noutliers = outliers.sum(axis=-1)  # Nb of outliers in each sample

rows = [ (i, SS.binom.pmf(i, n, eps), len(noutliers[noutliers == i]), len(noutliers[noutliers == i])/m)
          for i in range(n + 1) ]
table = AT.Table(rows=rows, names=('Outliers', 'Probability', '# of cases', 'Fraction'))
table.show_in_notebook()

Probability of outliers: 5.0%

Table length=5

idx	Outliers	Probability	# of cases	Fraction
0	0	0.81450625	8167	0.8167
1	1	0.17147500000000007	1699	0.1699
2	2	0.013537500000000008	128	0.0128
3	3	0.00047500000000000027	6	0.0006
4	4	6.250000000000004e-06	0	0.0

Outliers follow a log-normal distribution \(\ln\mathcal{N}(\mu=\ln(\mathtt{nsig}=10), \sigma^2=0.25^2)\):

nsig = 10
outliers_dist = SS.lognorm(shape, scale=nsig)

ax = plot_rv(dx_dist, label='StdDev')
ax.set(xscale='log', yscale='log')
ax = plot_rv(outliers_dist, label='Outliers', ax=ax)
ax.legend(fontsize='small')
ax.set_title("Input distributions");

Note: We assume here the outlying process does not significantly affect the measurement variance.

dx = dx_dist.rvs(size=size)
x = N.random.normal(loc=0, scale=dx, size=size)  # Pure gaussian errors
x += N.where(outliers, outliers_dist.rvs(size=size), 0)  # Add outliers

fig = plot_samples(x, dx)
fig.suptitle("Sample distribution", fontsize='large');

fig = plot_samples(x, outliers=outliers)
fig.suptitle("Input outliers", fontsize='large');

Median statistic

med, dmed = sample_median(x, dx)       # Median (unweighted)

Outlier detection (upper one-tailed)

Since the outliers are supposed to mimick additive cosmic rays, we are only looking for upper outliers.

side = +1  # Upper one-tailed outlier test
alpha = 0.01
zmax = grubbs_gmax(n, alpha, one_sided=True)
print("Grubb's critical value for an upper one-tailed test for outliers: zmax={}".format(zmax))

Grubb's critical value for an upper one-tailed test for outliers: zmax=1.492499999999993

𝜎-clipping (standard Grubbs’ test)

identified = outlier_clipping(x, dx, method='sigma', clip=zmax, side=side, verbose=True)
detection_rates(outliers, identified, verbose=True);

Outlier sigma-clipping: upper one-tailed, clip=1.492499999999993
  Iter #1/3: 644 outliers (total: 644)
Input outliers:      1973/40000 (4.93%)
Identified outliers: 644 (TP: 549, FP: 95, FN: 1424)
Sensitivity TPR=1-FNR: 27.83%
Specificity TNR=1-FPR: 99.75%
Precision PPV:         85.25%
Accuracy:              96.20%
Matthews Correlation Coefficient: 0.474

Standard Grubbs’ test has a precision of 84% but a low sensitivity of 28% because of too many false negatives. This method will not be studied further.

nMAD-clipping (robust Grubbs’s test)

identified = outlier_clipping(x, dx, method='robust', clip=zmax, side=side, maxiter=1, verbose=True)
detection_rates(outliers, identified, verbose=True);

Outlier robust-clipping: upper one-tailed, clip=1.492499999999993
  Iter #1/1: 4360 outliers (total: 4360)
Input outliers:      1973/40000 (4.93%)
Identified outliers: 4360 (TP: 1719, FP: 2641, FN: 254)
Sensitivity TPR=1-FNR: 87.13%
Specificity TNR=1-FPR: 93.05%
Precision PPV:         39.43%
Accuracy:              92.76%
Matthews Correlation Coefficient: 0.557

Modified Grubbs’ test has a reasonable accuracy of 93% but a very low precision of 40% because of too many false positives (even when stopping outlier detection after a single iteration). A higher clipping value improves the overall efficiency (as estimated from Matthews Correlation Coefficient), but not to the level of the methods described below:

identified = outlier_clipping(x, dx, method='robust', clip=3, side=side, maxiter=1, verbose=True)
detection_rates(outliers, identified, verbose=True);

Outlier robust-clipping: upper one-tailed, clip=3
  Iter #1/1: 2583 outliers (total: 2583)
Input outliers:      1973/40000 (4.93%)
Identified outliers: 2583 (TP: 1614, FP: 969, FN: 359)
Sensitivity TPR=1-FNR: 81.80%
Specificity TNR=1-FPR: 97.45%
Precision PPV:         62.49%
Accuracy:              96.68%
Matthews Correlation Coefficient: 0.698

Weighted 𝜎-clipping (weighted Grubbs’s test)

In absence of a theoretically justified clipping value, we use here numerically-optimized clip=3 (see Appendix).

clip = 3
identified = outlier_clipping(x, dx, method='weighted', clip=clip, side=side, verbose=True)
detection_rates(outliers, identified, verbose=True)

fig = plot_samples(x, outliers=outliers, identified=identified)
fig.suptitle("Identified outliers (wsigma-clipping)", fontsize='large')

mx = N.ma.masked_array(x, mask=identified)
wm, dwm = sample_weighted_mean(mx, dx)

Outlier weighted-clipping: upper one-tailed, clip=3
  Iter #1/3: 1823 outliers (total: 1823)
  Iter #2/3: 130 outliers (total: 1953)
Input outliers:      1973/40000 (4.93%)
Identified outliers: 1953 (TP: 1946, FP: 7, FN: 27)
Sensitivity TPR=1-FNR: 98.63%
Specificity TNR=1-FPR: 99.98%
Precision PPV:         99.64%
Accuracy:              99.91%
Matthews Correlation Coefficient: 0.991

Pull-clipping

In absence of a theoritically justified clipping value, we use here numerically-optimized clip=3.5 (see Appendix).

clip = 3.5
identified = outlier_clipping(x, dx, method='pull', clip=clip, side=side, verbose=True)
detection_rates(outliers, identified, verbose=True)

fig = plot_samples(x, outliers=outliers, identified=identified)
fig.suptitle("Identified outliers (pull-clipping)", fontsize='large')

mx = N.ma.masked_array(x, mask=identified)
pm, dpm = sample_weighted_mean(mx, dx)

Outlier pull-clipping: upper one-tailed, clip=3.5
  Iter #1/3: 1817 outliers (total: 1817)
  Iter #2/3: 125 outliers (total: 1942)
Input outliers:      1973/40000 (4.93%)
Identified outliers: 1942 (TP: 1937, FP: 5, FN: 36)
Sensitivity TPR=1-FNR: 98.18%
Specificity TNR=1-FPR: 99.99%
Precision PPV:         99.74%
Accuracy:              99.90%
Matthews Correlation Coefficient: 0.989

p = pull(x, dx)
fig = plot_samples(p, outliers=outliers)
fig.suptitle("Pull distribution on original samples", fontsize='large');

mp = pull(x, dx,
          fixed_mean=pm.reshape(-1, 1), fixed_mean_error=dpm.reshape(-1, 1))
fig = plot_samples(mp, outliers=outliers, identified=identified)

fig.suptitle("Pull distribution on identified samples", fontsize='large');

Analysis

We compare now the results from the 3 promising estimators:

• Median: median on input sample
• WSigmaC: weighted mean on (weighted) \(\sigma_w\)-clipped sample
• PullC: weighted mean on pull-clipped sample.

fig, ((axh, axp), (axn, axq)) = P.subplots(2, 2, figsize=(10,10))

_, bins = AS.freedman_bin_width(med, return_bins=True)
axh.hist(med, bins=bins, density=True, log=True, histtype='stepfilled', alpha=0.5,
         label=u"Median: µ={:.3f}, σ={:.3f}".format(med.mean(), med.std(ddof=1)))
axh.hist(wm, bins=bins, density=True, log=True, histtype='stepfilled', alpha=0.5,
         label=u"WSigmaC: µ={:.3f}, σ={:.3f}".format(wm.mean(), wm.std(ddof=1)))
axh.hist(pm, bins=bins, density=True, log=True, histtype='stepfilled', alpha=0.5,
         label=u"PullC: µ={:.3f}, σ={:.3f}".format(pm.mean(), pm.std(ddof=1)))
axh.set(title="Sample mean distribution")
axh.legend(loc='upper right', fontsize='small')

print("\nMedian:")
probplot(med, ax=axn, label='Median')
plot_pull(med, dmed, fixed_mean=0, ax=axp, alpha=0.5, label='Median', log=True, density=True)
probplot(pull(med, dmed, fixed_mean=0), ax=axq, label='Median')
print("\nWeighted sigma-clipping:")
probplot(wm, ax=axn, label='WSigmaC')
plot_pull(wm, dwm, fixed_mean=0, ax=axp, alpha=0.5, label='WSigmaC', log=True, density=True)
probplot(pull(wm, dwm, fixed_mean=0), ax=axq, label='WSigmaC')
print("\nPull-clipping:")
probplot(pm, ax=axn, label='PullC')
plot_pull(pm, dpm, fixed_mean=0, ax=axp, alpha=0.5, label='PullC', log=True, density=True)
probplot(pull(pm, dpm, fixed_mean=0), ax=axq, label='PullC')

axn.set(title="Normal probability plot",
        xlabel="Predicted quantiles", ylabel="Observed quantiles")
axn.legend(loc='upper left', fontsize='small')

plot_rv(SS.norm, ax=axp, label='$\mathcal{N}$(0, 1)', c='k', lw=2, ls='--')
axp.set_title("Pull distribution")
axp.legend(loc='upper right', fontsize='small')

axq.plot([-4, +4], [-4, +4], label='1:1', c='k', lw=2, ls='--')
axq.set(title="Normal probability plot (pull)",
        xlabel="Predicted quantiles")
axq.legend(loc='upper left', fontsize='small');

Median:
Probplot best-fit ax+b: a=0.729777, b=0.13532, R²=0.881908
Pull: mean=+0.234, std=1.430, normality p-value=0
Probplot best-fit ax+b: a=1.25895, b=0.233974, R²=0.879959

Weighted sigma-clipping:
Probplot best-fit ax+b: a=0.516514, b=0.00315546, R²=0.987245
Pull: mean=+0.005, std=1.024, normality p-value=2.14329e-201
Probplot best-fit ax+b: a=1.01755, b=0.00480122, R²=0.993391

Pull-clipping:
Probplot best-fit ax+b: a=0.517475, b=0.00387058, R²=0.988725
Pull: mean=+0.006, std=1.022, normality p-value=3.83791e-92
Probplot best-fit ax+b: a=1.01853, b=0.00588453, R²=0.995928

• The median distribution significantly deviates from the normal distribution starting at predicted quantile +2.1, i.e. for 1 - SS.norm.cdf(2.1) = 1.8% of the samples. This corresponds to the fraction of samples with more than a single outlier, for which the median estimator breaks down. On the other hand, except for very few cases (\(\lesssim 0.03\%\)), the pull-clipped weighted mean is the least sensitive to outliers.
• All three pull distributions are well approximated by the \(\mathcal{N}(0,1)\) distribution (except for the outlying fraction), from which one can conclude that all three sample mean estimators are not significantly biased, and mean errors are well estimated.

plot_comparison(x=med, dx=dmed, xlabel='Median',
                y=pm, dy=dpm, ylabel='Pull-clipping');

rvar = (dmed / dpm.filled())**2  # Variance ratio (converted from masked array)
_, bins = AS.freedman_bin_width(rvar, return_bins=True)

fig, ax = P.subplots()
ax.hist([ rvar[noutliers == i] for i in (0, 1) ], bins=bins,
        density=True, log=True, stacked=True, histtype='stepfilled', alpha=0.5,
        label=[ u"{}: µ={:.2f}".format(i, rvar[noutliers == i].mean()) for i in (0, 1) ])
ax.set(title='Variance ratio Median/PullC', xlabel='Variance ratio')
ax.legend(title="# of outliers", fontsize='small', loc='upper left');

• In absence of outliers, the median is, as expected, significantly less efficient than the weighted mean, by 44% on average.
• Even in presence of a single outlier, the pull-clipped weighted mean still overperforms the median by 8% on average.

Conclusions

We tested three different approachs to compute the mean of a small-size (\(n=4\)) sample in presence of a small fraction (5% probability) of outlying (10 \(\sigma\)-level) values:

• the median of the sample,
• the inverse-variance weighted mean on (weighted) \(\sigma\)-clipped sample,
• the inverse-variance weighted mean on pull-clipped sample.

Among these 3 estimators, the inverse-variance weighted mean on pull-clipped sample appears to have the best performances:

• it is always statistically more efficient than the median, its variance being 30 to 40% smaller than the one of the median in absence of outliers (which is expected to be the vast majority of the cases, 81% in our experiment), and still 10% smaller in presence of a single outlier;
• it is less impacted by multiple outliers, given the very good performance of the pull-clipping process (as measured by a Matthews Correlation Coefficient of 99% in our experiment).

The results that were obtained here are probably slightly dependant of the details of the numerical experiment, notably the distribution adopted for the outliers (fraction, significance, impact on the measurement variance), however, the general conlusion would hold:

• the usage of the significantly sub-efficient median estimator is overdone when samples are only rarely affected by outliers,
• the pull-clipping is a performant procedure to detect outliers.

Appendix

Clip study

We study now the dependency of the outlier detection efficiency (measured in terms of Matthews Correlation Coefficient) on the actual clipping value.

def study_clip(x, dx, clips, method='pull', side=+1):
    """Detection rates as function of clip."""

    rates = N.array([detection_rates(outliers,
                                     outlier_clipping(x, dx, method=method, clip=clip, side=side),
                                     rates=('TPR', 'TNR', 'MCC'), verbose=False)
                     for clip in clips ])

    print(method)
    table = AT.Table([clips, rates[:, 0], rates[:, 1], rates[:, 2]],
                     names='Clip,TPR,TNR,Matthews CC'.split(','))
    print(table)

    fig = P.figure()
    ax = fig.add_subplot(1, 1, 1,
                         xlabel='Clip',
                         ylabel='Fraction or Correlation coefficient',
                         title='{}-clipping'.format(method))
    ax.plot(clips, rates[:, 0], label='Sensitivity (TPR)')
    ax.plot(clips, rates[:, 1], label='Specificity (TNR)')
    ax.plot(clips, rates[:, 2], label='Matthews CC')
    ax.legend(loc='best', fontsize='small')
    ax.set_ylim(0.8, 1.01)

    return ax

clips = N.linspace(2, 6, 9)
study_clip(x, dx, clips, method='weighted');

weighted
Clip        TPR                TNR            Matthews CC
---- ------------------ ------------------ ------------------
 2.0 0.9964521033958439  0.989191889972914 0.9027933521789058
 2.5 0.9959452610238216 0.9976595576826991 0.9748550441978189
 3.0 0.9863152559553978 0.9998159202671786 0.9909070840782379
 3.5 0.9756715661429296 0.9999737028953112 0.9868681293438823
 4.0  0.947795235681703 0.9999737028953112 0.9719578797035292
 4.5 0.9178915357323872                1.0 0.9560323140299808
 5.0  0.863659401926001                1.0 0.9260630959402388
 5.5 0.7901672579827673                1.0  0.884113880294137
 6.0  0.712113532691333                1.0 0.8376356045918112

study_clip(x, dx, clips, method='pull');

pull
Clip        TPR                TNR            Matthews CC
---- ------------------ ------------------ ------------------
 2.0 0.9964521033958439 0.9759118521050832 0.8143276001406025
 2.5 0.9959452610238216 0.9937675861887606   0.93967076450641
 3.0 0.9918905220476432 0.9986851447655614 0.9825886661625981
 3.5 0.9817536746071972 0.9998685144765561 0.9890229015978553
 4.0 0.9706031424227065 0.9999474057906225 0.9839003190634793
 4.5 0.9447541814495691 0.9999737028953112 0.9703200727600765
 5.0 0.9199189052204765                1.0  0.957137665322455
 5.5   0.87886467308667                1.0 0.9345456120797493
 6.0 0.8286872782564623                1.0 0.9063037290967367