Grubbs, Dixon, and Generalized ESD: Turning Outlier Detection into a Hypothesis Test

import numpy as np
from scipy import stats
 
def grubbs_test(x, alpha=0.05):
    """
    Grubbs test — tek outlier için (two-sided).
 
    Returns:
        is_outlier_present: True/False
        outlier_index: int or None
        G: test istatistiği
        G_critical: kritik değer
    """
    x = np.asarray(x, dtype=float)
    n = len(x)
    if n < 3:
        raise ValueError("Grubbs için en az 3 gözlem gerek")
 
    mu = np.mean(x)
    s = np.std(x, ddof=1)  # sample std
    abs_dev = np.abs(x - mu)
    max_idx = np.argmax(abs_dev)
    G = abs_dev[max_idx] / s
 
    # Kritik değer (two-sided)
    t_val = stats.t.ppf(1 - alpha / (2 * n), n - 2)
    G_crit = ((n - 1) / np.sqrt(n)) * np.sqrt(
        t_val**2 / (n - 2 + t_val**2)
    )
 
    is_outlier = G > G_crit
    return {
        'is_outlier_present': is_outlier,
        'outlier_index': int(max_idx) if is_outlier else None,
        'outlier_value': float(x[max_idx]) if is_outlier else None,
        'G': float(G),
        'G_critical': float(G_crit),
    }
 
# Örnek
np.random.seed(42)
data = np.random.normal(loc=50, scale=5, size=30)
data = np.append(data, 100)  # bir outlier ekle
 
result = grubbs_test(data, alpha=0.05)
print(f"Outlier var mı: {result['is_outlier_present']}")
print(f"G = {result['G']:.3f}, G_crit = {result['G_critical']:.3f}")
print(f"Outlier değer: {result['outlier_value']}")

def dixon_q_test(x, alpha=0.05):
    """
    Dixon Q-test (Q10) — small sample outlier test.
    Sadece n=3-7 için lookup table dahili; n=8-30 için kabaca extend ettim.
    """
    Q_TABLE = {
        3: 0.941, 4: 0.765, 5: 0.642, 6: 0.560,
        7: 0.507, 8: 0.468, 9: 0.437, 10: 0.412,
        15: 0.338, 20: 0.300, 25: 0.277, 30: 0.260,
    }
    x_sorted = np.sort(x)
    n = len(x_sorted)
    if n < 3 or n > 30:
        raise ValueError("Dixon Q-test için 3 <= n <= 30 önerilir")
 
    q_lower = (x_sorted[1] - x_sorted[0]) / (x_sorted[-1] - x_sorted[0])
    q_upper = (x_sorted[-1] - x_sorted[-2]) / (x_sorted[-1] - x_sorted[0])
 
    # En yakın tabloya hangisi n'e karşılık?
    closest_n = min(Q_TABLE.keys(), key=lambda k: abs(k - n))
    q_crit = Q_TABLE[closest_n]
 
    return {
        'q_lower': float(q_lower),
        'q_upper': float(q_upper),
        'q_critical': float(q_crit),
        'lower_outlier': q_lower > q_crit,
        'upper_outlier': q_upper > q_crit,
        'lower_value': float(x_sorted[0]) if q_lower > q_crit else None,
        'upper_value': float(x_sorted[-1]) if q_upper > q_crit else None,
    }
 
# Kimya lab örneği — 7 ölçüm
ölçümler = [9.95, 10.01, 10.03, 10.05, 10.07, 10.10, 11.50]
result = dixon_q_test(ölçümler)
print(f"Q (üst): {result['q_upper']:.3f}, Q_crit: {result['q_critical']:.3f}")
print(f"Üst outlier: {result['upper_outlier']}")
print(f"Outlier değer: {result['upper_value']}")

def generalized_esd(x, max_outliers=10, alpha=0.05):
    """
    Generalized ESD — multiple outliers in normal data.
    """
    x = np.asarray(x, dtype=float).copy()
    n_init = len(x)
    outlier_indices = []  # orijinal indislerle
    original_indices = np.arange(n_init)
 
    R_values = []
    lambda_values = []
 
    for k in range(1, max_outliers + 1):
        if len(x) < 3:
            break
        mu = np.mean(x)
        s = np.std(x, ddof=1)
        abs_dev = np.abs(x - mu)
        max_idx_local = np.argmax(abs_dev)
        R_k = abs_dev[max_idx_local] / s
 
        n = len(x)
        p = 1 - alpha / (2 * (n - k + 1))
        t_val = stats.t.ppf(p, n - k - 1)
        lambda_k = ((n - k) * t_val) / np.sqrt(
            (n - k - 1 + t_val**2) * (n - k + 1)
        )
 
        R_values.append(R_k)
        lambda_values.append(lambda_k)
 
        if R_k > lambda_k:
            # Bu gözlem outlier — orijinal indis ekle
            outlier_indices.append(int(original_indices[max_idx_local]))
 
        # Bu gözlemi çıkar
        mask = np.arange(len(x)) != max_idx_local
        x = x[mask]
        original_indices = original_indices[mask]
 
    return {
        'outlier_indices': outlier_indices,
        'R_values': R_values,
        'lambda_values': lambda_values,
        'n_outliers': len(outlier_indices),
    }
 
# Örnek
np.random.seed(42)
data = np.concatenate([
    np.random.normal(50, 5, 100),
    [80, 90, 100, 15, 10],  # 5 outlier
])
 
result = generalized_esd(data, max_outliers=10, alpha=0.05)
print(f"Toplam outlier: {result['n_outliers']}")
print(f"Outlier indisleri: {result['outlier_indices']}")
print(f"Outlier değerler: {data[result['outlier_indices']]}")