[핸즈온 머신러닝] chapter 9. 비지도 학습

 

 

 

In [2]:

# 파이썬 ≥3.5 필수
import sys
assert sys.version_info >= (3, 5)

# 사이킷런 ≥0.20 필수
import sklearn
assert sklearn.__version__ >= "0.20"

# 공통 모듈 임포트
import numpy as np
import os

# 노트북 실행 결과를 동일하게 유지하기 위해
np.random.seed(42)

# 깔끔한 그래프 출력을 위해
%matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt
mpl.rc('axes', labelsize=14)
mpl.rc('xtick', labelsize=12)
mpl.rc('ytick', labelsize=12)

plt.style.use('default')

 

 

비지도 학습¶

예산을 많이 받아먹고 각광받는 것은 지도학습이지만 현실은 데이터에 레이블이 없는 경우가 허다하다.
이 경우에 모든 샘플에 레이블을 일일이 붙여야 하지만, 바로 이용할 수 있는 방법도 있다.

• 군집 (clustering)
  비슷한 샘플을 클러스터로 모아준다.
• 이상치 탐지 (outlier detection)
  정상 데이터가 어떻게 보이는지를 학습하고 이 데이터와 멀리 떨어진 경우 이상치로 판단한다.
• 밀도 추정
  데이터셋 생성 확률 과정 (random process)의 확률 밀도 함수를 추정한다.
  이 함수는 이상치 탐지에 사용되고 밀도가 낮은 곳에 위치한 데이터일수록 이상치일 확률이 올라간다

 

 

군집 (clustering)¶

레이블을 몰라도 그냥 비슷해보이는 샘플끼리 묶어 주는 것은 가능하다.
묶어주는 단위를 클러스터라고 부르고 분류와 마찬가지로 각 샘플은 하나의 클러스터에 묶인다.

k-평균¶

센트로이드라는 기준 점이 되는 데이터 포인트를 랜덤하게 선정한다. (우선 무작위로 k개의 샘플을 뽑는다)
그 다음 모든 샘플에 가까이 있는 센트로이드의 클러스터를 할당하고, 새롭게 센트로이드를 업데이트 한다.
이 과정을 반복하면 제한된 횟수 안에 수렴하는 것을 보장한다.
(계산복잡도는 샘플, 클러스터, 차원 개수에 선형적이다. 단, 군집할 수 없는 데이터 구조라면 폭증한다)
그런데 센트로이드 초기화를 할 때 잘못하면 최적의 결과에 도달하지 못할 수 있다.
센트로이드 초기화를 최적화하여 해결할 수 있는데 그냥 랜덤 초기화를 여러 번 시행하여 가장 좋은 솔루션을 얻는다.
최선의 솔루션을 얻을 수 있는 지표는 가장 가까운 센트로이드 사이의 평균제곱 거리이며 이를 이너셔(inertia)라고 부른다.
사이킷런의 score()메서드는 이너셔의 음숫값을 반환하는데, 큰것이 좋다는 사이킷런의 큐칙을 따라야하기 때문이다.

 

In [3]:

from sklearn.datasets import make_blobs

blob_centers = np.array([[0.2, 2.3], [-1.5, 2.3], [-2.8, 1.8], [-2.8, 2.8],
                         [-2.8, 1.3]])
blob_std = np.array([0.4, 0.3, 0.1, 0.1, 0.1])
X, y = make_blobs(n_samples=2000,
                  centers=blob_centers,
                  cluster_std=blob_std,
                  random_state=7)


def plot_clusters(X, y=None):
    plt.scatter(X[:, 0], X[:, 1], c=y, s=1)
    plt.xlabel("$x_1$", fontsize=14)
    plt.ylabel("$x_2$", fontsize=14, rotation=0)


plt.figure(figsize=(8, 4))
plot_clusters(X)
plt.show()

 

 

 

In [4]:

from sklearn.cluster import KMeans

k = 5
kmeans = KMeans(n_clusters=k)
y_pred = kmeans.fit_predict(X)
print(y_pred)
print(kmeans.labels_) # KMeans 클래스의 인스턴스는 labels_ 인스턴스 변수에 훈련된 샘플의 레이블을 가짐
print(kmeans.cluster_centers_) # 찾은 센트로이드 (k=5)를 저장한다.
plot_clusters(X)
for cent in kmeans.cluster_centers_:
    plt.plot(cent[0],cent[1],"ro")
plt.show()

 

 

[4 0 1 ... 2 1 0]
[4 0 1 ... 2 1 0]
[[-2.80389616  1.80117999]
 [ 0.20876306  2.25551336]
 [-2.79290307  2.79641063]
 [-1.46679593  2.28585348]
 [-2.80037642  1.30082566]]

 

 

 

k-평균++¶

k-평균++ 알고리즘이 2006년에 발표되었는데 최적이 아닌 솔루션으로 수렴할 가능성을 크게 낮춰주었다.
우선 데이터셋에서 무작위로 균등하게 하나의 센트로이드를 선택하고 현재 샘플의 가장 가까운 센트로이드까지의 거리 제곱이 모든 샘플에 대하여 샘플과 가장 가까운 센트로이드와의 거리의 제곱 중에서 차지하는 비율을 확률로 이용하여
하나의 샘플을 새로운 센트로이드로 선택한다.
이 확률 분포는 이미 선택한 센트로이드와 멀리 떨어진 샘플을 센트로이드로 선정하게 도와준다.

k- 평균 속도 개선과 미니배치¶

불필요한 거리 계산을 많이 피함으로써 알고리즘의 속도를 높일 수 있다.
이는 삼각 부등식과 샘플과 센트로이드 사이의 거리에 대한 하한과 상한선을 정해줌으로써 구현할 수 있다.
이 방식은 KMeans에서 기본값으로 이용된다.

또는 데이터를 전체를 한꺼번에 보내주는 것보다 미니배치를 이용하여 센트로이드를 조금씩 이동시킨다.
이를 이용하면 알고리즘의 속도를 3배에서 4배 정도 높일 수 있으며 메모리에 들어가지 않는 대량의 데이터셋에
클러스터링 알고리즘을 적용할 수 있다.
사이킷런에는 MiniBatchKMeans에 구현되어 있으며 KMeans와 같은 방법으로 이용이 가능하다.

 

In [5]:

from sklearn.cluster import MiniBatchKMeans

minibatch_kmeans = MiniBatchKMeans(n_clusters=5)
minibatch_kmeans.fit(X)

 

 

C:\Users\hesh0\anaconda3\lib\site-packages\sklearn\cluster\_kmeans.py:887: UserWarning: MiniBatchKMeans is known to have a memory leak on Windows with MKL, when there are less chunks than available threads. You can prevent it by setting batch_size >= 5120 or by setting the environment variable OMP_NUM_THREADS=1
  warnings.warn(

Out[5]:

MiniBatchKMeans(n_clusters=5)

 

 

최적의 클러스터 개수(k) 찾기¶

클러스터의 개수가 적으면 별개의 클러스터를 합치고 많을수록 같은 클러스터를 나눠버리게 된다.
이너셔는 k가 증가할 수록 감소하므로 손실함수처럼 이용하기에는 좋지않다.
물론 k를 늘려가면서 이너셔의 감소율을 보다가 느리게 감소하기 시작하는 지점인 엘보를 찾으면 되지만 너무 엉성하다.
더 정확한 방법은 실루엣 점수를 구하는 것이다. 이 값은 모든 샘플에 대한 실루엣 계수의 평균이다.
샘플의 실루엣 계수는 다음과 같은 식으로 구해진다.
$$\begin{equation} (b-a)\; / \max{(a,b)} \end{equation}$$ (단, a는 동일한 클러스터에 있는 다른 샘플까지 평균 거리이고 b는 가장 가까운 다른 클러스터까지 평균 거리)
계수가 1에 가까울수록 자신이 속한 클러스터에 잘 속해있으며 0에 가까울 수록 경계에 위치해있으며
-1에 가깝다면 샘플이 잘못된 클러스터에 위치해있음을 의미한다.

 

In [6]:

from sklearn.metrics import silhouette_score
silhouette_score(X,kmeans.labels_)

 

Out[6]:

0.655517642572828

 

 

모든 모든 샘플의 실루엣 계수를 클러스터와 계수값으로 정렬한 실루엣 다이어그램을 살펴보자.
높이는 클러스터가 포함한 샘플의 개수를, 너비는 샘플의 실루엣 계수를 나타낸다.
수직 파선은 클러스터 개수에 해당하는 실루엣 점수를 나타내는데, 이 점수보다 낮은 계수를 가지면
클러스터의 샘플이 다른 클러스터랑 가깝다는 것을 의미한다.
k값을 정할 때 모든 샘플이 해당 파선을 넘기고 굵기가 동일한 것을 골라야 한다.

 

In [7]:

from sklearn.metrics import silhouette_samples
from matplotlib.ticker import FixedLocator, FixedFormatter

plt.figure(figsize=(11, 9))
kmeans_per_k = [
    KMeans(n_clusters=k, random_state=42).fit(X) for k in range(1, 10)
]
silhouette_scores = [
    silhouette_score(X, model.labels_) for model in kmeans_per_k[1:]
]

for k in (3, 4, 5, 6):
    plt.subplot(2, 2, k - 2)

    y_pred = kmeans_per_k[k - 1].labels_
    silhouette_coefficients = silhouette_samples(X, y_pred)

    padding = len(X) // 30
    pos = padding
    ticks = []
    for i in range(k):
        coeffs = silhouette_coefficients[y_pred == i]
        coeffs.sort()

        color = mpl.cm.Spectral(i / k)
        plt.fill_betweenx(np.arange(pos, pos + len(coeffs)),
                          0,
                          coeffs,
                          facecolor=color,
                          edgecolor=color,
                          alpha=0.7)
        ticks.append(pos + len(coeffs) // 2)
        pos += len(coeffs) + padding

    plt.gca().yaxis.set_major_locator(FixedLocator(ticks))
    plt.gca().yaxis.set_major_formatter(FixedFormatter(range(k)))
    if k in (3, 5):
        plt.ylabel("Cluster")

    if k in (5, 6):
        plt.gca().set_xticks([-0.1, 0, 0.2, 0.4, 0.6, 0.8, 1])
        plt.xlabel("Silhouette Coefficient")
    else:
        plt.tick_params(labelbottom=False)

    plt.axvline(x=silhouette_scores[k - 2], color="red", linestyle="--")
    plt.title("$k={}$".format(k), fontsize=16)

plt.show()

 

 

C:\Users\hesh0\anaconda3\lib\site-packages\sklearn\cluster\_kmeans.py:881: UserWarning: KMeans is known to have a memory leak on Windows with MKL, when there are less chunks than available threads. You can avoid it by setting the environment variable OMP_NUM_THREADS=8.
  warnings.warn(

 

 

 

k평균의 한계¶

속도도 빠르고 확장도 용이한 k평균에도 단점은 존재한다.

1. 최적이 아닌 솔루션을 피하기 위해 알고리즘을 여러 번 실행시켜야 한다.
2. 클러스터 개수를 직접 지정해야 한다.
3. 클러스터의 크기나 밀집도가 서로 다르고 원형이 아니라면 잘 작동하지 않는다.

군집을 이용한 이미지 분할¶

말 그대로 이미지를 여러 세그먼트로 나눠버리는 작업이다.
시맥틱 분할에서는 동일한 종류의 물체에 속한 모든 픽셀을 같은 세그먼트에 할당시킨다.
좋은 성능을 내기 위해서는 CNN을 사용해야 하지만 우선 색상 분할로 진행해보자.
이 알고리즘에서는 동일한 색상을 가진 픽셀을 같은 세그먼트에 할당시킨다.

 

In [8]:

import urllib.request

PROJECT_ROOT_DIR = "."
images_path = os.path.join(PROJECT_ROOT_DIR, "images", "unsupervised_learning")
os.makedirs(images_path, exist_ok=True)
DOWNLOAD_ROOT = "https://raw.githubusercontent.com/rickiepark/handson-ml2/master/"
filename = "ladybug.png"
print("Downloading", filename)
url = DOWNLOAD_ROOT + "images/unsupervised_learning/" + filename
urllib.request.urlretrieve(url, os.path.join(images_path, filename))

from matplotlib.image import imread

image = imread(os.path.join(images_path, filename))
image.shape

 

 

Downloading ladybug.png

Out[8]:

(533, 800, 3)

 

 

각 픽셀에 대한 RGB값(0~1)이 담긴 3D 벡터에 대하여 RGB색상의 리스트로 변환한다.
그 다음 이를 KMeans를 이용하여 색상을 클러스터로 모아준다.

 

In [8]:

X = image.reshape(-1,3)
segmented_imgs = []
n_colors = (10, 8, 6, 4, 2)
for n_clusters in n_colors:
    kmeans = KMeans(n_clusters=n_clusters, random_state=42).fit(X)
    segmented_img = kmeans.cluster_centers_[kmeans.labels_]
    segmented_imgs.append(segmented_img.reshape(image.shape))
plt.figure(figsize=(10,5))
plt.subplots_adjust(wspace=0.05, hspace=0.1)

plt.subplot(231)
plt.imshow(image)
plt.title("Original image")
plt.axis('off')

for idx, n_clusters in enumerate(n_colors):
    plt.subplot(232 + idx)
    plt.imshow(segmented_imgs[idx])
    plt.title("{} colors".format(n_clusters))
    plt.axis('off')

plt.show()

 

 

 

 

비슷한 색들을 하나의 클러스터로 묶어줄 수 있지만 묶을 수 있는 색의 개수가 줄어듬에 따라 주위 색에 합쳐진다.
이는 KMeans가 비슷한 크기의 클러스터를 만들려는 경향이 있기 때문이다.

 

 

군집을 사용한 전처리¶

군집은 지도 학습을 적용하기 전에 전처리 단계에 이용될 수 있다. 예를 들어 차원 축소에 군집을 사용할 수 있다.

 

In [9]:

from sklearn.datasets import load_digits

X_digits, y_digits = load_digits(return_X_y=True)

from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(X_digits,
                                                    y_digits,
                                                    random_state=42)

 

In [10]:

from sklearn.linear_model import LogisticRegression

log_reg = LogisticRegression(multi_class="ovr",
                             solver="lbfgs",
                             max_iter=5000,
                             random_state=42)
log_reg.fit(X_train, y_train)
log_reg.score(X_test, y_test)

 

Out[10]:

0.9688888888888889

 

In [11]:

from sklearn.pipeline import Pipeline

pipeline = Pipeline([("kmeans", KMeans(n_clusters=50)),
                     ("log_reg",
                      LogisticRegression(multi_class="ovr",
                                         solver="lbfgs",
                                         max_iter=5000,
                                         random_state=42))])
pipeline.fit(X_train, y_train)
pipeline.score(X_test, y_test)

 

Out[11]:

0.9777777777777777

 

In [ ]:

from sklearn.model_selection import GridSearchCV

param_grid = dict(kmeans__n_clusters=range(2,100))
grid_clf = GridSearchCV(pipeline, param_grid, cv=3)
grid_clf.fit(X_train,y_train)
grid_clf.best_params_

 

In [ ]:

print(grid_clf.score(X_test,y_test))

 

 

군집을 사용한 준지도 학습¶

레이블이 없는 데이터가 많고 레이블이 있는 데이터가 적을 때 이용 가능한 방법인 준지도 학습.
훈련 세트를 k개의 클러스터로 모으고 각 클러스터에서 센트로이드에 가장 가까운 대표 이미지를 찾는다.
대표 이미지만으로 구성된 적은 수의 샘플들이더라도 정확도가 비약적으로 상승하게 된다.

 

In [15]:

n_labeled = 50
log_reg = LogisticRegression(max_iter = 1000)
log_reg.fit(X_train[:n_labeled],y_train[:n_labeled])
log_reg.score(X_test,y_test)

 

Out[15]:

0.8222222222222222

 

In [16]:

k = 50
kmeans = KMeans(n_clusters=k, random_state=42)
X_digits_dist = kmeans.fit_transform(X_train)

# 각 클러스터마다 센트로이드에 가장 가까운 대표 이미지를 찾는다.
representative_digit_idx = np.argmin(X_digits_dist, axis=0)
X_representative_digits = X_train[representative_digit_idx]

plt.figure(figsize=(8, 2))
for index, X_representative_digit in enumerate(X_representative_digits):
    plt.subplot(k // 10, 10, index + 1)
    plt.imshow(X_representative_digit.reshape(8, 8), cmap="binary", interpolation="bilinear")
    plt.axis('off')

plt.show()

 

 

 

In [17]:

y_representative_digits = y_train[representative_digit_idx]
print(y_representative_digits.shape)

 

 

(50,)

 

In [18]:

log_reg = LogisticRegression(multi_class="ovr",
                             solver="lbfgs",
                             max_iter=5000,
                             random_state=42)
# 단 50개의 이미지로 정확도 92퍼를 달성했다
log_reg.fit(X_representative_digits, y_representative_digits)
log_reg.score(X_test, y_test)

 

Out[18]:

0.9222222222222223

 

 

그런데 여기서 더 나아가 대표 이미지의 레이블을 해당 클러스터에 속한 모든 샘플에 전파하면 어떨까?
이를 레이블 전파 (label propagation) 이라고 부른다.

 

In [19]:

y_train_propagated = np.empty(len(X_train), dtype=np.int32)
for i in range(k):
    y_train_propagated[kmeans.labels_ == i] = y_representative_digits[i]

log_reg = LogisticRegression(multi_class="ovr",
                             solver="lbfgs",
                             max_iter=5000,
                             random_state=42)
log_reg.fit(X_train, y_train_propagated)
log_reg.score(X_test, y_test)

 

Out[19]:

0.9333333333333333

 

 

성능이 좋아지긴 했지만, 경계에 있는 레이블들의 경우 에러가 발생했을 확률이 높다.
이번에는 센트로이드와 가까운 샘플의 15%에만 레이블을 전파해보자.

 

In [20]:

percentile_closest = 15

# 각 샘플마다 속한 클러스터에서의 센트로이드 까지의 거리가 담긴다.
X_cluster_dist = X_digits_dist[np.arange(len(X_train)), kmeans.labels_]

# 클러스터안에 속해있지만 15% 밖에 있는 샘플의 거리를 -1로 설정한다.
for i in range(k):
    in_cluster = (kmeans.labels_ == i)
    cluster_dist = X_cluster_dist[in_cluster]
    cutoff_distance = np.percentile(cluster_dist, percentile_closest)
    above_cutoff = (X_cluster_dist > cutoff_distance)
    X_cluster_dist[in_cluster & above_cutoff] = -1

# 거리가 -1을 제외한 모든 샘플에 레이블 전파를 적용한다.
partially_propagated = (X_cluster_dist != -1)
X_train_partially_propagated = X_train[partially_propagated]
y_train_partially_propagated = y_train_propagated[partially_propagated]

log_reg = LogisticRegression(multi_class="ovr",
                             solver="lbfgs",
                             max_iter=5000,
                             random_state=42)
log_reg.fit(X_train_partially_propagated, y_train_partially_propagated)
log_reg.score(X_test, y_test)

 

Out[20]:

0.9466666666666667

 

 

성능이 좋아짐을 알 수 있는데 이는 실제로 부분 전파된 샘플과 실제 데이터가 거의 일치하기 때문이다.

 

In [21]:

np.mean(y_train_partially_propagated==y_train[partially_propagated])

 

Out[21]:

0.9908675799086758

 

 

DBSCAN¶

밀집된 연속적인 지역을 하나의 클러스터로 정의하는 알고리즘이다.
우선 각 샘플에서 작은 거리인 $\epsilon$ 내에 샘플이 몇 개 놓여 있는지 카운팅 하고 이를 샘플의 $\epsilon -이웃$ 이라고 부른다.
$\epsilon -이웃$ 내에 적어도 min-samples개 샘플이 존재한다면 이를 핵심 샘플로 간주한다.
따라서 핵심 샘플은 밀집된 지역 내의 샘플이 된다.
핵심 샘플의 이웃에 있는 모든 샘플은 동일한 클러스터에 속한다. 따라서 이웃에는 다른 핵심 샘플이 포함될 수 있다.
그런데 핵심 샘플도, 이웃 샘플도 아닌 샘플은 이상치로 판단한다.

 

In [9]:

from sklearn.cluster import DBSCAN
from sklearn.datasets import make_moons

X, y = make_moons(n_samples=1000, noise=0.05)
dbscan = DBSCAN(eps=0.2, min_samples=5)
dbscan.fit(X)
print(dbscan.labels_)  # -1인 경우에는 이상치로 판단한 경우
print(dbscan.core_sample_indices_)  # 핵심 샘플의 인덱스
print(dbscan.components_)  # 핵심 샘플 자체가 담긴 인스턴스 변수

 

 

[0 0 0 1 1 0 1 1 0 0 1 1 0 1 0 1 1 1 1 1 0 0 1 1 1 1 0 0 1 0 0 0 0 0 0 1 1
 1 1 1 0 0 0 0 0 1 0 0 1 0 1 1 0 0 0 0 0 0 1 0 1 0 1 0 1 1 0 1 1 0 1 1 1 0
 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 0 0 0 1 1 1 1 0 1 0 1 0 0 0 1 1 0
 1 0 0 0 1 1 1 0 1 1 1 0 0 0 0 1 1 1 1 0 0 1 1 1 0 0 0 1 1 1 1 1 1 1 0 1 1
 0 1 1 1 1 0 1 1 0 0 1 0 1 1 0 1 0 0 1 0 1 1 1 1 1 1 1 0 1 0 1 0 0 0 1 1 0
 1 0 1 1 0 0 1 0 0 0 1 1 1 0 0 1 1 1 1 0 0 0 1 1 1 1 1 0 1 0 0 0 1 0 0 0 0
 1 0 0 0 1 0 1 0 0 0 0 0 0 0 1 1 0 1 1 0 1 0 0 0 0 1 0 1 0 0 1 0 0 0 1 0 0
 1 0 1 1 1 0 0 0 0 1 1 0 0 0 1 1 1 1 1 1 0 0 0 0 0 1 1 1 1 0 0 1 1 0 1 1 0
 0 0 0 1 1 0 0 0 1 0 1 0 0 0 1 1 0 1 1 1 1 1 0 0 0 1 1 1 0 1 1 1 0 1 1 0 0
 1 1 1 1 1 0 1 1 0 1 1 0 0 0 0 1 0 1 0 0 1 0 1 1 1 0 1 0 0 0 0 0 1 0 0 1 0
 0 1 1 0 1 1 1 0 0 0 0 0 1 0 0 0 0 0 1 0 1 0 0 0 1 1 0 1 1 0 1 0 0 0 0 1 1
 1 1 1 0 1 0 1 0 1 0 1 0 1 0 0 1 0 1 0 0 1 1 1 0 1 1 0 1 1 1 0 1 1 0 0 1 1
 0 1 0 1 1 0 1 1 0 1 0 1 0 1 0 1 1 0 0 1 1 0 1 1 1 1 0 0 1 1 0 0 1 1 0 0 0
 0 0 1 0 0 1 1 0 1 0 0 1 1 1 0 0 0 1 1 0 0 0 1 1 1 0 1 1 0 1 1 0 0 0 0 0 0
 0 0 0 1 1 1 1 0 1 1 0 0 0 0 1 0 0 0 1 1 1 1 1 0 0 1 1 0 0 0 1 1 0 0 0 0 1
 0 1 0 0 0 0 0 0 1 1 1 1 1 0 1 1 1 1 1 1 1 0 0 1 1 1 1 1 0 0 0 1 0 0 0 0 0
 0 0 1 1 1 1 0 0 1 0 0 0 1 0 1 1 0 1 1 0 0 0 1 0 0 0 1 1 0 0 1 0 1 1 0 1 1
 0 1 0 1 1 1 1 0 0 0 0 0 0 1 1 1 1 1 0 0 1 1 0 1 1 1 0 1 0 0 1 1 1 0 0 0 1
 1 0 1 1 0 0 0 0 1 1 0 1 0 0 0 0 1 0 1 0 1 0 0 0 0 1 1 1 0 1 1 1 0 0 0 1 0
 1 1 1 1 1 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1 0 1 0 1 1 0 0 0 1 0 0 1 1 1 1
 1 1 0 0 0 0 0 1 0 0 0 1 1 0 0 0 0 1 0 0 1 0 0 0 0 0 1 1 0 1 1 1 1 1 1 1 0
 0 0 0 1 0 0 0 1 1 1 1 1 1 0 1 1 1 1 1 0 0 0 1 0 1 0 0 0 0 1 0 1 0 1 1 0 1
 1 0 0 0 0 1 0 1 0 0 1 0 1 0 1 1 1 0 1 1 1 0 1 0 1 1 0 0 1 1 1 0 1 0 0 1 0
 0 1 0 0 0 0 1 1 0 0 1 0 0 1 1 0 1 1 0 1 0 0 1 0 1 0 1 1 1 1 0 1 0 1 1 0 1
 1 0 1 0 1 0 0 1 0 0 1 1 0 0 1 0 1 0 0 0 0 1 0 1 0 1 0 0 1 1 1 0 1 0 1 1 1
 0 0 0 1 0 1 1 1 0 0 0 0 0 1 0 1 0 1 0 0 1 0 1 1 1 1 1 0 1 0 1 0 0 1 1 0 1
 1 0 1 0 0 1 0 1 0 1 1 1 1 0 1 0 0 1 0 1 1 0 1 0 0 0 0 0 0 0 1 1 0 1 0 1 0
 0]
[  0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17
  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35
  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53
  54  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71
  72  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89
  90  91  92  93  94  95  96  97  98  99 100 101 102 103 104 105 106 107
 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125
 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143
 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161
 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179
 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197
 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215
 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233
 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251
 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269
 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287
 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305
 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323
 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341
 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359
 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377
 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395
 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413
 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431
 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449
 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467
 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485
 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503
 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521
 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539
 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557
 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575
 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593
 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611
 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629
 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647
 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665
 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683
 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701
 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719
 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737
 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755
 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773
 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791
 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809
 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827
 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845
 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863
 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881
 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899
 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917
 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935
 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953
 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971
 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989
 990 991 992 993 994 995 996 997 998 999]
[[ 0.80689722 -0.62325446]
 [ 0.50340641 -0.22042752]
 [ 0.11044464  0.08785997]
 ...
 [ 0.84847754  0.52598046]
 [ 1.76493483 -0.07692777]
 [ 0.61368154 -0.414085  ]]

 

 

불행히도 DBSCAN은 predict 매서드가 아닌 fit_predict 메서드를 제공한다.
즉 새로운 샘플에 대한 예측이 불가능하다. 따라서 예측을 진행할 예측기를 따로 선택해야 한다.

 

In [10]:

from sklearn.neighbors import KNeighborsClassifier

knn = KNeighborsClassifier(n_neighbors=50)
knn.fit(dbscan.components_, dbscan.labels_[dbscan.core_sample_indices_])
X_new = np.array([[-0.5, 0], [0, 0.5], [1, -0.1], [2, 1]])
print(knn.predict(X_new))
print(knn.predict_proba(X_new))
# k최근접이웃으로 핵심 샘플과 레이블을 학습시켜준 뒤 predict를 하면 된다.

 

 

[1 0 1 0]
[[0.14 0.86]
 [1.   0.  ]
 [0.22 0.78]
 [1.   0.  ]]

 

 

가우시안 혼합 모델 (GMM)¶

Gaussian Mixture Model은 샘플들이 파라미터가 알려지지 않은 여러 개의 혼합된 가우시안 분포에서 생성되었다고 가정하는 확률모델이다. 하나의 가우시안 분포에서 생성된 샘플은 같은 클러스터에 속하게 된다.
클러스터는 일반적으로 타원형이고 샘플이 주어지면 어떤 분포인지는 모르지만 가우시안 분포에서 나왔음을 알 수 있다.
(다만 어떤 분포이고 파라미터가 무엇인지는 모른다)

가우시안 혼합 모델이 작동되는 원리는 다음과 같다.
우선 샘플마다 k개의 클러스터에서 랜덤한 하나의 클러스터가 선택 된다.
j번째 클러스터를 선택할 확률은 $\Phi^{(j)}$로 정의되며 i번째 샘플이 선택한 레이블을 $z^{(i)}$라 하자.
$z^{(i)} = j$일 때, 즉 i번째 샘플이 j번째 클러스터에 할당됬다면 이 샘플의 위치 $x^{(i)}$ 는
평균이 $\mu^{(j)}$ 이고 공분산 행렬을 $\sum^{(j)}$ 로 가지는 가우시안 분포에서 랜덤하게 샘플링 된다.

이 모델을 이용한 추정은 데이터셋 X가 주어진다면 가중치 $\Phi$ 와 전체 분포의 파라미터인 $\mu, \sum$ 을 추정한다.

 

In [11]:

from sklearn.datasets import make_blobs
import numpy as np

X1, y1 = make_blobs(n_samples=1000, centers=((4, -4), (0, 0)), random_state=42)
X1 = X1.dot(np.array([[0.374, 0.95], [0.732, 0.598]]))
X2, y2 = make_blobs(n_samples=250, centers=1, random_state=42)
X2 = X2 + [6, -8]
X = np.r_[X1, X2]
y = np.r_[y1, y2]

 

In [12]:

from sklearn.mixture import GaussianMixture

gm = GaussianMixture(n_components=3, n_init=10)
gm.fit(X)
print(gm.weights_, '\n')
print(gm.means_, '\n')
print(gm.covariances_)

 

 

[0.3902064  0.20970929 0.40008431] 

[[ 0.0512132   0.07512971]
 [ 3.39854326  1.05926051]
 [-1.4076241   1.42704927]] 

[[[ 0.68780591  0.79597839]
  [ 0.79597839  1.21232694]]

 [[ 1.14901133 -0.03241901]
  [-0.03241901  0.95484861]]

 [[ 0.63480459  0.72971575]
  [ 0.72971575  1.16110086]]]

 

 

사이킷런의 GaussianMixture 클래스는 클러스터를 찾기 위하여 기댓값-최대화 (Expectation-Maximization)을 이용한다.
k-평균과 비슷하게 우선 샘플을 클러스터에 할당하고 각 클러스터 파라미터에 대하여 클러스터에 속할 확률을 예측한다.
(기댓값 단계, Expectation Step)
그런 다음 각 클러스터가 데이터셋에 있는 모든 샘플을 이용하여 업데이트된다. (최대화 단계, Maximization Step)
이때, 기댓값 단계에서 구한 클러스터에 속할 확률이 최대화 단계에서 샘플의 가중치로 적용되고
이를 샘플에 대한 클러스터에 대한 책임이라고 부른다.
따라서 최대화 단계에서 클러스터의 업데이트는 책임이 가장 큰 샘플에게 크게 영향을 받게 된다.
이런 식으로 클러스터의 중심, 모양, 상대적 가중치를 찾아가는 k-평균의 일종이라고 생각 가능하다.
이렇게 클러스터의 평균, 공분산, 각 샘플에 대한 확률을 점진적으로 개선해가며 수렴할 때 까지 진행한다.
만약 수렴이 된다면 각 샘플이 클러스터에 속할 확률이 나올 것이고 클러스터링 해주는 것이 가능하다.
다만, 횟수가 적을 수록 k-평균과 마찬가지로 나쁜 솔루션으로 수렴 가능하니 n_init을 늘려서 실행회수를 늘리자.

 

In [13]:

# 우선 가우시안 혼합 모델이 수렴됬는지 확인하자
print(gm.converged_,'\n')

# 언제 수렴 됬는지?
print(gm.n_iter_,'\n')

# predict도 당연히 있다.
print(gm.predict(X),'\n')

# predict_proba도 당연히 있다.
print(gm.predict_proba(X),'\n')

# 심지어 가우시안 모델은 생성 모델이기에 이 모델에서 새로운 샘플을 만들 수 있다.
X_new, y_new = gm.sample(6)
print(X_new,'\n')
print(y_new)

 

 

True 

4 

[0 0 2 ... 1 1 1] 

[[9.76688618e-01 2.33107018e-02 6.80230786e-07]
 [9.82774531e-01 1.65481663e-02 6.77302883e-04]
 [7.42085667e-05 2.04358925e-06 9.99923748e-01]
 ...
 [4.21582331e-07 9.99999578e-01 2.17682412e-26]
 [4.93027253e-16 1.00000000e+00 1.50345348e-41]
 [2.19825399e-15 1.00000000e+00 8.24833927e-41]] 

[[-0.24935178 -0.8012166 ]
 [-0.76823089 -0.98808634]
 [ 2.13733683  3.64772362]
 [ 3.83013019  0.56550826]
 [ 3.75281167  1.58978153]
 [-2.20685621 -0.10786302]] 

[0 0 1 1 1 2]

 

 

주어진 샘플에 대하여 해당 위치에서의 모델의 밀도를 계산하는 것이 가능하다.
이를 위하여 score_samples() 메서드를 이용하면 그 위치에 대한
확률 밀도 함수 (Probability Density Function)의 로그를 예측한다.
실제 값을 얻고 싶다면 지숫값으로 계산하면 되고 이 값으로 등고선을 그리는 것이 가능하다.

 

In [14]:

gm.score_samples(X)

 

Out[14]:

array([-2.60753797, -3.57117632, -3.32962239, ..., -3.51337454,
       -4.39800533, -3.80760349])

 

In [15]:

from matplotlib.colors import LogNorm

def plot_centroids(centroids, weights=None, circle_color='w', cross_color='k'):
    if weights is not None:
        centroids = centroids[weights > weights.max() / 10]
    plt.scatter(centroids[:, 0], centroids[:, 1],
                marker='o', s=35, linewidths=8,
                color=circle_color, zorder=10, alpha=0.9)
    plt.scatter(centroids[:, 0], centroids[:, 1],
                marker='x', s=2, linewidths=12, 
                color=cross_color, zorder=11, alpha=1)
    
def plot_gaussian_mixture(clusterer, X, resolution=1000, show_ylabels=True):
    mins = X.min(axis=0) - 0.1
    maxs = X.max(axis=0) + 0.1
    xx, yy = np.meshgrid(np.linspace(mins[0], maxs[0], resolution),
                         np.linspace(mins[1], maxs[1], resolution))
    Z = -clusterer.score_samples(np.c_[xx.ravel(), yy.ravel()])
    Z = Z.reshape(xx.shape)

    plt.contourf(xx, yy, Z,
                 norm=LogNorm(vmin=1.0, vmax=30.0),
                 levels=np.logspace(0, 2, 12))
    plt.contour(xx, yy, Z,
                norm=LogNorm(vmin=1.0, vmax=30.0),
                levels=np.logspace(0, 2, 12),
                linewidths=1, colors='k')

    Z = clusterer.predict(np.c_[xx.ravel(), yy.ravel()])
    Z = Z.reshape(xx.shape)
    plt.contour(xx, yy, Z,
                linewidths=2, colors='r', linestyles='dashed')
    
    plt.plot(X[:, 0], X[:, 1], 'k.', markersize=2)
    plot_centroids(clusterer.means_, clusterer.weights_)

    plt.xlabel("$x_1$", fontsize=14)
    if show_ylabels:
        plt.ylabel("$x_2$", fontsize=14, rotation=0)
    else:
        plt.tick_params(labelleft=False)
plt.figure(figsize=(8, 4))

plot_gaussian_mixture(gm, X)

plt.show()

 

 

 

 

실제 데이터셋에서는 가우시안 분포나 저차원이 아닌 경우가 많고, 정확한 값들을 지정해주었다.
이런 경우 EM 알고리즘이 최적의 솔루션을 찾기가 어려우므로 파라미터의 개수를 제한해야 한다.
예를 들면 클러스터의 모양과 방향, 범위를 공분산 행렬을 통해 제약을 걸 수 있다.

 

 

가우시안 혼합을 이용한 이상치 탐지¶

이상치 탐지는 기존과 다른 샘플을 탐지하는 과정이다. 이를 이상치라고 부르고 보톰 샘플을 정상치라고 부른다.
가우시안 혼합 모델에서는 밀도가 낮은 지역에 있는 모든 샘플을 모두 이상치로 볼 수 있긴 하다.
(즉, 임계 밀도 이하에 존재하는 모든 샘플을 이상치로 판단하는 것이다.)
너무 거짓양성이 많이 나온다면 임계밀도를 낮추고 거짓음성이 많이 나온다면 임계밀도를 임계밀도를 높여주자.

 

In [16]:

densities = gm.score_samples(X)
density_threshold = np.percentile(densities, 4) # 밀도가 4퍼보다 낮다면 이상치가 된다.
anomalies = X[densities < density_threshold]
anomalies.shape # 총 50개의 샘플이 이상치로 판명되었다.

 

Out[16]:

(50, 2)

 

 

이와 유사하게 특이치 탐지가 있다. 특이치는 지금 까지 훈련된 데이터에 이상치가 없다는 것을 전제로 한다.
그 이후, 들어오는 데이터가 지금까지의 데이터와 유사한지 (같은 클러스터일지)를 따지는 문제이다.
반면 이상치 탐지는 이런 전제 없이 데이터셋에 이상한 값이 있는 지를 판단하고 정제하는 것을 목표로 한다.

 

 

베이즈 가우시안 혼합 모델¶

최적의 클러스터 개수를 알아서 찾아주는 베이즈 가우시안 모델.
이 알고리즘은 불필요한 클러스터를 제거해주고 클러스터 파라미터를 고정된 모델 파라미터가 아닌 확률 변수로 취급한다.

 

In [17]:

from sklearn.mixture import BayesianGaussianMixture
bgm = BayesianGaussianMixture(n_components=10, n_init=10, random_state=42)
bgm.fit(X)

np.round(bgm.weights_, 2) # 3개의 클러스터가 필요함이 나왔다.

 

Out[17]:

array([0.4 , 0.21, 0.4 , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  ])