Outlier Detection: Modified Thompson Test, Second Version¶
example of outlier test using modified Thompson technique¶
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import t
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 = np.array([28, 31, 27, 28, 29, 25, 29, 28, 18, 27])
plt.plot(1+np.arange(len(x)),x,'o',markersize=10),
plt.xlim([0,11])
plt.xlabel('$n$')
plt.ylabel('$T\; [^\circ C \,]$');
is there an outlier?¶
Sort, look at the lowest and largest and plot to visualize
[18 25 27 27 28 28 28 29 29 31]
plt.plot(np.arange(2, 10), x_1[1:-1], 'o', markersize=10)
plt.plot(1, x_1[0], 'rs', markersize=10, linewidth=2)
plt.plot(10, x_1[-1], 'rs', markersize=10)
plt.xlim([0, 11])
plt.xlabel('$\\hat{n}$')
plt.ylabel('Sorted $T\\; [^\\circ C \\,]$')
The sample mean and standard deviation, \(\bar{x}\), \(S_x\)¶
mean = 27.00, std = 3.53
Defivations for each suspicious point, take the largest \(\delta_i = |x_i - \bar{x}|\)¶
delta = np.abs(x_1 - meanx)
# print ('suspicious points are with \\delta %.1f %.1f' % (delta[0],delta[-1]))
print('suspicious point is:', x_1[np.argmax(delta)], 'deviation is = %.1f' % np.max(delta))
suspicious point is: 18 deviation is = 9.0
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.
_n = len(x_1)
_tv = t.isf(0.05 / 2, _n - 2)
_tau = _tv * (_n - 1) / (np.sqrt(_n) * np.sqrt(_n - 2 + _tv ** 2))
print('n = %d, t = %.2f, tau = %.2f' % (_n, _tv, _tau))
print('compare: %.2f to %.2f ' % (np.max(delta), _tau * stdx))
print('Is max() above $t_{\\nu,95}S$? %s ' % (np.max(delta) > _tau * stdx))
n = 10, t = 2.31, tau = 1.80
compare: 9.00 to 6.34
Is max() above $t_{\nu,95}S$? True
So we remove the outlier and repeat the test (from the beginning)¶
x_2 = np.array([28, 31, 27, 28, 29, 25, 29, 28, 27])
x_2 = np.sort(x_2)
meanx_1 = np.mean(x_2)
stdx_1 = np.std(x_2, ddof=1)
delta_1 = np.abs(x_2 - meanx_1)
print('suspicious point is:', x_2[np.argmax(delta_1)], 'deviation is = %.1f' % np.max(delta_1))
_n = len(x_2)
_tv = t.isf(0.05 / 2, _n - 2)
_tau = _tv * (_n - 1) / (np.sqrt(_n) * np.sqrt(_n - 2 + _tv ** 2))
print('n = %d, t = %.2f, tau = %.2f' % (_n, _tv, _tau))
print('compare: %.2f to %.2f ' % (np.max(delta_1), _tau * stdx_1))
print('Is max() above $t_{\\nu,95}S$? %s ' % (np.max(delta_1) > _tau * stdx_1))
suspicious point is: 25 deviation is = 3.0
n = 9, t = 2.36, tau = 1.78
compare: 3.00 to 2.95
Is max() above $t_{\nu,95}S$? True
x_3 = np.array([28, 31, 27, 28, 29, 29, 28, 27])
x_3 = np.sort(x_3)
(print('Sorted x:'), print(x_3))
meanx_2 = np.mean(x_3)
stdx_2 = np.std(x_3, ddof=1)
delta_2 = np.abs(x_3 - meanx_2)
print('suspicious point is:', x_3[np.argmax(delta_2)], 'deviation is = %.1f' % np.max(delta_2))
_n = len(x_3)
_tv = t.isf(0.05 / 2, _n - 2)
_tau = _tv * (_n - 1) / (np.sqrt(_n) * np.sqrt(_n - 2 + _tv ** 2))
print('n = %d, t = %.2f, tau = %.2f' % (_n, _tv, _tau))
print('compare: %.2f to %.2f ' % (np.max(delta_2), _tau * stdx_2))
print('Is max() above $t_{\\nu,95}S$? %s ' % (np.max(delta_2) > _tau * stdx_2))
Sorted x:
[27 27 28 28 28 29 29 31]
suspicious point is: 31 deviation is = 2.6
n = 8, t = 2.45, tau = 1.75
compare: 2.62 to 2.28
Is max() above $t_{\nu,95}S$? True
x_4 = np.array([28, 27, 28, 29, 29, 28, 27])
x_4 = np.sort(x_4)
(print('Sorted x:'), print(x_4))
meanx_3 = np.mean(x_4)
stdx_3 = np.std(x_4, ddof=1)
delta_3 = np.abs(x_4 - meanx_3)
print('suspicious point is:', x_4[np.argmax(delta_3)], 'deviation is = %.1f' % np.max(delta_3))
_n = len(x_4)
_tv = t.isf(0.05 / 2, _n - 2)
_tau = _tv * (_n - 1) / (np.sqrt(_n) * np.sqrt(_n - 2 + _tv ** 2))
print('n = %d, t = %.2f, tau = %.2f' % (_n, _tv, _tau))
print('compare: %.2f to %.2f ' % (np.max(delta_3), _tau * stdx_3))
print('Is max() above $t_{\\nu,95}S$? %s ' % (np.max(delta_3) > _tau * stdx_3))
Sorted x:
[27 27 28 28 28 29 29]
suspicious point is: 27 deviation is = 1.0
n = 7, t = 2.57, tau = 1.71
compare: 1.00 to 1.40
Is max() above $t_{\nu,95}S$? False