Outlier Detection: Modified Thompson Test¶
example of outlier test using modified Thompson technique¶
from pylab import arange, argmax, array, np, plot, sqrt, xlabel, xlim, ylabel
import matplotlib as mpl
mpl.rcParams['lines.linewidth']=2
mpl.rcParams['lines.color']='r'
mpl.rcParams['figure.figsize']=(8,6)
mpl.rcParams['font.size']=14
mpl.rcParams['axes.labelsize']=20
x = array([48.9, 49.2, 49.2, 49.3, 49.3, 49.8, 49.9, 50.1, 50.2, 50.7])
plot(arange(1,11),x,'o'),xlim([0,11]),xlabel('$n$'),ylabel('$T\; [^\circ C \,]$');
is there an outlier?¶
Sort, look at the lowest and largest and plot to visualize
plot(arange(2,10),x[1:-1],'o'), plot(1,x[0],'rs',markersize=10,linewidth=2),plot(10,x[-1],'gs',markersize=10)
xlim([0,11])
xlabel('$\hat{n}$'),ylabel('Sorted $T\; [^\circ C \,]$');
The sample mean and standard deviation, \(\bar{x}\), \(S_x\)¶
mean = 49.66, std = 0.57
Defivations for each suspicious point, take the largest \(\delta_i = |x_i - \bar{x}|\)¶
delta = abs(x - meanx)
print ('suspicious points are first and last:')
print ('%4.3f %4.3f' % (delta[0],delta[-1]))
print ('suspicious point is:', argmax(delta), 'deviation is = %4.3f' % max(delta))
suspicious points are first and last: 0.760 1.040 suspicious point is: 9 deviation is = 1.040
Modified Thompson technique, \(\tau\)¶
Define the confidence level (95%), therefore \(\alpha = 5\%\). Because we deal with outliers, the DOF is smaller:
\[ \tau = \frac{t_{\alpha/2} (n-1)}{\sqrt{n} \sqrt{n-2+t_{\alpha/2}^2}}, \qquad \alpha = 0.05, \quad df = n - 2 \]
if:
\[\delta_i > \tau S\]
then the point is an outlier
remove it, estimate new \(\bar{x}, S, \delta_i\) and repeat the test. until there is no outlier in the set.
from scipy.stats import t
_n = len(x)
_tv = t.isf(0.05 / 2, _n - 2)
_tau = _tv * (_n - 1) / (sqrt(_n) * sqrt(_n - 2 + _tv ** 2))
print('n = %d, t = %6.4f, tau = %6.4f' % (_n, _tv, _tau))
print('compare: %6.3f vs. %6.3f ' % (max(delta), _tau * stdx))
print('Is max() above $t_{\\nu,95}S$? %s ' % (max(delta) > _tau * stdx))
n = 10, t = 2.3060, tau = 1.7984
compare: 1.040 vs. 1.022
Is max() above $t_{\nu,95}S$? True
So we remove the outlier and repeat the test (from the beginning)¶
xnew = x[:-1].copy()
meanx_1 = np.mean(xnew)
stdx_1 = np.std(xnew, ddof=1)
print('x: ', xnew)
print('mean = %6.2f, std = %6.2f\n' % (meanx_1, stdx_1))
delta_1 = abs(xnew - meanx_1)
print('deviations: ', delta_1)
print('\n')
print('suspicious point is: %f, its deviation is = %f \n' % (argmax(delta_1), max(delta_1)))
_n = len(xnew)
_tv = t.isf(0.05 / 2, _n - 2)
_tau = _tv * (_n - 1) / (sqrt(_n) * sqrt(_n - 2 + _tv ** 2))
print('n = %d, t = %6.4f, tau = %6.4f\n' % (_n, _tv, _tau))
print('compare: %6.3f vs. %6.3f \n' % (max(delta_1), _tau * stdx_1))
print('Is it outlier? :', max(delta_1) > _tau * stdx_1)
x: [48.9 49.2 49.2 49.3 49.3 49.8 49.9 50.1 50.2] mean = 49.54, std = 0.46 deviations: [0.64444444 0.34444444 0.34444444 0.24444444 0.24444444 0.25555556 0.35555556 0.55555556 0.65555556] suspicious point is: 8.000000, its deviation is = 0.655556 n = 9, t = 2.3646, tau = 1.7770 compare: 0.656 vs. 0.820 Is it outlier? : False