Data Science/데이터 분석 📊
[데이터 분석] 20. 차원 축소 (Dimensional Reduction) Ⅱ : 커널 PCA, LLE, LDA
SLYK1D
2024. 10. 7. 14:51
728x90
반응형
1. 커널 PCA
이전 SVM에 대한 내용 중 샘플을 매우 높은 고차원의 특성 공간에서 암묵적으로 매핑해 SVM의 비선형 분류와 회귀를 가능하게 하는 커널트릭에 대해서 다뤘적이 있는데, 고차원의 특성 공간에서의 선형 결정경계는 원본 공간에서의 복잡한 비선형 결정경계에 해당하는데, 이를 나눌 수 있는 기법을 커널 트릭이라고 설명했다.
PCA에서도 커널트릭을 적용해 복잡한 비선형 투영으로의 차원축소가 가능하며, 이를 가리켜, kPCA(커널 PCA, kernel PCA) 라고 부른다. 이 기법은 이전 장에서 봤던 스위스롤과 같은 매니폴드 형태의 데이터 셋을 대상으로 투영한 뒤, 샘플의 군집을 유지하거나, 펼칠 때 주로 사용된다. 이해를 돕기 위해 아래 코드를 실행해서 살펴보자.
...
from sklearn.decomposition import KernelPCA
...
rbf_pca = KernelPCA(n_components=2, kernel='rbf', gamma=0.04)
X_reduced = rbf_pca.fit_transform(X)위와 같은 방법으로 앞서 본 스위스롤 데이터에 비선형 투영을 적용하면, 아래와 같은 결과 그래프를 얻을 수 있다. 아래 그래프의 왼쪽에서부터 각각 선형 커널, RBF 커널, 시그모이드 커널을 적용한 결과이다.
커널 PCA 는 비지도학습으로 분류되기 때문에, 좋은 커널과 하이퍼파라미터를 선택하기 위해 명확한 성능 측정 기준이 없다. 하지만, 지도학습의 전처리 과정으로 사용되기 때문에 그리드 탐색을 통해 최적의 커널과 하이퍼파라미터를 선택할 수 있다. 해당 내용을 확인해보기 위해 앞서 사용했던 스위스롤 데이터를 다시 사용할 예정이며, 작업은 크게 2단계의 파이프라인으로 구성된다.
먼저, 1단계에서는 커널 PCA를 사용해 3차원의 원본데이터를 2차원으로 축소하고, 분류를 하기 위해 로지스틱 회귀를 선택한다. 2단계에서는 가장 높은 분류 정확도를 얻기 위해서 GridSearchCV 알고리즘을 사용하며, 커널 PCA의 가장 좋은 커널과 파라미터 중 하나인 Gamma 값을 산출한다. 끝으로, 가장 성능이 좋은 커널 및 하이퍼파라미터는 GridSearch 모델에서 제공하는 최적값으로 변수에 저장한다.
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
from sklearn.model_selection import GridSearchCV
from sklearn.linear_model import LogisticRegression
from sklearn.pipeline import Pipeline
from sklearn.decomposition import KernelPCA
from sklearn.datasets import make_swiss_roll
X, t = make_swiss_roll(n_samples=1000, noise=0.2, random_state=42)
y = t > 6.9
clf = Pipeline([
("kpca", KernelPCA(n_components=2)),
("log_reg", LogisticRegression(solver='liblinear'))
])
param_grid = [{
"kpca__gamma": np.linspace(0.03, 0.05, 10),
"kpca__kernel": ["rbf", "sigmoid"]
}]
grid_search = GridSearchCV(clf, param_grid, cv=3)
grid_search.fit(X, y)
print(grid_search.best_params_)# 필요한 라이브러리
library(ggplot2)
library(caret)
library(kernlab) # For KernelPCA
library(MASS) # For logistic regression
library(gridExtra)
# 스위스 롤 데이터셋 생성 (Python의 make_swiss_roll에 해당하는 코드)
set.seed(42)
n_samples <- 1000
t <- 3 * pi / 2 * (1 + 2 * runif(n_samples))
height <- 21 * runif(n_samples)
X1 <- t * cos(t)
X2 <- height
X3 <- t * sin(t)
X <- data.frame(X1 = X1, X2 = X2, X3 = X3)
y <- t > 6.9
# 1. KernelPCA와 Logistic Regression 파이프라인 구축
# 커널 PCA 수행 함수 정의
kernel_pca_logistic_pipeline <- function(gamma, kernel_function) {
kpca_result <- kpca(~., data = X, kernel = kernel_function, kpar = list(sigma = gamma), features = 2)
X_kpca <- as.data.frame(pcv(kpca_result)) # 변환된 데이터를 데이터프레임으로 변환
# Kernel PCA로 변환된 데이터에 y 추가
X_kpca$y <- as.factor(y)
# Logistic Regression
log_reg_model <- glm(y ~ ., data = X_kpca, family = binomial)
return(list(model = log_reg_model, X_kpca = X_kpca))
}
# 2. 그리드 서치 준비
grid <- expand.grid(gamma = seq(0.03, 0.05, length.out = 10), kernel = c("rbfdot", "tanhdot"))
# 3. 그리드 서치를 통해 최적의 하이퍼파라미터 찾기
best_model <- NULL
best_accuracy <- -Inf
best_params <- NULL
for (i in 1:nrow(grid)) {
gamma <- grid$gamma[i]
kernel_name <- as.character(grid$kernel[i]) # factor를 문자열로 변환
# 커널 함수를 이름에 맞게 가져옴
kernel_function <- switch(kernel_name,
"rbfdot" = rbfdot(),
"tanhdot" = tanhdot())
# 모델 훈련
result <- kernel_pca_logistic_pipeline(gamma, kernel_function)
log_reg_model <- result$model
X_kpca <- result$X_kpca
# 교차 검증
cv_results <- train(y ~ .,
data = X_kpca,
method = "glm",
trControl = trainControl(method = "cv", number = 3),
family = "binomial"
)
accuracy <- max(cv_results$results$Accuracy)
if (accuracy > best_accuracy) {
best_accuracy <- accuracy
best_model <- log_reg_model
best_params <- list(gamma = gamma, kernel = kernel_name)
}
}
# 최적의 하이퍼파라미터 출력
print(best_params)
[실행 결과]
$gamma
[1] 0.03222222
$kernel
[1] "rbfdot"위에서 설명한 과정 이외에도, 완전 비지도학습 방법으로 가장 낮은 재구성 오차를 만드는 커널과 하이퍼파라미터를 선택하는 방법도 있으나, PCA 만큼 쉽지 않다. 축소된 공간에 있는 샘플에 대해 선형 PCA를 역전시키면, 재구성된 데이터 포인트들은 원본공간이 아닌 특성공간에 놓이게 된다.
특성공간은 무한 차원이기에, 재구성된 데이터 포인트들을 계산할 수 없고, 재구성에 따른 실제 오차를 계산할 수 없다. 하지만, 재구성된 포인트에 가깝게 매핑된 원본공간의 포인트를 찾을 수 있으며, 이러한 포인트들을 가리켜 재구성 원상이라고 부른다. 원상을 산출하게되면, 원본 샘플과의 제곱거리를 측정할 수 있으며, 결과적으로 재구성 원상과의 오차를 최소화하는 커널과 하이퍼파라미터를 선택할 수 있게 된다.
재구성을 하는 방법으로는 투영된 샘플을 학습 데이터 셋으로, 원본 샘플을 타겟으로 하여, 지도학습 중 회귀모델을 훈련시키는 방법이다.
위의 내용을 코드로 구현해보자면, 아래와 같이 작성할 수 있다.
rbf_pca = KernelPCA(n_components=2, kernel="rbf", gamma=0.0433, fit_inverse_transform=True)
X_reduced = rbf_pca.fit_transform(X)
X_preimage = rbf_pca.inverse_transform(X_reduced)# 필요한 라이브러리
library(kernlab)
rbf_kpca <- kpca(~., data = X, kernel = "rbfdot", kpar = list(sigma = 0.0433), features = 2, fit.inverse = TRUE)
X_reduced <- pcv(rbf_kpca)
print("Reduced data (2D):")
print(head(X_reduced))
X_preimage <- fitted(rbf_kpca)
print("Preimage (Recovered data):")
print(head(X_preimage))위의 코드에서 커널PCA 모델을 생성할 때, scikit-learn 의 경우, fit_inverse_transform 파라미터 값을 True로 설정해야 모델의 inverse_transform() 메소드를 사용할 수 있으며, 만약 설정 없이 실행하면, 아래와 같은 에러 문구가 나오니 주의바란다. 참고로, fit_inverse_transform 옵션의 기본값을 false 이다.
sklearn.exceptions.NotFittedError: The fit_inverse_transform parameter was not set to True when instantiating and hence the inverse transform is not available.
끝으로 재구성 원상 오차를 계산하여 확인해보자.
from sklearn.metrics import mean_squared_error
print("MSE between original to prediction : ", mean_squared_error(X, X_preimage))
[실행 결과]
MSE between original to prediction : 32.78630879576615
2. 지역 선형 임베딩(LLE, Locally Linear Embedding)
비선형 차원축소 기법으로 투영에 의존하지 않는 매니폴드 학습 기법이다. 이전 알고리즘들과의 차이는 각 학습 데이터 샘플이 가장 가까이 위치한 이웃과 얼마나 선형적으로 연관있는가를 측정한다. 특히, 노이즈가 심하지 않은 경우의 꼬인 매니폴드 형태를 펼치는 데 잘 작동한다.구체적인 과정은 다음과 같다.
① 각 학습 데이터 샘플 $ {X}^{(i)} $ 에 대해 가장 가까운 샘플 K개를 찾는다.
② 이웃에 대한 선형함수로 $ {X}^{(i)} $ 를 재구성한다.
③ 이 때, $ {X}^{(i)} $ 와 $ \sum_{j=1}^{m} {W_{i,j} X^{(i)}} $ 사이의 제곱거리가 최소가 되도록 하는 $ W_{i,j} $ 를 찾는다.
만약, ${X}^{(j)} $ 가 $ {X}^{(i)} $ 의 가장 가까운 이웃 k 개에 포함되지 않으면, $ W_{i,j} = 0 $ 이 된다. 또한 각 학습
데이터 샘플 $ {X}^{(i)} $ 에 대한 가중치는 정규화가 되어있어야 한다.
④ 앞선 단계를 거치면, 가중치 행렬 W 는 학습 데이터 샘플 사이에 있는 지역 선형 관계를 갖게 된다. 따라서 관계가 보존되도록
학습 데이터 샘플을 d차원 공간으로 매핑한다.
위의 알고리즘 과정을 아래 코드를 통해서 확인해보도록 하자.
from sklearn.manifold import LocallyLinearEmbedding
lle = LocallyLinearEmbedding(n_components=2, n_neighbors=10)
X_reduced = lle.fit_transform(X)library(lle)
X_reduced <- lle(X, m = 2, k = 10)
print("Reduced data (2D):")
print(head(X_reduced$Y))3. LDA (Linear Discriminant Analysis)
선형 판별 분석이라고도 부르는 LDA 기법은 주로 규제가 없는 모델에서 차원의 저주로 인해 과대적합을 줄이고, 계산 효율성을 높이기 위해서 사용되는 특성 추출 기법이다. PCA와 매우 유사하지만, 지도학습의 일종이며, 분류 작업에서 뛰어난 성능을 보여준다. 또한 데이터가 정규분포를 가진다고 가정하여, 클래스가 동일한 공분산 행렬을 갖는 샘플은 서로 통계적으로 독립적이라고 가정한다. 구체적인 알고리즘의 동작과정은 아래와 같다.
① d 차원의 데이터 셋을 표준화 전처리한다.
② 각 클래스에 대해 d 차원의 평균 벡터를 계산한다.
③ 클래스 간의 산포행렬인 $ S_b $ 와 클래스 내 산포 행렬인 $ S_w $ 로 구성된다.
④ $ S_b $ 와 $ S_w $ 행렬의 고유 벡터와 고유 값을 계산한다.
⑤ 고유 값을 내림차순으로 정렬하여, 고유 벡터의 순서를 매긴다.
⑥ 고유 값이 가장 큰 k개의 고유 벡터를 선택하여 d x k 차원의 변환 행렬인 W 를 구성한다.
⑦ 변환 행렬 W 를 사용해 샘플의 새로운 특성 부분 공간으로 투영한다.
3.1 산포 행렬 계산
앞서 알고리즘 과정에서 설명했듯이, 선형 판별 분석에서 중요한 부분 중 하나인 산포 행렬을 어떻게 생성하는지 좀 더 살펴보도록 하자. 과정에서 이야기했듯, 클래스 간의 산포행렬과 클래스 내의 산포행렬은 모두 평균 벡터를 통해서 계산되는데, 이 때 평균 벡터는 클래스 i 의 샘플에 대한 특성의 평균값을 저장한다.
위의 내용을 구체적인 예시를 통해서 확인해보도록 하자. 실습에 사용할 데이터는 앞선 예시로 사용된 와인 데이터를 사용할 것이며, 데이터 표준화 부분까지 동일한 과정으로 진행되기에 생략하도록 하겠다. 클래스는 1~3까지 총 3개로 구성되며, 각 클래스별 평균 벡터는 아래 코드와 같이 계산할 수 있다.
# 실행전 y 의 label 갯수 확인
y_label = set(y_train)
print(y_label) # 1 ~ 3
np.set_printoptions(precision=4) # 부동소수점, Array, 기타 numpy 객체가 표시되는 방식을 설정함
mean_vecs =[]
for label in range(1, 4):
mean_vecs.append(np.mean(x_train_std[y_train == label], axis=0)) # 각 label 별 평균값을 계산
print("Mean Vector %s : %s\n" % (label, mean_vecs[label-1]))
[실행 결과]
Mean Vector 1 : [ 0.9066 -0.3497 0.3201 -0.7189 0.5056 0.8807 0.9589 -0.5516 0.5416
0.2338 0.5897 0.6563 1.2075]
Mean Vector 2 : [-0.8749 -0.2848 -0.3735 0.3157 -0.3848 -0.0433 0.0635 -0.0946 0.0703
-0.8286 0.3144 0.3608 -0.7253]
Mean Vector 3 : [ 0.1992 0.866 0.1682 0.4148 -0.0451 -1.0286 -1.2876 0.8287 -0.7795
0.9649 -1.209 -1.3622 -0.4013]library(dplyr)
# y_train의 레이블 갯수 확인
y_label <- unique(y_train)
print(y_label) # 1, 2, 3
# 소수점 출력 설정 (Python의 np.set_printoptions(precision=4)와 동일)
options(digits = 4)
# 각 클래스별 평균 벡터 계산
mean_vecs <- list()
for (label in 1:3) {
# 각 클래스별 평균 계산
mean_vec <- colMeans(x_train_std[y_train == label, ])
mean_vecs[[label]] <- mean_vec
cat(sprintf("Mean Vector %d : %s\n", label, toString(mean_vecs[[label]])))
}산포 행렬은 평균 벡터를 이용해 아래의 식으로 표현할 수 있다.
$ S_W = \sum_{i=1}^{C} {S_i} $
위의 식을 코드로 구현하면 다음과 같이 표현할 수 있다.
d = 13 # 특성 계수
Sw = np.zeros((d, d))
for label, mv in zip(range(1, 4), mean_vecs):
class_scatter = np.zeros((d, d))
for row in x_train_std[y_train == label]:
row, mv = row.reshape(d, 1), mv.reshape(d, 1)
class_scatter += (row - mv).dot((row - mv).T)
Sw += class_scatter
print("클래스 내의 산포 행렬 : %s x %s" % (Sw.shape[0], Sw.shape[1]))
print('클래스 레이블 분포 : %s' % np.bincount(y_train)[1:])
[실행 결과]
클래스 내의 산포 행렬 : 13 x 13
클래스 레이블 분포 : [41 50 33]# 특성 계수 (d)
d <- ncol(x_train_std)
# 클래스 내의 산포 행렬 초기화
Sw <- matrix(0, nrow = d, ncol = d)
# 클래스별 산포 행렬 계산
for (label in 1:3) {
mv <- mean_vecs[[label]]
class_scatter <- matrix(0, nrow = d, ncol = d)
# 해당 레이블에 대한 행만 추출
label_rows <- x_train_std[y_train == label, ]
# 각 샘플에 대해 산포 행렬 계산
for (i in 1:nrow(label_rows)) {
row <- as.matrix(label_rows[i, ]) # 1 x d 형태로 변환
mv_matrix <- matrix(mv, ncol = 1) # d x 1 형태로 변환
class_scatter <- class_scatter + (row - mv_matrix) %*% t(row - mv_matrix)
}
Sw <- Sw + class_scatter
}
# 클래스 내의 산포 행렬 크기 출력
cat(sprintf("클래스 내의 산포 행렬 : %d x %d\n", nrow(Sw), ncol(Sw)))
# 클래스 레이블 분포 확인
class_distribution <- table(y_train)
cat(sprintf("클래스 레이블 분포 : %s\n", toString(class_distribution)))개별 산포행렬인 $ S_i $ 를 산포행렬 $ S_W $ 로 모두 더하기 전에 스케일 조정을 해야한다. 산포행렬을 클래스 샘플 개수로 나누면, 산포 행렬을 계산하는 것이 공분산 행렬을 계산하는 것과 동일해진다. 결과적으로 공분산 행렬은 산포행렬의 정규화한 결과라고 표현할 수 있다.
Sw = np.zeros(((d, d)))
for label, mv in zip(range(1, 4), mean_vecs):
class_scatter = np.cov(x_train_std[y_train == label].T, bias=True)
Sw += class_scatter
print("스케일 조정된 클래스 내의 산포 행렬 : %s x %s" % (Sw.shape[0], Sw.shape[1]))
[실행 결과]
스케일 조정된 클래스 내의 산포 행렬 : 13 x 13# 특성 계수 (d)
d <- ncol(x_train_std)
# 클래스 내의 산포 행렬 초기화
Sw <- matrix(0, nrow = d, ncol = d)
# 클래스별 공분산 행렬 계산 및 합산
for (label in 1:3) {
mv <- mean_vecs[[label]]
# 해당 레이블에 대한 데이터 추출
label_rows <- x_train_std[y_train == label, ]
# 클래스별 공분산 행렬 계산 (bias = TRUE는 편향된 추정치로 계산하는 옵션)
class_scatter <- cov(label_rows, use = "everything") * (nrow(label_rows) - 1) / nrow(label_rows)
# 클래스 산포 행렬 합산
Sw <- Sw + class_scatter
}
# 스케일 조정된 클래스 내의 산포 행렬 크기 출력
cat(sprintf("스케일 조정된 클래스 내의 산포 행렬 : %d x %d\n", nrow(Sw), ncol(Sw)))다음으로 클래스 간 산포행렬 인 $S_B$ 를 계산한다.
mean_class = np.mean(x_train_std, axis=0)
mean_class = mean_class.reshape(d, 1)
S_B = np.zeros((d, d))
for i, mean_vec in enumerate(mean_vecs):
n = x_train[y_train == i + 1, :].shape[0]
mean_vec = mean_vec.reshape(d, 1)
S_B += n * (mean_vec - mean_class).dot((mean_vec - mean_class).T)
print("클래스 간 산포행렬 : %s x %s" % (S_B.shape[0], S_B.shape[1]))
[실행결과]
클래스 간 산포행렬 : 13 x 13# 특성 계수 (d)
d <- ncol(x_train_std)
# 전체 데이터의 평균 계산
mean_class <- colMeans(x_train_std)
mean_class <- matrix(mean_class, nrow = d, ncol = 1) # d x 1 형태로 변환
# 클래스 간 산포 행렬 초기화
S_B <- matrix(0, nrow = d, ncol = d)
# 클래스 간 산포 행렬 계산
for (i in 1:3) {
mean_vec <- matrix(mean_vecs[[i]], nrow = d, ncol = 1) # 클래스별 평균 벡터를 d x 1 형태로 변환
n <- sum(y_train == i) # 해당 클래스에 속하는 샘플의 개수
# S_B 계산
S_B <- S_B + n * (mean_vec - mean_class) %*% t(mean_vec - mean_class)
}
# 클래스 간 산포 행렬 크기 출력
cat(sprintf("클래스 간 산포 행렬 : %d x %d\n", nrow(S_B), ncol(S_B)))3.2 선형 판별 벡터 생성 및 투영
이 후 과정은 PCA에서와 유사하게 행렬 $ {S}_{W}^{-1} \cdot {S}_{B} $ 를 계산해주면 된다.
eigen_vals, eigen_vecs = np.linalg.eig(np.linalg.inv(Sw).dot(S_B))
eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:, i]) for i in range(len(eigen_vals))]
eigen_pairs = sorted(eigen_pairs, key=lambda k: k[0], reverse=True)
print("내림차순 고유값\n")
for eigen_val in eigen_pairs:
print(eigen_val[0])
[실행결과]
내림차순 고유값
358.00420701336594
177.07768640666228
2.897788744761761e-14
2.842170943040401e-14
1.9861923711452444e-14
1.969800906168211e-14
1.969800906168211e-14
1.8105958343606454e-14
1.4708272757021412e-14
8.653798256529554e-15
8.653798256529554e-15
3.4503200980472656e-15
8.922969632820571e-17# 1. 고유값 및 고유벡터 계산
eigen_result <- eigen(solve(Sw) %*% S_B)
eigen_vals <- eigen_result$values # 고유값
eigen_vecs <- eigen_result$vectors # 고유벡터
# 2. 고유값과 고유벡터를 쌍으로 묶기
eigen_pairs <- lapply(1:length(eigen_vals), function(i) {
list(abs(eigen_vals[i]), eigen_vecs[, i])
})
# 3. 고유값을 기준으로 내림차순 정렬
eigen_pairs <- eigen_pairs[order(sapply(eigen_pairs, function(x) x[[1]]), decreasing = TRUE)]
# 4. 내림차순 고유값 출력
cat("내림차순 고유값\n")
for (eigen_val in eigen_pairs) {
cat(eigen_val[[1]], "\n")
}LDA에서 선형판별벡터의 길이는 "클래스 수 - 1" 개 이다. 클래스 내 산포행렬인 $ S_B $ 가 랭크 1 또는 그 이하인 클래스의 수만큼의 행렬을 합한 것이기 때문이다. 이를 그래프로 표현해보자면, 아래와 같은 형태를 볼 수 있다.
끝으로, 2개의 판별 고유 벡터를 열로 쌓은 변환 벡터 W를 생성하고, 이를 학습 데이터를 곱해 데이터를 변환시켜보자.
W = np.hstack((eigen_pairs[0][1][:, np.newaxis].real,
eigen_pairs[1][1][:, np.newaxis].real))
print("행렬 W:\n", W)
x_train_lda = x_train_std.dot(W)
colors = ['r', 'b', 'g']
markers = ['s', 'x', 'o']
for l, c, m in zip(np.unique(y_train), colors, markers):
plt.scatter(x_train_lda[y_train==l, 0],
x_train_lda[y_train==l, 1] * (-1),
c=c, label=l, marker=m)
plt.xlabel("LD 1")
plt.ylabel("LD 2")
plt.legend(loc="lower right")
plt.tight_layout()
plt.show()# 필요한 라이브러리
library(ggplot2)
# 고유값과 고유벡터가 이미 계산되었다고 가정 (eigen_pairs 리스트 사용)
# 1. W 행렬 생성 (두 개의 고유벡터를 선택)
W <- cbind(eigen_pairs[[1]][[2]], eigen_pairs[[2]][[2]])
cat("행렬 W:\n")
print(W)
# 2. LDA로 차원 축소
x_train_lda <- as.matrix(x_train_std) %*% W
# 3. y_train을 factor로 변환 (R에서는 범주형 데이터로 취급해야 함)
y_train_factor <- as.factor(y_train)
# 4. 차원 축소된 데이터 시각화
lda_data <- data.frame(LD1 = x_train_lda[, 1],
LD2 = x_train_lda[, 2] * (-1), # 축소된 두 번째 축을 반전
Label = y_train_factor)
# 색상과 마커 설정
colors <- c('r', 'b', 'g')
markers <- c(16, 17, 18)
# 시각화
ggplot(lda_data, aes(x = LD1, y = LD2, color = Label, shape = Label)) +
geom_point(size = 3) +
scale_color_manual(values = colors) +
scale_shape_manual(values = markers) +
labs(x = "LD 1", y = "LD 2") +
theme_minimal() +
theme(legend.position = "lower right")
추가적으로 위의 모든 과정에 대해 파이썬의 scikit-learn 에서는 lda() 함수를 통해 단순하면서도 쉽게 계산하는 것이 가능하다. 확인을 위해 앞서 변환한 데이터를 이용해 로지스틱 회귀로 예측해보도록 하자.
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
from sklearn.linear_model import LogisticRegression
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.font_manager as fm
from matplotlib.colors import ListedColormap
lda = LDA(n_components=2)
x_train_lda = lda.fit_transform(x_train_std, y_train)
def plot_decision_regions(x, y, classifier, test_idx=None, resolution=0.02):
markers = ('s', 'x', 'o', '^', 'v')
colors = ("red", "blue", "lightgreen", "gray", "cyan")
cmap = ListedColormap(colors[:len(np.unique(y))])
x1_min, x1_max = x[:, 0].min() - 1, x[:, 0].max() + 1
x2_min, x2_max = x[:, 1].min() - 1, x[:, 1].max() + 1
xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, resolution),
np.arange(x2_min, x2_max, resolution))
z = classifier.predict(np.array([xx1.ravel(), xx2.ravel()]).T)
z = z.reshape(xx1.shape)
plt.contourf(xx1, xx2, z, alpha=.3, cmap=cmap)
plt.xlim(xx1.min(), xx1.max())
plt.ylim(xx2.min(), xx2.max())
for idx, cl in enumerate(np.unique(y)):
plt.scatter(x = x[y == cl, 0], y = x[y == cl, 1],
alpha=0.8, c=colors[idx], marker=markers[idx],
label=cl, edgecolors="black")
if test_idx:
x_test, y_test = x[test_idx, :], y[test_idx]
plt.scatter(x_test[:, 0], x_test[:, 1],
c='', edgecolors="black", alpha=1.0,
linewidth=1, marker="o",
s=100, label="test set")
lr = LogisticRegression(solver="liblinear", multi_class="auto")
lr = lr.fit(x_train_lda, y_train)
plot_decision_regions(x_train_lda, y_train, classifier=lr)# 필요한 라이브러리
library(MASS) # LDA
library(ggplot2)
library(caret) # Logistic regression
library(RColorBrewer)
# 1. LDA 수행
lda_model <- lda(y_train ~ ., data = as.data.frame(x_train_std)) # y_train을 종속 변수로 하여 LDA 수행
x_train_lda <- as.matrix(predict(lda_model)$x) # 차원 축소된 LDA 결과
# 2. 로지스틱 회귀 모델 학습
lr_model <- train(as.factor(y_train) ~ ., data = as.data.frame(x_train_lda), method = "glm", family = "binomial")
# 3. 결정 경계 시각화 함수 정의
plot_decision_regions <- function(x, y, model, resolution = 0.02) {
colors <- c("red", "blue", "lightgreen", "gray", "cyan")
markers <- c(16, 17, 18, 19, 20)
# 결정 경계 설정
x1_min <- min(x[, 1]) - 1
x1_max <- max(x[, 1]) + 1
x2_min <- min(x[, 2]) - 1
x2_max <- max(x[, 2]) + 1
# 그리드 생성
xx1 <- seq(x1_min, x1_max, by = resolution)
xx2 <- seq(x2_min, x2_max, by = resolution)
grid <- expand.grid(xx1, xx2)
colnames(grid) <- c("LD1", "LD2")
# 예측
z <- predict(model, newdata = grid, type = "raw")
# 예측 결과를 그리드 형태로 변환
grid$z <- as.numeric(as.factor(z))
# 결정 경계 플롯 생성
ggplot() +
geom_tile(data = grid, aes(x = LD1, y = LD2, fill = as.factor(z)), alpha = 0.3) +
scale_fill_manual(values = colors[1:length(unique(y))]) +
geom_point(data = as.data.frame(x), aes(x = LD1, y = LD2, color = as.factor(y), shape = as.factor(y)), size = 3) +
scale_color_manual(values = colors) +
scale_shape_manual(values = markers) +
labs(x = "LD 1", y = "LD 2", title = "LDA Decision Regions") +
theme_minimal()
}
# 4. 결정 경계 시각화
plot_decision_regions(x_train_lda, y_train, model = lr_model)
확인해본 결과, 클래스 2로 분류된 1개 데이터만 오분류했다는 점을 알 수 있으며, 규제 강도를 낮추게 되면 보다 정확히 분류할 수 있을 것이다. 위의 모델에 테스트 데이터를 넣고 예측한 결과를 확인해보자.
x_test_lda = lda.transform(x_test_std)
plot_decision_regions(x_test_lda, y_test, classifier=lr)
plt.xlabel("LD 1")
plt.ylabel("LD 2")
plt.legend(loc="lower left")
plt.tight_layout()
plt.show()x_test_lda <- as.matrix(predict(lda_model, newdata = as.data.frame(x_test_std))$x)
plot_decision_regions(x_test_lda, y_test, model = lr_model)
테스트 데이터를 적용해본 결과, 13개의 와인 데이터 특성 대신 2차원의 특성부분공간을 이용해 정확히 분류할 수 있다는 것을 확인할 수 있었다.
728x90
반응형