上週我們提到了使用 DBSCAN 找離群值，可以解決特徵超過 3 個因而無法視覺化的問題。可是 DBSCAN 真的這麼好用嗎？其實 DBSCAN 暗藏 2 個問題：
1、使用 DBSCAN 需要設定超參數eps：2 筆資料的最大距離，如果是要處理一個較陌生的資料集，這時候可能就不知道 eps 要設定多少，只能一直 trial and error。
2、如果是超高維度的資料，DBSCAN 可能會遇到維度的詛咒（curse of dimensionality），詳細說明請看這裡。
這時候該怎麼辦呢？有一篇論文（Wilkinson, 2017）提供了解決方案：HDOutliers 演算法，我們就來試試看吧。
步驟 1、將資料做 Normalise
Python 已經有提供 Normalise 的函式庫，所以我們可以一行程式就解決了。
from sklearn import preprocessing
def func_normalise(x):
return preprocessing.normalize(x, axis = 0)x = func_normalise(x)
步驟 2、計算 delta
這個步驟也很單純，因為論文已經定義好計算 delta 的公式了，我們只要先找出資料筆數跟特徵個數，就能算出 delta。
n = len(x)
p = len(x[0])
delta = 0.1 / (np.log(n)**(1.0/p))步驟 3、初始化子集代表（exemplars）以及子集資料索引（members）
HDOutliers 演算法會將所有資料分成好幾個子集，每 1 個子集都有 1 個代表的資料，我們將每 1 個子集的代表資料存在 exemplars 這個 list 裡面。
另外，我們也需要知道每 1 個子集有哪幾筆資料，這時候需要把每 1 個子集的資料索引都存在 members 裡面。因此，members 是 2 維的 list：第 1 個維度是子集 list，第 2 個維度是每一個子集所含的資料索引 list。
初始化的時候，exemplars 裡面只有首筆資料，此外 members 裡面只有1個 list，此 list 只含索引 0 。
步驟 4、將所有資料分配在不同的子集
我們取出每 1 筆資料，找到此資料跟 exemplars 裡面所有資料的最短距離（Minimum Euclidean Distance）。如果最短距離小於 delta ，我們就將此資料加入此子集；反之，我們就將這筆資料獨自建立 1 個新的子集。
請注意，這個步驟是整個演算法執行最久的部分，因為需要檢視每筆資料，這也是處理大數據的過程中最麻煩的事情。
def func_min_dist(x, exemplars):
diff = np.array(exemplars) - np.array(x)
diff = diff ** 2
dist = np.sqrt(np.sum(diff, axis = 1))
return np.min(dist), np.argmin(dist)for index in range(1, n):
min_dist, min_member = func_min_dist(x[index], exemplars)
if(min_dist < delta):
members[min_member].append(index)
else:
exemplars.append(x[index])
members.append([index])
步驟 5、計算所有子集代表之間的距離
當我們將資料分成子集之後，我們就可以計算子集代表之間的距離。看起來我們需要取出 exemplars 裡頭每筆資料，計算這筆資料跟 exemplars 裡面其他資料的距離，有點麻煩。但是，其實 Python 已經有函式可以很快算出答案：distance_matrix()。我們來試試看吧。
from scipy.spatial import distance_matrix
dist = distance_matrix(np.array(exemplars), np.array(exemplars))如果 exemplars 裡面有 k 筆資料，則計算出來的結果會得到 k by k 的矩陣，矩陣第 1 行第 2 列代表第 1 個子集代表到第 2 個子集代表的距離，以此類推。
步驟 6、使用 Exponential Distribution 擬合子集代表距離的尾巴
首先，我們把子集代表之間的距離做排序。接著取出距離較大的尾巴幾個距離資料。最後找一組參數，可以讓 Exponential Distribution貼合資料。
import copy
from scipy.stats import expondist_upper = copy.deepcopy(dist[0])
dist_upper = sorted(dist[0])dist_upper = dist_upper[len(dist_upper)//2:]loc, scale = expon.fit(dist_upper)
步驟 7、找出 Exponential Distribution累積分布到 t% 的點
舉例來說，如果我們認為 99% 的 Exponential Distribution 都屬於正常的數值，只有最後 1 % 才可能是屬於離群值。這時候我們就找分佈中，哪一個點以下的面積，累積機率剛好是 99 %。這個點稱為 cutpoint。
以我們的範例來說，超過 99.95% 才考慮是離群值，因此可以得知 cutpoint是 0.8561392559906366。
cutpoint = expon.ppf(0.9995, loc = loc, scale = scale)步驟 8、找出離群值
如果有一個子集代表到所有其他子集代表的距離大於 cutpoint，則該子集的所有資料都是屬於離群值。
dist_cutpoint = (dist > cutpoint)
dist_cutpoint_n = np.count_nonzero(dist_cutpoint, axis = 1)
outlier_index = [index for index, element in enumerate(dist_cutpoint_n) if (element == (len(exemplars) - 1))]
for sample in outlier_index:
print(members[sample])這樣我們就完成的 HDOutliers 了。這個演算法可以比較有效處理高維度的資料，而且只需要掃過全部資料1次，執行速度上也是很有優勢。
HDOutliers 演算法中，我們要設定的參數有 2 個：
1、要用 Exponential Distribution 擬合多少筆距離資料。
2、Cutpoint 的百分比應該設定多少。
讀者可能會疑惑：這樣 HDOutliers 跟 DBSCAN 都需要 trial and error 了呀！不對喔！請回頭看一下演算法，我們可能需要嘗試幾次，但是，我們依然只需要掃過全部資料 1 次就好，因為嘗試不同的參數時，我們只需要使用 exemplars 就好了，通常 exemplars 裡面的資料筆數會遠遠小於原始資料集！
文章最後，我們附上 1 個範例。此程式是要找出 MNIST 手寫數字辨識資料集，有 2 筆參有雜訊的圖片。
參考資料
Wilkinson L. (2018). Visualizing Big Data Outliers Through Distributed Aggregation. IEEE Transactions on Visualization and Computer Graphics, 24(1), pp. 256–266.
關於作者
Chia-Hao Li received the M.S. degree in computer science from Durham University, United Kingdom. He engages in computer algorithm, machine learning, and hardware/software codesign. He was former senior engineer in Mediatek, Taiwan. His currently research topic is the application of machine learning techniques for fault detection in the high-performance computing systems.
範例程式
from sklearn import preprocessing
from scipy.spatial import distance_matrix
from scipy.stats import expon
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import sys
import copydef func_normalise(x):
return preprocessing.normalize(x, axis = 0)
def func_dimension(x):
n = len(x)
p = len(x[0])
return n, p
def func_delta(n, p):
return 0.1 / (np.log(n)**(1.0/p))
def func_min_dist(x, exemplars):
diff = np.array(exemplars) - np.array(x)
diff = diff ** 2
dist = np.sqrt(np.sum(diff, axis = 1))
return np.min(dist), np.argmin(dist)
def func_exp_distribution(x, lamb):
return lamb * np.exp(-lamb * x)train = pd.read_csv('/kaggle/input/digit-recognizer/train.csv')
train_x = train.drop(['label'], axis = 1)
x = train_x.to_numpy()
for i in range(784):
if(np.random.choice([1, 0], p = [0.9, 0.1]) == 1):
x[2][i] = 255
if(np.random.choice([1, 0], p = [0.9, 0.1]) == 1):
x[3][i] = 255x = func_normalise(x)
n, p = func_dimension(x)
delta = func_delta(n, p)exemplars = [x[0]]
members = [[0]]for index in range(1, n):
min_dist, min_member = func_min_dist(x[index], exemplars)
if(min_dist < delta):
members[min_member].append(index)
else:
exemplars.append(x[index])
members.append([index])dist = distance_matrix(np.array(exemplars), np.array(exemplars))dist_upper = copy.deepcopy(dist[0])
dist_upper = sorted(dist[0])
dist_upper = dist_upper[len(dist_upper)//2:]loc, scale = expon.fit(dist_upper)cutpoint = expon.ppf(0.9995, loc = loc, scale = scale)dist_cutpoint = (dist > cutpoint)
dist_cutpoint_n = np.count_nonzero(dist_cutpoint, axis = 1)
outlier_index = [index for index, element in enumerate(dist_cutpoint_n) if (element == (len(exemplars) - 1))]
for sample in outlier_index:
print(members[sample])z = np.linspace(expon.ppf(0.0001, loc = loc, scale = scale),
expon.ppf(0.9999, loc = loc, scale = scale), 1000)plt.figure(figsize=(10, 10))
plt.plot(z, expon.pdf(z, loc = loc, scale = scale))
plt.plot([cutpoint, cutpoint],
[0, expon.pdf(z[0], loc = loc, scale = scale)],
color = 'red')
plt.show()
plt.close()
