비모수 밀도 추정을 통한 클러스터링
기본 개념
라벨링되어 있지 않은 데이터 포인트들을 비슷한 특성을 가진 한 개 이상의 군집으로 묶는 것을 클러스터링이라고 부르고, 이러한 군집화에는 다양한 접근법이 있습니다. 그 중 하나의 아이디어는 데이터 공간에서 데이터 포인트가 가깝게 모여 있는, 굉장히 밀도가 높은 어떤 특정 지역이 있어서 다른 지역들과 구분될 것이고 이 지역에 모인 포인트들은 하나의 군집이라고 할 수 있다는 것입니다. 그러면 얼마나 밀도가 높으면 하나의 성립된 군집이라고 할 수 있을까요?
$d$ 차원의 확률밀도함수 $f$ 에 대해 다음과 같이 $R(c)$ 와 $p_c$ 를 정의합니다.
\[R(c) = \{ x : x \in \mathbb{R}^d, f(x) \ge c \}, \quad (0 \le c \le \max f )\] \[p_c = \int_{R(c)} f(x) dx\]여기서 $c$ 는 ‘얼마나 밀도가 높은지’의 ‘얼마나’를 뜻하는 값입니다. 즉 $R(c)$는 특정 기준 이상의 밀도를 지니는 지역 혹은 구간(region)이고, $p_c$ 는 해당 구간에 대해 확률밀도함수를 적분한 것= 즉 어떤 데이터포인트가 그 지역에 속할 확률입니다. 만약에 $f$ 가 단봉형(unimodal)이라면 $R (c)$은 연결된 구간일 것이고, 그렇지 않다면 연결되어 있을 수도 아닐 수도 있죠.
이 그림은 d 가 1일 때, 즉 1차원일 때의 예시입니다. 좌측 그림이 확률밀도함수 $f$ 인데, 양봉형(bimodal)이고, 우리가 특정 값의 $c$를 어떻게 정하느냐에 따라 1개의 혹은 2개의 $R(c)$ 를 얻습니다. 또 우측 그림처럼 $p_c$와 연결된 지역의 수의 관계를 일종의 계단 함수 $m(p_c)$로 생각해볼 수 있습니다. $c$의 값이 작아진다고 가정할 때 (확률밀도함수에서 확률이란 적분값, 그림에서 쉽게 파악하기로는 면적이므로) $c$가 낮아지면 노란색 면적은 커지고, 따라서 그 면적에 속할 확률인 $p_c$는 커지고, 연결된 지역의 수는 특정 포인트를 지나면 2개로 늘어났다가, 더 내려가면 하나의 큰 지역으로 바뀝니다(우측 그림에서 쭉 오른쪽으로 이동하는 겁니다).
이게 기본 개념들이지만, 우리는 대부분의 경우 진짜 분포 $f$를 모릅니다. 관찰된 샘플 데이터 $S = x_1, \cdots\, x_n$ 를 가지고 변수 $X$의 분포를 추정했다고 합시다. 어떤 분포를 따른다고 가정하지 않고 비모수적으로 추정된 밀도 $\hat{ f}(x)$ 가 있다고 한다면, 우리가 가지게 되는 것도 $R(c)$ 대신 $\hat {R}(c)$일 것입니다. 우리의 목적은 가지고 있는 이 샘플들을 클러스터에 할당하는 것이므로, 다음과 같이 다시 씁니다.
\[S(c) = \{ x_i : x_i \in S, \hat {f }(x_i) \ge c \}, \quad (0 \le c \le \max \hat{f})\] \[\hat {p_c} = \frac {\vert S(c) \vert }{n} \]이 문제는 데이터가 1차원일 때는 매우 간단해 보이지만, 차원이 높아질수록 해결하기 복잡해질 것입니다. 조금 더 자세히는, 각 데이터포인트 x가 아무 주어진 집합 $R(c)$ 에 속하는지 마는지를 파악하는 건 쉽지만, 그 $R(c)$가 몇 개의 연결된 집합(클러스터)로 이루어져 있는지를 알아내는 게 어렵다는 겁니다. 사실 그게 클러스터링의 전부이기도 한데요. 이 접근에서 제시된 해결 방법은 두 가지입니다.
방법1. 그래프 분할
첫번째 방법을 이해하려면 들로네 삼각 분할(Delaunay triangulation)과 보로노이 모자이크(Voronoi tessellation)의 개념을 알아야 합니다. 저는 여기의 도움을 많이 받았습니다! 제일 먼저 해야하는 것은 샘플 S의 들로네 삼각 분할을 만든 다음, 모든 $x_i \in S$ 를 가지고 비모수 밀도 추정을 해서 $\hat {f} (x_i)$ 를 얻습니다. 그리고 모든 가능한 $p_c \in (0,1)$ 값에 대해, $\hat{f} (x_i) <c$ 인 포인트들을 제거합니다(쉽게 생각하면 길이가 긴 변은 제거해서 모여 있는 애들만 남긴다고 보면 됩니다). 이 과정을 시각화한다면 아래 그림처럼 보일 겁니다.
2차원 데이터에서 예시입니다. 좌측 그림에서 점선이 voronoi tessellation이고, 실선이 superimposed Delaunay triangulation입니다. 밀도가 낮은 지점들을 모두 제외하고 나면 우측 그림에서 2개의 연결 집합들만 남은 것이 보이죠?
이 단계들을 모든 가능한 $p_c$ 값에 대해 반복하고, 그 과정에서 샘플 내 원소들이 할당되는 집단의 정보를 계속 기록해 두겠습니다. 그러면 일종의 연결집합들의 트리를 만들 수 있는데, 이 트리의 잎은 밀도 함수의 최빈값이고, $p_c$ 의 수준에 따라 그 개수가 달라지게 됩니다. 이 최빈값이 바로 클러스터 중심이며, 각 데이터포인트들은 이 중심에 가까우면 해당 클러스터로 결정됩니다.
이 그림은 2차원 데이터 기준으로 데이터포인트들을 각 클러스터 중심에 할당한 것입니다. 트리 구조로 보면 오른쪽 그림과 같습니다. y축을 밀도를 거꾸로 표시한 것으로 생각하면 됩니다. 트리의 루트는 밀도 0에 대응하고, 그룹 1로 마킹 되어있는 부분이 제일 높은 최빈값입니다. 밀도 수준을 어디서 끊을지에 따라 실제로 결정되는 클러스터의 개수가 달라집니다. 좌측 그림을 보면 아직 색깔이 칠해지지 않은 데이터포인트들이 보일 것입니다. 다음 단계는 이 모든 남은 포인트들을 정해진 M개의 클러스터 중심에 할당하는 것입니다.
방법2. 계곡이 없는 쌍 연결
두번째 방법은 특정 구간 $[ x_1, x_2 ]$ 을 따라 이동할 때 $\hat{f} (x)$ 의 변화를 생각해보자는 아이디어에서 출발합니다. 다시 1차원의 밀도함수 그림으로 돌아가서, 만약 형성된 2개 이상의 연결집합으로 이루어진 $\hat{R}(c)$ 에 대해 $x_1$ 와 $x_2$ 가 같은 연결집합에 속한다면, 그 구간 내 밀도함수는 줄어들었다가 다시 증가하는 부분, 즉 local minimum을 가지고 있지 않을 겁니다. 그러나 만약 다른 연결집합에 속하면, 그 사이에 분명 일종의 계곡 같은 것이 존재하게 될 것입니다.
이 관점으로 그래프를 만들어볼게요. 모든 샘플 포인트들을 노드로 하되, 각 포인트들과의 사이에 ‘계곡이 없는’ 다른 한 쌍의 포인트를 에지로 연결하는 그래프라고 합시다. 이때 밀도함수는 어느 구간에서든 계속 변동할 수 있으므로, 얼만큼 떨어졌다 올라와야 할 건지, 즉 그 정도가 어느 수준 이상이어야 계곡이 있는 것으로 볼 것인지를 정해야 할 겁니다. 즉 계곡의 강도를 나타내는 정도 $V$를 수치화할 필요가 있습니다.
좌측 그림에서 $x_1$ 와 $x_2$ 는 다른 클러스터에 속하고, 그럼 분명 둘 사이에는 어떤 밀도의 계곡이 존재할 것 같습니다. 우측 그림처럼요. 여기서 보조 함수 $\varPhi$ 가 필요한데요, 이 함수는 계곡에 넘치지 않을 정도로만 물을 채우는 함수입니다. 물을 채운 부분을 $A_1$ 이라 하고, 이 부분의 밀도함수 면적을 A_2라 할 때, 계곡의 강도는 다음과 같습니다.
\[V = \frac{A_1}{A_1+A_2} \in [ 0, 1) \]따라서 정해진 파라미터 $\lambda \in (0,1)$ 에 대해 $V < \lambda$ 이면 계곡이 존재하지 않는다, 즉 $x_1$ 와 $x_2$ 는 같은 연결집합에 속한다고 판단하도록 하겠습니다. 이렇게 모든 데이터포인드들의 쌍에 대해 연결 여부를 판단하고 나면, 그 다음 단계는 이전 그래프 분할 방식과 동일합니다.
클러스터링 결과 평가하기
그러면 클러스터링이 잘 되었는지 어떻게 알 수 있을까요? 일반적으로는 실루엣 점수라는 것을 사용합니다. 동일 클러스터로 할당한 데이터 포인트들 간 거리가 가깝고, 다른 클러스터에 속한 포인트들과의 거리가 멀수록 실루엣 점수가 높고 클러스터링이 잘 된 것이라고 봅니다. 밀도 기반 실루엣(dbs; density-based silhouette)은 유사한 개념이지만, 포인트들 간 거리를 확률 로그비로 바꾼 것입니다. 우선 $x_i$ 에 대해 $\pi_m$ 은 $m$번째 클러스터에 대한 사전(prior) 확률입니다. 사전에 클러스터 구성에 대한 정보가 없다면 균등분포에 따라 설정하게 되겠죠.
\[\hat{\tau_m} (x_i) = \frac{ \pi m \hat {f}m (x_i)}{ \sum ^M _{k=1} \pi_k \hat{f} _k (x_i)}\]이걸 가지고 다음과 같이 dbs를 계산합니다.
\[dbs_i = \frac {log(\frac{\hat{\tau_m0} (x_i) }{ \hat{\tau_m1} (x_i)} ) } { \max \vert log(\frac{\hat{\tau_m0} (x_i) }{ \hat{\tau_m1} (x_i)}) \vert } \]이때 $m_0$ 은 $x_i$ 가 할당된 클러스터고, $m_1$ 은 $\hat{\tau_m(x_i)}$ 가 두번째로 높은 값이 되는 클러스터입니다. 즉 어떤 한 포인트의 dbs 값은, 그 포인트가 할당된 클러스터에 속할 사후 확률과, 다른 클러스터에 속할 사후 확률의 최대값 간의 로그비와 비례합니다. 이게 클수록 해당 클러스터에 속할 확률이 다른 클러스터에 속할 확률에 비해 매우 크고, 따라서 그 해당 클러스터에 확실히 할당된 것이므로 군집이 잘 나누어졌다고 해석할 수 있습니다.
R로 해보기
R에서는 이러한 과정을 pdfCluster 패키지를 이용하면 간단하게 해볼 수 있습니다. 오늘의 예제 데이터는 와인 데이터입니다. 이 데이터는 178종의 와인에 대해 13개 변수를 가지고 있으며, 3종류로 라벨링이 되어 있습니다. 그 중 (용이한 시각화를 위해) 3개의 변수만 이용해 밀도 추정부터 먼저 해보겠습니다.
1
2
3
4
5
> library("pdfCluster")
> data("wine", package = "pdfCluster")
> winesub <- wine[,c(2,5,8)]
> pdf <- kepdf(winesub)
> plot(pdf)
다음과 같이 클러스터링 후 원래의 라벨과 비교해 봅니다.
1
2
3
4
5
6
7
> cl.winesub <- pdfCluster(winesub)
#두번째 방법을 사용하고 싶다면 pdfCluster(winesub, graphtype = "pairs")
> table(wine[, 1], groups(cl.winesub))
1 2 3
Barolo 58 1 0
Grignolino 4 62 5
Barbera 0 0 48
plot 메서드로 최빈값 함수와 pairwise 스캐터 플랏, 클러스터 트리와 마지막으로 dbs까지 그려볼 수 있습니다. 정오표뿐 아니라 dbs 플랏을 봐도 (거의 모든 관측치의 dbs 값이 0보다 큼) 데이터 포인트들이 꽤 잘 묶인 것으로 보이네요.
사족
가끔 뭔가 검색어를 잘못 입력 넣거나 이상한 데로 빠져서 원래 검색하려던 게 아닌 것을 한참 보게 되는 경우가 있는 것 같습니다. 사실 저는 이번주 업무 중에 각 유저들의 특정 행동에 대한 분포를 추정한 다음 그 분포(유저 수만큼 추정된) 클러스터링하고 싶어서, 그러니까 추정된 밀도들을 각각 하나의 데이터포인트로 하고 클러스터링하는 방법을 알고 싶어서 검색했다가 위 내용을 발견했습니다. 이건가? 하고 한참 읽어보니까 제가 원했던 내용이 아니더라고요. 그래도 재미있게 읽었기 때문에 블로그에 글을 남겨 둡니다.