Density-based Novelty Detection
안녕하세요. 고려대학교 LTIS 연구실 석사과정 음수민입니다.
이번 포스팅에서는 Density-based Novelty Detection에 대해 살펴보도록 하겠습니다.
해당 포스트는 고려대 강필성 교수님의 강의자료를 바탕으로 정리되었습니다.
Novelty Detection
우리는 새로운 관측치가 기존 관측치와 동일한 분포에 속하는지 혹은 다른 분포에 속하는지 여부를 결정할 수 있어야 합니다. 이러한 기능은 주로 두가지로 구분할 수 있습니다. 첫째, Outlier Detection 입니다. 이는 일반적으로 가장 집중되어있는 데이터와 멀리 떨어져 있는 관측치를 특이치로 정의하는 것입니다. 둘째, Novelty Detection 입니다. 이는 새로운 관측지를 학습데이터로 사용하지 않으며 Outlier Detection이전에 제거될 수도 있습니다.
Novel data란?
“Observations that deviate so much from other observations as to arouse suspicions that they were generated by a different mechanism(Hawkins, 1980)” “Instances that their true probability density is very low(Harmelinget al., 2006)”
이상치에 대한 정확환 정의는 위와 같습니다. 쉽게 말해 다른 관측치보다 많이 벗어나 있는 관측치가 이상치인 것입니다. 여기서 주의해야할 점은 novel data와 noise data는 다르다는 것입니다. Noise는 랜덤 에러로서 데이터 처리과정에서 제거해야하는 부분이라면 novel data는 우리가 찾고자 하는 관측치인 것입니다.
Density-based Novelty Detection
밀도에 기반한 Novelty Detection(이하 이상치 탐지)의 목적은 데이터의 분포를 추정하고 훈련된 밀도 함수에 따라 새로운 관측치가 낮은 확률을 취한다면 특이값으로 정의를 내리는 것입니다.
밀도기반 이상치 탐지법에는 일반적으로 가우시간 밀도 추정법(Gaussian density estimation), 혼합 가우시간 밀도 추정법(Mixture of Gassian estimation), 커널 밀도 추정법(Kernel density estimation),파젠 윈도우 밀도 추정법(Parzen window density estimation), 로컬 아웃라이어 팩터(Local Outlier factors)가 있습니다. 가장 먼저 가우시안 밀도 추정법에 대해 설명 드리겠습니다.
1. Gaussian Density Estimation
가우시안 밀도 추정법은 데이터가 하나의 정규분포를 따른다고 가정합니다.
위 식에서 X+는 정상영역을 뜻합니다.
가우시간 밀도 추정법에는 먼저 데이터를 스케일링하는데 민감하지 않다는 장점이 있습니다. 데이터의 각 변수별로 분산을 고려하기 때문에 따로 스케일링 하지 않아도 분산까지 고려하게 됩니다. 또한 식 자체가 주어져 있기 때문에 추가적인 분석이 가능하다는 장점이 있습니다. 함수를 미분하는 등 최적의 임계값을 계산하는 것이 가능합니다.
Maximum Likelihood Estimation (최대 우도 추정법)
최대 우도 추정법으로 이를 증명할수 있습니다.
위 그림처럼 두가지의 정규분포가 있다고 가정할 경우 x로 표시된 곳이 데이터의 위치입니다. 해당 데이터들이 어떤 분포를 따른다고 보는 것이 좋을지 판별하는 방법이 바로 최대 우도 추정법입니다.
최대 우도 추정법은 어떤 확률변수에서 표집한 값들을 토대로 그 확률변수의 모수를 구하는 방법입니다. 즉, 어떤 모수가 주어졌을 때 원하는 값들이 나올 가능성을 최대로 만드는 모수를 선택하는 방법으로 점추정 방식에 속합니다. 방정식은 아래와 같습니다.
Gamma를 1/분산로 치환하여 표현하면 아래와 같습니다.
Code Implementation
다음은 파이썬으로 작성된 코드 예시를 참조하겠습니다. 데이터는 아이리스 데이터로 사용되었습니다.
import pandas as pd
df = pd.read_csv("C:/Users/Sooomin Eum/Desktop/iris_virginica.csv")
#단별량 데이터 구현
y = pd.DataFrame.as_matrix(df[['sepallength']])
sns.distplot(y, bins=20, kde = False)
#package seabon을 이용하여 gaussian분포 추정
sns.distplot(y, fit=stats.norm, bins=20, kde=False,)
class Gaussian:
def __init__(self, mu, sigma):
self.mu = mu
self.sigma = sigma
def pdf(self, data): # 가우시안 분포 pdf 값 return
u = (data - self.mu) / abs(self.sigma)
y = (1 / (sqrt(2 * pi) * abs(self.sigma))) * exp(-u * u / 2)
return y
best= Gaussian(np.mean(y), np.std(y))
print("best mu=" , best.mu)
print("best sigma=", best.sigma)
# 가우시안 분포 그리기
x = np.linspace(3,9,200)
g_single = stats.norm(best.mu, best.sigma).pdf(x)
sns.distplot(y, bins=20, kde = False, norm_hist= True)
plt.plot(x,g_single, label = 'Single Gaussian')
plt.legend()
print(y[0:5])
#정상 boundary 외 outlier filtering
n = 0
b=0
for i in range(0,y.shape[0]):
if (stats.norm(best.mu, best.sigma).pdf(y[i])) >0.05 and (stats.norm(best.mu, best.sigma).pdf(y[i])) < 0.995:
print(y[i],"= normal")
n=n+1
else:
print(y[i],"=abnormal")
b=b+1
print("normal=",n)
print("abnormal=",b)
Result
아이리스 데이터를 토대로 가우시간 분포를 그린 결과입니다.
가우시간 분포에 의해 정상 바운더리(boundary)와 이상치를 탐지한 결과 정상 데이터는 총 145개, 비정상 데이터는 총 5개로 탐지 된것을 알 수 있습니다.
2. Mixture of Gaussian Density Estimation(MoG)
혼합 가우시안 밀도 추정법은 멀티모달 분포를 따르며, 정규분포의 선형 결합으로 이루어져 있습니다. 이는 단일 가우시안 분포보다는 작은 분산을 가지지만 학습을 위해서는 보다 많은 데이터를 요구한다는 특성을 지니고 있습니다.
다음은 혼합 가우시간 분포의 예시입니다.
위 그림은 하나의 데이터가 5개의 혼합 가우시안 분포를 따르는 것을 보여주는 예시입니다.
Component of MoG
- 정상 데이터를 구하는 확률
- 혼합 가우시안 모델의 분포
Expectation-Maximization Algorithm(E-M Algorithm)
혼합 가우시안 모델은 단일 가우시안 모델과는 달리 convex하지 않기 때문에 MLE로 쉽게 최적값을 찾아낼 수 없습니다. 따라서 휴리스틱 기법들을 적용해야 합니다. 대표적인 알고리즘인 E-M Algorithm은 Expectation을 나타내는 E단계와 Maximization을 나타내는 M단계가 존재합니다.
E단계는 로그 우도의 기댓값을 계산하는 단계이며, M단계는 기댓값을 최대화하는 단계로 이를 반복수행합니다.
다음 그림은 데이터의 분포를 추정하는 과정을 나타내는 그림입니다.
- Expectation
- Maximization
Covariance Matrix Shape
혼한 가우시안 밀도 추정법에서도 단일 가우시안 밀도 추정법과 같이 공분산 행렬이 3가지의 형태에 따라 달라지게 됩니다.
Code Implementation
다음은 파이썬으로 작성된 코드 예시입니다. 데이터는 아이리스 데이터로 사용되었습니다.
데이터의 1열과 2열을 토대로 2차원 그래프 생성
# 1열과 2열 데이터만 가지고 2차원 그래프 생성
plotsize = 8
sizeMean =10
text_size = 16
axis_font = {'fontname':'Arial', 'size':'24'}
Title_font= {'fontname':'Arial', 'size':'28'}
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.plot(y[:,0], y[:,1], 'k.', markersize= sizeMean)
for label in (ax.get_xticklabels() + ax.get_yticklabels()):
label.set_fontname('Arial')
label.set_fontsize(text_size)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.xaxis.set_ticks_position('bottom')
ax.yaxis.set_ticks_position('left')
ax.set_ylabel('Property 2', **axis_font)
ax.set_xlabel('Property 1', **axis_font)
ax.set_title('Classes based on the first two properties', y=1.08, **Title_font)
ax.figure.set_size_inches(plotsize, plotsize)
plt.show()
시각화 결과
#parameters
NProperties = y.shape[1] #variable의 갯수
NClasses = 6 # 몇개의 class로 나눌것인가 , 가우시간 분포 갯수, Hyper parameter
NObjects = y.shape[0] # 데이터의 행 데이터 갯수
# 초기 Mu, Cov 값 크기 설정
initMu = np.empty([NClasses, NProperties])
initCov = np.empty([NProperties, NProperties, NClasses])
print("initMu size =" ,initMu.shape)
print("initCov size =", initCov.shape)
#sd 를 가지는 n*n 의 양의 matrix를 만들어 내는 함수 (초기값을 주기 위해 생성)
def PosSymDefMatrix(n,sd):
M = np.matrix(np.random.rand(n,n))
M = 0.5*(M + M.T)
M = M + sd*np.eye(n)
return M
Cov = [PosSymDefMatrix(NProperties,i) for i in range(0,NClasses)]
for j in range(NClasses):
initMu[j, :] = np.random.random(NProperties) * np.amax(y, axis=0)
initCov[:, :, j] = np.mean(np.array(Cov), axis=0) + 2
print("initMu=" ,initMu) # average
print("initCov=" ,initCov) # sigma
print("Cov=", Cov)
theta = np.repeat(1.0/NClasses,NClasses)
#EM Algorithm 구현 Expectation
def EStep(y, w, mu, cov):
r_ij = np.zeros((y.shape[0], mu.shape[0]))
for Object in range(y.shape[0]):
r_ij_Sumj = np.zeros(mu.shape[0])
for jClass in range(mu.shape[0]):
r_ij_Sumj[jClass] = w[jClass] * sc.stats.multivariate_normal.pdf(y[Object, :], mu[jClass, :],cov[:, :, jClass])
for jClass in range(r_ij_Sumj.shape[0]):
r_ij[Object, jClass] = r_ij_Sumj[jClass] / np.sum(r_ij_Sumj)
return r_ij
r_n = EStep(y, initW, initMu, initCov)
print (r_n[0:5,:])
#EM Algorithm 구현 Maximazation
def MStep(r, y, mu, cov):
N = y.shape[0]
# the weigths , Estep의 평균
w_j = np.sum(r, axis=0) / N
Allmu_j = np.zeros((N, mu.shape[0], mu.shape[1]))
Allcov_j = np.zeros((N, cov.shape[0], cov.shape[1], cov.shape[2]))
# mean
for Object in range(N):
Allmu_j[Object, :, :] = np.outer(r[Object, :], y[Object, :])
mu_j = np.zeros((mu.shape[0], mu.shape[1]))
for j in range(cov.shape[2]):
mu_j[j, :] = (1 / np.sum(r, axis=0)[j]) * np.sum(Allmu_j, axis=0)[j, :]
# sd
for Object in range(N):
for j in range(cov.shape[2]):
Allcov_j[Object, :, :, j] = r[Object, j] * np.outer((y[Object, :] - mu_j[j, :]),
(y[Object, :] - mu_j[j, :]))
cov_j = np.zeros((cov.shape[0], cov.shape[1], cov.shape[2]))
for j in range(cov.shape[2]):
cov_j[:, :, j] = (1 / np.sum(r, axis=0)[j]) * np.sum(Allcov_j, axis=0)[:, :, j]
return w_j, mu_j, cov_j
w_n,mu_n,cov_n = MStep(r_n, y, initMu, initCov)
print(w_n)
print(mu_n)
print(cov_n)
# implement EM algorithm
Inititeration = 100
EMiteration = 40
lookLH = 20
for init in range(Inititeration):
# starting values
initMu = np.empty([NClasses, NProperties])
for j in range(NClasses):
initMu[j, :] = np.random.random(NProperties) * np.amax(y, axis=0)
r_n = EStep(y, initW, initMu, initCov)
w_n, mu_n, cov_n = MStep(r_n, y, initMu, initCov)
if init == 0:
logLH = -1000000000000
for i in range(EMiteration):
# E step
r_n = EStep(y, w_n, mu_n, cov_n)
# M step
w_n, mu_n, sigma_n = MStep(r_n, y, mu_n, cov_n)
# log likelihood를 계산
logLall = np.zeros((y.shape[0]))
for Object in range(y.shape[0]):
LH = np.zeros(NClasses)
for jClass in range(NClasses):
LH[jClass] = w_n[jClass] * sc.stats.multivariate_normal.pdf(y[Object, :], mu_n[jClass, :],
cov_n[:, :, jClass])
logLall[Object] = np.log(np.sum(LH))
logL = np.sum(logLall)
if i > EMiteration - lookLH:
print(logL)
if logL > logLH:
logLH = logL
print('found larger: ', logLH)
w_p = w_n
mu_p = mu_n
sigma_p = sigma_n
r_p = r_n
# 판정
mul_pdf =np.zeros(NClasses)
tot=0
for i in range(0, y.shape[0]):
n = 0
b = 0
for jClass in range(0,NClasses):
mul_pdf[jClass] = sc.stats.multivariate_normal.pdf(y[i, :], mu_p[jClass, :],sigma_p[:, :, jClass])
if mul_pdf[jClass] > 0.025 and mul_pdf[jClass] < 0.975:
n=n+1
else:
b =b+1
if n==0:
tot = tot+1
print("abnormal =", y[i,:])
print("abnormal count=", tot)
print("normal count =" ,y.shape[0]-tot)
#gaussian plot
# plot parameters
plotsize = 8
sizeMean = 10
text_size = 16
axis_font = {'fontname': 'Arial', 'size': '24'}
Title_font = {'fontname': 'Arial', 'size': '28'}
color = ['b', 'g', 'r', 'c', 'm', 'y', 'k']
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.plot(y[:, 0], y[:, 1], 'k.', markersize=sizeMean)
#gaussian plot
# plot parameters
plotsize = 8
sizeMean = 10
text_size = 16
axis_font = {'fontname': 'Arial', 'size': '24'}
Title_font = {'fontname': 'Arial', 'size': '28'}
color = ['b', 'g', 'r', 'c', 'm', 'y', 'k']
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.plot(y[:, 0], y[:, 1], 'k.', markersize=sizeMean)
for i in range(NClasses):
# the sd with ellipses
# central point of the error ellipse
pos = [mu_p[i, 0], mu_p[i, 1]]
# for the angle we need the eigenvectors of the covariance matrix
w, ve = np.linalg.eig(sigma_p[0:2, 0:2, i])
# We pick the largest eigen value
order = w.argsort()[::-1]
w = w[order]
ve = ve[:, order]
# we compute the angle towards the eigen vector with the largest eigen value
thetaO = np.degrees(np.arctan(ve[1, 0] / ve[0, 0]))
# Compute the width and height of the ellipse based on the eigen values (ie the length of the vectors)
width, height = 2 * np.sqrt(w)
# making the ellipse
ellip = Ellipse(xy=pos, width=width, height=height, angle=thetaO)
ellip.set_alpha(0.5)
ellip.set_facecolor(color[i])
ax.add_artist(ellip)
for label in (ax.get_xticklabels() + ax.get_yticklabels()):
label.set_fontname('Arial')
label.set_fontsize(text_size)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.xaxis.set_ticks_position('bottom')
ax.yaxis.set_ticks_position('left')
ax.set_ylabel('Property 2', **axis_font)
ax.set_xlabel('Property 1', **axis_font)
ax.set_title('The inferred classes based on the first two properties', y=1.08, **Title_font)
ax.figure.set_size_inches(plotsize, plotsize)
plt.show()
Result
3. Kernal-density Estimation(KDE)
커널 밀도 추정법은 데이터가 특정 분포를 따르지 않는다는 가정 하에 밀도를 추정하는 방법입니다. 모수의 분포를 모를때 밀도를 추정하는 가장 단순한 방법은 히스토그램이지만, bin의 경계에서 불연속성이 나타는 점, bin의 크기 및 시작위치에 따라 히스토그램이 달라진다는 문제 등의 한계점이 존재합니다.
이러한 연속성의 문제를 해결하기 위한 방법이 바로 커널 함수 기반 방법론입니다.
커널 함수(kernel function)은 원점을 중심으로 대칭이면서 넓이가 1인 함수로 정의되며 일반적으로 가우시안, 유니폼 함수가 있습니다.
대표적인 커널 함수는 다음과 같이 나타낼 수 있습니다.
예를들어 어떠한 분포가 있을 때 region R에 해당하는 부분 내에서 벡터 x가 나타날 확률을 P라고 하면 식은 다음과 같이 나타낼 수 있습니다.
N개의 벡터 {x1,x2,…,xn}을 가정할 때, N개의 벡터 중 k개의 벡터가 region R에 들어갈 확률은 이항분포로 나타낼 수 있습니다.
즉, 이항분포를 통해 평균과 분산을 알 수 있으며 N이 무한대로 갈 때 분산은 0에 수렴한다는 것을 알 수 있습니다.
만약 region R이 굉장히 작아서 p(x)가 그 안에서 눈에 띄는 차이를 갖지 않는다면, 다음과 같은 결과를 나타냅니다.
(단, V는 region R로 둘러쌓여있습니다.)
K를 고정시키고 V를 결정한다면 이는 K-최근접 이웃 밀도 추정법(K-nearest neighbor density estimation)이 될 것이고,
V를 고정시키고 K를 결정하면 커널 밀도 추정법이 될 것입니다.
4. Parzen window density estimation
위 그림은 데이터가 3차원의 큐브 공간안에 있을때를 나타내며 이 경우 부피는 h의 세제곱이 됩니다. 어떠한 점 x가 이 큐브공간의 정 중앙에 있다고 가정할때 확률밀도를 결정하고자 합니다.
이를 다차원의 입방체로 일반화하여 갯수를 세는 식을 만들면
이렇게 나타낼 수 있습니다. k는 입방체 안의 샘플 갯수를 나타내며, x가 주어졌을 때 이를 중심으로 입방체 안에 몇개의 샘플이 있는지 확인하는 것입니다. 위와 같은 함수는 커널 함수(kernel function)의 일종이며 파젠 윈도우(Parzen window)라고 합니다.
K함수를 p(x)에 대입하면 아래와 같습니다.
Smooth kernel function
K(U) 함수는 물연속적이며 모든 x마다 같은 가중치를 갖는다는 한계점이 있습니다. 따라서 이를 완화해주는 “스무딩(Smoothing)”이라는 기법을 사용합니다.
이렇게 표현하면 가장자리의 값도 표현할 수 있습니다.
- Smooth kernel 예시
- Uniform
- Triangular
- Epanechnikov
- Quartic
- Triweight
- Tricube
- Gaussian
- Cosine
- Logistic
- Silverman kernel
Smoothing parameter (bandwidth) h
- Large h는 밀도가 크고 완만한 분포로 추정됩니다.
- Small h는 뾰족뾰족한 밀도의 분포로 추정됩니다.
이는 h를 잘못 지정할 경우 생기는 oversmooth 현상 혹은 undersmooth 현상을 초래할 수 있음을 나타냅니다.
Code Implementation
다음은 파이썬으로 작성된 코드 예시입니다.
Window function 구현
def window_function(x_vec, unit_len=1):
"""
Implementation of the window function. Returns 1 if 3x1-sample vector
lies within a origin-centered hypercube, 0 otherwise.
"""
for row in x_vec:
if np.abs(row) > (unit_len/2):
return 0
return 1
3D-하이퍼큐브 내의 샘플링 포인트 정량화
앞에서 구현한 윈도우 함수를 사용하며, 실제로 얼마나 많은 포인트가 하이퍼큐브 내부와 외부에 있는지 정량화 한다.
X_all = np.vstack((X_inside,X_outside))
assert(X_all.shape == (10,3))
k_n = 0
for row in X_all:
k_n += window_function(row.reshape(3,1))
print('Points inside the hypercube:', k_n)
print('Points outside the hybercube:', len(X_all) - k_n)
Points inside the hypercube: 3
Points outside the hypercube: 7
Parzen-window 추정법 구현
하이퍼큐브 커널에 대한 코드 실행
def parzen_window_est(x_samples, h=1, center=[0,0,0]):
'''
Implementation of the Parzen-window estimation for hypercubes.
Keyword arguments:
x_samples: A 'n x d'-dimensional numpy array, where each sample
is stored in a separate row.
h: The length of the hypercube.
center: The coordinate center of the hypercube
Returns the probability density for observing k samples inside the hypercube.
'''
dimensions = x_samples.shape[1]
assert (len(center) == dimensions),
'Number of center coordinates have to match sample dimensions'
k = 0
for x in x_samples:
is_inside = 1
for axis,center_point in zip(x, center):
if np.abs(axis-center_point) > (h/2):
is_inside = 0
k += is_inside
return (k / len(x_samples)) / (h**dimensions)
print('p(x) =', parzen_window_est(X_all, h=1))
p(x) = 0.3
5. Local Outlier Factors(LOF)
LOF는 오브젝트(object) 근처에 있는 데이터들의 밀도까지 고려하는 추정법입니다. 해당 데이터만 고려하는 다른 방법론들과는 달리 근처에 있는 다른 데이터들의 거리와 밀도까지 상대적으로 고려한다는 특성이 있습니다.
- Definition 1 : k- distance of an object p
k-distance란 A로부터 k번재 근접 이웃까지의 거리를 뜻합니다.
- Definition 2 : k- distance neiborhood of an object p
N_k(A)는 k-distance안에 들어오는 오브젝트의 집합을 나타내며 이는 k-distance)보다 작거나 같은 거리를 가집니다.
- Definition 3 : reachability distance
A와 B까지의 거리와 k-distance중 큰 값을 사용하는 것입니다.
- Example
- reachability distance의 분포 확률
- Definition 4 : local reachability density of an object p
k- distance neiborhood of an object p를 분자로 하고 A에서 다른 오브젝트들까지의 모든 reachability distance들을 분모로 하여 나누는 값을 나타냅니다. 이는 A로부터 다른 오브젝트들 까지의 reachability distance들의 평균을 뒤집은 값과 같습니다.
- Case 1 : 만약 A가 밀도높은(denser area)지역에 위치한다면 분모는 작아질 것이고 이에따라 lrd의 값은 커지게 됩니다.
- Case 2 : 만약 A가 밀도가 높지 않은(sparse area)지역에 위치한다면 분모는 거질 것이고 이에따라 lrd의 값은 작아지게 됩니다.
- Definition 5 : local outlier factor of an object p
A에 속한 B의 lrd의 평균을 k- distance neiborhood of an object A로 나눈것입니다.
다시말해 위 그림을 예시로 들면 파란색 점이 A이고 초록색 점이 B라 하였을 때, LOF(A)값이 크다는 것은 초록색 점들의 lrd가 높고 파란색점의 lrd가 낮다는 것을 뜻합니다. 이는 초록색 점들은 밀도가 높은 지역에, 파란색 점은 밀도가 낮은 지역에 위치함을 나타냅니다.
즉, A는 밀도가 낮은지역에 있을 수록, B는 밀도가 높은 지역에 있을 수록 A의 LOF는 커지게 됩니다.
LOF는 밀도가 높은 클러스터에서 조금만 떨어져 있어도 이상치로 탐지된다는 장점이 있습니다.
Code Implementation
다음은 파이썬으로 작성된 코드 예시입니다.
import matplotlib.pyplot as plt
import numpy as np
import scipy.spatial as sc
from scipy import stats
import seaborn as sns
sns.set_style("white")
import pandas as pd
import math as mt
from random import uniform
import random
#Data Generation (IRIS, Sepal length, Sepal Width)
data = [[5.1, 3.5], [4.9, 3], [4.7, 3.2], [4.6, 3.1], [5, 3.6], [5.4, 3.9], [4.6, 3.4], [5, 3.4], [4.4, 2.9], [4.9, 3.1]]
idx = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
df = pd.DataFrame(data, columns=['X', 'Y'], index=idx)
#Array, Variables
lrd = []
inpointidx = []
knnidx = []
LOF = []
lrdko = []
a = 0
#Euclidean distance matrix
distance = pd.DataFrame(sc.distance_matrix(df.values, df.values), index=df.index, columns=df.index)
#k-distance
for i in range(1,11):
dis = np.array(distance[i])
dis = np.delete(dis, [i-1])
sortdis = np.array(distance[i])
sortdis = np.delete(sortdis, [i-1])
sortdis.sort()
kdis = sortdis[2]
#N(k)p
inpointidx1 = np.where(dis<=kdis)
knnidx.append(inpointidx1)
inpointidx = inpointidx1[0]
nkp = np.size(inpointidx)
#Rechability distance
for j in inpointidx:
dis[j] = kdis
lrdp = (nkp) / sum(dis[inpointidx])
lrd.append(lrdp)
#LOF
for k in range(0, 10):
nkpo = knnidx[k]
for l in range(0, len(nkpo[0])):
nkpoe = nkpo[0]
a = a + lrd[nkpoe[l]]
lrdko.append(a)
a = 0
for t in range(0, 10):
lof = (lrdko[t]/lrd[t]) / nkp
LOF.append(lof)
print(LOF)
[1.1619346643522603, 1.011462411773047, 0.5959147191794347, 0.9196151234054012, 0.9388214281731583, 1.488115090671161, 1.3740625717777226, 0.6781290784397118, 1.8641585236851919, 1.096774694850131]
Generation (IRIS, Sepal length, Sepal Width)
data = [[5.1, 3.5], [4.9, 3], [4.7, 3.2], [4.6, 3.1], [5, 3.6], [5.4, 3.9], [4.6, 3.4], [5, 3.4], [4.4, 2.9], [4.9, 3.1]]
idx = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
df = pd.DataFrame(data, columns=['X', 'Y'], index=idx)
# 유클리디언 거리 매트릭스 생성
distance = pd.DataFrame(sc.distance_matrix(df.values, df.values), index=df.index, columns=df.index)
# k-distance
for i in range(1,11):
dis = np.array(distance[i])
dis = np.delete(dis, [i-1])
sortdis = np.array(distance[i])
sortdis = np.delete(sortdis, [i-1])
sortdis.sort()
kdis = sortdis[2]
#N(k)p
inpointidx1 = np.where(dis<=kdis)
knnidx.append(inpointidx1)
inpointidx = inpointidx1[0]
nkp = np.size(inpointidx)
#Rechability distance
for j in inpointidx:
dis[j] = kdis
lrdp = (nkp) / sum(dis[inpointidx])
lrd.append(lrdp)
kdis
0.3000000000000007
dis
array([0.4472136 , 0.3 , 0.3 , 0.3 , 0.50990195, 0.94339811, 0.42426407, 0.31622777, 0.53851648])
idx
[(array([3, 6, 8], dtype=int64),),
(array([1, 2, 8], dtype=int64),),
(array([2, 5, 8], dtype=int64),),
(array([2, 5, 7], dtype=int64),),
(array([0, 5, 6], dtype=int64),),
(array([0, 4, 6], dtype=int64),),
(array([2, 3, 6], dtype=int64),),
(array([0, 4, 8], dtype=int64),),
(array([1, 2, 3], dtype=int64),),
(array([1, 2, 3], dtype=int64),),
(array([3, 6, 8], dtype=int64),),
(array([1, 2, 8], dtype=int64),),
(array([2, 5, 8], dtype=int64),),
(array([2, 5, 7], dtype=int64),),
(array([0, 5, 6], dtype=int64),),
(array([0, 4, 6], dtype=int64),),
(array([2, 3, 6], dtype=int64),),
(array([0, 4, 8], dtype=int64),),
(array([1, 2, 3], dtype=int64),),
(array([1, 2, 3], dtype=int64),)]
nkp
3
lrd
[2.2360679774997916,
3.162277660168372,
4.472135954999576,
3.3333333333333353,
2.236067977499788,
1.56173761888606,
2.499999999999998,
3.162277660168382,
1.9611613513818404,
3.3333333333333255,
2.2360679774997916,
3.162277660168372,
4.472135954999576,
3.3333333333333353,
2.236067977499788,
1.56173761888606,
2.499999999999998,
3.162277660168382,
1.9611613513818404,
3.3333333333333255]
#LOF
for k in range(0, 10):
nkpo = knnidx[k]
for l in range(0, len(nkpo[0])):
nkpoe = nkpo[0]
a = a + lrd[nkpoe[l]]
lrdko.append(a)
a = 0
for t in range(0, 10):
lof = (lrdko[t]/lrd[t]) / nkp
LOF.append(lof)
print(LOF)
result
[1.1619346643522603,
1.011462411773047,
0.5959147191794347,
0.9196151234054012,
0.9388214281731583,
1.488115090671161,
1.3740625717777226,
0.6781290784397118,
1.8641585236851919,
1.096774694850131,
1.1619346643522603,
1.011462411773047,
0.5959147191794347,
0.9196151234054012,
0.9388214281731583,
1.488115090671161,
1.3740625717777226,
0.6781290784397118,
1.8641585236851919,
1.096774694850131]
이상 본 포스팅을 마치겠습니다. 부족하지만 긴 글 읽어주셔서 감사합니다.
References
[1] P.S. Kang, “Novelty Detection.” Business Analytics. Korea University, Fall 2018. Lecture.
[2] https://jayhey.github.io/novelty%20detection/2017/10/18/Novelty_detection_overview/
[3] https://jayhey.github.io/novelty%20detection/2017/11/02/Novelty_detection_Gaussian/
[4] https://scikit-learn.org/stable/modules/outlier_detection.html
[5] https://sebastianraschka.com/Articles/2014_kernel_density_est.html#232-quantifying-the-sample-points-inside-the-3d-hypercube