异常值？
首先，什么是异常值，某个特别的值真的是异常值吗？
在统计学中，异常值（outlier）是与其他观察值显著不同的数据点。它们可能来源于测量方法的变异性，可能指示出新的数据情况，或者也可能是实验误差的产物。然而，判断某个值是否为异常值并不轻松，也没有一劳永逸的方法或是唯一可信的标准，因此这里整理几种常见的单变量异常值检测方法，以供自己或他人将来参考。
至于如何处理异常值，是原样保留还是剔除它们，往往会随着分析目的等因素发生变化，但请一定三思，avoid cherry-picking。
异常值检测
three-sigma rule
首先计算数据的平均值mu，与标准差std
three-sigma 法则将任何 $abs(x - µ) > 3 * std$ 的值作为异常值，这里的3也可以根据情况进行调整
这种方式与另一种Z-score方法实际是一致的：对数据进行Z-score变换，就是在计算数据的Z统计量
$$ z = (x-µ)/σ $$
对应于three-sigma方法中，就是这里的 $abs(z) > 3$会被视为异常值。
这里划定异常值的系数取3，实际已经是一个非常极端的值了，我们假设Z-score变换后的数据符合标准正态分布，此时置信区间内的概率$p(|z| > 3) = 0.9973$，即数据落在这个区间外的概率只有0.0027。
我们假设随机来一组数据：
set.seed(2501)
# 生成一组10个符合正态分布的随机数，µ = 100，σ = 10
random_numbers_1 <- round(rnorm(10, mean = 100, sd = 10))
# 99 104 103 91 106 100 96 90 110 99
# 再添加几个值，分别为68,220
data_1 <- c(random_numbers_1, 68, 220)
# 99 104 103 91 106 100 96 90 110 99 68 220
此时样本均值和标准差分别为
mu_1 = mean(data_1)
# 107.17
sigma_1 = sd(data_1)
# 37.13
此时 three-sigma的区间为$mu \pm 3*sigma = (-4.23, 218.57)$，是个大到离谱的区间了。上面的数据中只有220会被这个方式定义为异常值。这个范围对应初始的 $N(100,10^2)$ 中的概率为 $1 - 9.7^{-26}$，即有真实值落在区间外的概率为$9.7^{-26}$，因此是一个区间范围非常大，非常宽松的方法。其原因一方面是选择 3 作为系数，因此可以视情况将其调整为其他值，比如 2.5 等。
但是，计算置信区间的方式相对比较繁琐，一种更直观的是对所有数据直接计算 z 统计量
$$ |z| = \left|\frac{X - \bar{X}}{s}\right| $$
abs(round(scale(data_1)[,1],2))
# 0.22 0.09 0.11 0.44 0.03 0.19 0.30 0.46 0.08 0.22 1.05 3.04
其中68的z变量为1.05，220的z变量为3.04，如果以3为标准就可以将其作为异常值。
另一方面，该方法本身涉及到计算平均值和标准偏差，这两个值都受异常值影响很大，尤其是标准偏差，假设上面额外增加的 220 改成 200，这种方法的区间就会成为 $(10.54, 200.45)$，这时 200 落在了边界附近，却刚好不会被视为异常值。
IQR 方法
有时也称为箱线图方法，通过计算数据的分位数判断异常值，大概方式如下
Q25 = 25th_percentile
Q75 = 75th_percentile
IQR = Q75 - Q25 // inter-quartile range
IF (x < Q25 - 1.5*IQR) OR (Q75 + 1.5*IQR < x) THEN x is a mild outlier
IF (x < Q25 - 3.0*IQR) OR (Q75 + 3.0*IQR < x) THEN x is an extreme outlier
IQR 方法的好处在于，去掉异常值前后，其他值的中位数和 IQR 基本是不变的。
同样用上面的数据，计算一下这种方式的区间：
quantile(data_1, probs = c(0.25,0.75))
# Q25 = 94.75
# Q75 = 104.50
# IQR = 9.75
此时1.5倍IQR的区间为，相对还是个比较窄的范围，后添加进去的 68 和 220 在此情况下都会被视为异常值。
$$ (Q25-1.5\cdot IQR, Q75+1.5\cdot IQR) = (80.12, 119.12) $$
3 倍 IQR 的范围则为 $(65.5, 133.75)$ 相对会宽一些，此时只有220会被作为异常值。
$$ (Q25-3\cdot IQR, Q75+3\cdot IQR) = (65.5, 133.75) $$
这种方式的一个好处就是分位数对异常值不太敏感，如果把220变为200，这个区间依然是不会变的。
1.5 倍 IQR 的区间$(80.12, 119.12)$放到原始的分布中去大致对应于概率 95% 的正态区间；3 倍 IQR 的区间对应 99.94% 的概率，真实值落入这个区间外的概率也足够低了；当然相比于 three-sigma，还是严格了一些。
因此，IQR 方法更容易检测出异常值，当然也更有可能把真实值视为异常值，因此需要谨慎选择 IQR 区间的倍数。
中位数绝对偏差 MAD
上面的方法可能会遇到一些问题，包括
- 均值和标准差的计算依赖于所有数值，因此离群值可能会对均值和标准差产生较大影响
- three-sigma、IQR等方法对于小样本不一定合适（？
还有另一种可以衡量变异程度的统计量，中位数绝对偏差（Median absolute deviation，MAD）
其计算也相当简便：
对于一组数据$X_1, X_2, X_3, …, X_n$ ，其中位数为$\tilde{X} = median(X)$
MAD是各个数据与其中位数绝对偏差的中位数：
$$ MAD = median(|X_i - \tilde{X}|) $$
MAD 相比方差更能适应数据集中的异常值：在标准偏差中，与均值的距离是平方的，因此大偏差的权重更大，因此离群值会严重影响它；在 MAD 中，由于中位数天生受异常值影响较小，少数异常值的偏差几乎不影响MAD。
上面 three-sigma 使用标准方差 σ 衡量变异程度，IQR 方法使用 IQR 衡量变异，类似的，在这里使用 MAD 衡量数据的变异幅度。MAD 的使用方式类似Z-score变换中将表标准差应用于平均值的过程，为了将 MAD 作为与标准偏差 σ 一致的估计量，需要有
$$ \hat{\sigma} = k \cdot MAD $$
其中 k 是一个比例系数，取决于分布。对于正态分布的数据，k 被认为是
$$ k = 1/\Phi^{-1}(3/4) \approx 1.4826 $$
即标准正态分布$Z = (X-\mu)/\sigma$的分位数函数$\Phi^{-1}$（也称累积分布函数，R中计算使用 qnorm）的倒数。
用一开始的数据计算
# 99 104 103 91 106 100 96 90 110 99 68 220
mad_1 <- median(abs(data_1 - median(data_1)))
# mad_1 = 5.5
round(abs((data_1 - median(data_1)) / (mad_1*1.4826)), 2)
0.06 0.55 0.43 1.04 0.80 0.06 0.43 1.17 1.29 0.06 3.86 14.78
类似，视情况选择合适的值作为阈值，即可检出异常值。如果选择 3，对应的区间为 (75.04, 123.96)，概率范围为98.5%；如果选择 4，对应的区间为 (58.73, 140.27)，概率范围为99.89%。
k 的推导过程如下
为了使 $\pm k*MAD$ 的区间落在概率为 $p_0$ 的标准正态分布中区间内，有以下公式
$$ p_0 = P(|X - µ| \le k \cdot MAD) = P\left( \left|\frac{X - µ}{σ}\right| \le \frac{k \cdot MAD}{σ}\right) = P\left( |Z| \le \frac{k \cdot MAD}{σ} \right) $$
因此就有
$$ \Phi(k \cdot MAD/\sigma) - \Phi(-k \cdot MAD/\sigma) = p_0 $$
又因为
$$ \Phi(-k\cdot MAD/\sigma)=1-\Phi(k \cdot MAD/\sigma) $$
所以
$$ \Phi(k \cdot MAD/\sigma)= \frac{1 + p_0}{2} $$
令 $p_0 = 0.5$，即 50% 正态分布区间时，
$$ k = \frac{\sigma}{MAD} = \frac{1}{\Phi^{-1}(3/4)} \approx 1.4826 $$
这里选择 0.5 是因为，对于均值为零的对称分布，总体 MAD 是分布的第 75 个百分位数，因此标准正态分布的总体 MAD 刚好等于$\Phi^{-1}(3/4)$。