나이브 베이즈 분류모형

분류모형

• 현실적인 문제로 바꾸어 말하면 어떤 표본에 대한 데이터가 주어졌을 때 그 표본이 어떤 카테고리 혹은 클래스에 속하는지를 알아내는 문제이기도 하다.

• 흔히들 많이 사용하는 모형들의 분류모형 종류를 아래의 표로 기재해 놓았다. 이전에 기재했던 로지스틱 회귀 분석도 실질적으론 분류문제에 많이 사용된다.

모형	방법론
나이브 베이지안	확률적 생성모형
LDA/QDA	확률적 생성모형
로지스틱회귀	확률적 판별모형
의사결정나무	확률적 판별모형
퍼셉트론	판별함수 모형
서포트벡터머신	판별함수 모형
인공신경망	판별함수 모형

나이브 베이즈 분류 모형

• 아래 그림과 같은 질문에 대한 확률을 계산하려면, 아래 결합확률을 계산해야 할 것이다. 독립변수가 2개 즉, Feature의 차원이 2개인 지금은 크게 어렵진 않겠지만, 현실에서는 Feature의 개수가 훨씬 더 많을 것이다. 차원이 커진다면 이에따라 해당 조건부 함수의 가능도함수를 추정하는데 어려움을 겪을 것이다.

$$\begin{align} P( X = weather state, Humidity | Y = playing Tennis) \\ = \frac{P( X_{1} = weather state, Humidity \cap Y = playing Tennis)}{P(Y = playing Tennis)} \end{align}$$

• 위의 수식에서 분자부분의 확률인 결합확률(Joint probability)을 구하는 것은 쉬운 문제가 아니다. 결합되어있는 상황을 나누어서 확률을 계산할 순 없을까라는 가정에서 나이브 베이지안 모형은 출발한다. 그렇다면, 어떻게 나눌수 있을까? 통계를 조금이라도 아시는 분들은 독립이라는 개념을 알고계실 것이다. 예를 들어 A: 코로나가 집단 발병한다, B: 2주간 재택근무를 시행한다. 이 두가지 사건은 서로 발현될 확률에 영향을 미치기 때문에 두 사건은 종속(dependent) 관계라고 할 수 있다. 허나, A: 비가온다, B: 대출심사에서 승인을 받는다 라는 두 사건은 서로 둘 중 어떤 사건에도 발현될 확률에 영향을 주지 않는다. 이러한, 경우에 우리는 두 사건이 확률적으로 독립이라고 정의한다. 두 사건 A,B가 독립일 필요충분조건은 아래와 같다.

$$P(A \cap B) = P(A)P(B)$$

• 이렇게 두 사건이 독립이라면 좀 더 추정하기 쉬울 것이다. 허나, 위의 상황에선 조건부이기 때문에 조건부 독립의 가정이 맞을 것이다. 조건부 독립(conditional independence)은 일반적인 독립과 달리 조건이 되는 별개의 확률변수 C가 존재해야한다. 아래 수식이 성립되면 조건부 독립이라고 한다.

$$A \text{⫫} B \;\vert\; C \Longrightarrow P(A, B | C) = P(A|C)P(B|C)$$

• 조건부독립과 비교하여 일반적인 독립은 무조건부독립이라고 한다. 무조건부 독립은 다음과 같이 표기하기도 한다.

$$A \text{⫫} B \;\vert\; \emptyset$$

• A, B가 C에 대해 조건부 독립이면 다음도 만족한다.

$$P(A \;\vert\; B, C) = P(A \vert C)$$ $$P(B \;\vert\; A, C) = P(B \vert C)$$ $$\begin{align} P(A|B,C) &= \frac{P(A,B,C)}{P(B,C)} \\ &= \frac{P(A,B|C)P(C)}{P(B|C)P(C)} \\ &= \frac{P(A|C)P(B|C)}{P(B|C)} \\ &= P(A|C) \end{align}$$

• 주의할 점은 조건부 독립과 무조건부 독립은 관계가 없다는 점이다. 즉, 두 확률변수가 독립이라고 항상 조건부 독립이 되는 것도 아니고 조건부 독립이라고 꼭 독립이 되는 것도 아니다.

$$P(A,B) = P(A)P(B) \;\; \nRightarrow \;\; P(A,B|C) = P(A|C)P(B|C)$$ $$P(A,B|C) = P(A|C)P(B|C) \;\; \nRightarrow \;\; P(A,B) = P(A)P(B)$$

베이즈 정리(Bayes’ rule)

• 베이즈 정리는 사건 B가 발생함으로써 사건 A의 확률이 어떻게 변화하는지를 표현한 정리이다. 따라서 베이즈 정리는 새로운 정보가 기존의 추론에 어떻게 영향을 미치는지를 나타낸다.

$$P(A|B) = \frac{P(B|A) P(A)}{P(B)}$$

• $P(A|B)$ : 사후확률(Posterior). 사건 B가 발생한 후 갱신된 사건 A의 확률
• $P(A)$ : 사전확률(Prior). 사건 B가 발생하기 전에 가지고 있던 사건 A의 확률
• $P(B|A)$ : 가능도(likelihood). 사건 A가 발생한 경우 사건 B의 확률
• $P(B)$ : 정규화 상수(normalizing constant) 또는 증거(evidence). 확률의 크기 조정역할을 함.

이렇게 독립변수(Feature)들간에 서로 조건부독립임을 가정하여 조건을 나이브하게 하였고, 베이즈정리를 사용하여 MLE를 통해 가장 큰 확률값을 갖는 모수를 추정해내는 모형이 나이브 베이지안 모형이다.

$$\begin{align} P(y = k \mid x) &= \dfrac{ P(x_1, \ldots, x_D \mid y = k) P(y = k) }{P(x)} \\ &= \dfrac{ \left( \prod_{d=1}^D P(x_{d} \mid y = k) \right) P(y = k) }{P(x)} \end{align}$$

• 제일 처음 예시로 들었던 날씨의 상태와 습도를 독립변수로 하고, 종속변수를 테니스를 치는지에 대한 여부에 대한 문제를 푼다고 가정해보자. 먼저 테니스를 치는 경우에 대한 확률을 계산할 경우 테니스를 치는 사건에 대한 확률을 높게 해주는 두가지 경우의 수를 아래와 같이 생각해 볼 수 있다. 아래 두 그림을 수식으로 정리해보면 다음과 같다는 사실을 쉽게 생각해 볼 수 있다. 수식의 마지막에서 분자부분에서 likelihood와 Prior 확률값을 살펴보면 그림에서 언급한 두 가지 경우의 사건들의 확률 값의 곱으로 되어있음을 볼 수 있다. 즉, 사전확률에 새롭게 추가되는 likelihood가 추론의 어떤 영향을 미치는지를 나타내는 베이즈 정리의 의미를 생각해 볼 수 있다.

$$\begin{align} P(weather state, Humidity | playing Tennis) \\ &= \frac{P(weather state)}{P(playing Tennis)} \frac{P(Humidity)}{P(playing Tennis)} \\ &= \frac{P(playing Tennis|weather state)P(weather state)}{P(playing Tennis)} \frac{P(playing Tennis|Humidity)P(Humidity)}{P(playing Tennis)} \end{align}$$

• 아래 그림과 같이 데이터가 주어졌을 때의 베이즈 정리를 이용해 계산하면 다음과 같다.

• 예시로 계산해 본 것과 같이 아래의 그림처럼 독립변수들이 조건부 독립이라는 조건을 통해 곱으로 표현 될 수 있다. 예측값은 마지막 수식에서 보는 것과 같이 해당 가능도함수의 확률값을 최대로 해주는 class를 선택하여 예측값으로 출력해준다.

참조 : MLE 개념

• 이제 하나씩 앞에서 언급했던 순서대로 계산하여 최종적으로 출력을 내는 단계까지 그림을 통해 살펴볼 것이다.

나이브 베이지안 모형의 종류

나이브 가정

• 독립변수 $x$가 $D$차원이라고 가정하자.

$$x = (x_1, \ldots, x_D)$$

• 가능도 함수는 $x_{1}, \ldots, x_{D}$의 결합확률이 된다.

$$P(x \mid y = k) = P(x_1, \ldots, x_D \mid y = k)$$

• 원리상으로는 $y=k$인 데이터만 모아서 이 가능도함수의 모양을 추정할 수 있다. 하지만 차원 $D$가 커지면 가능도함수의 추정이 현실적으로 어려워진다. 따라서 나이브베이즈 분류모형(Naive Bayes classification model)에서는 모든 차원의 개별 독립변수가 서로 조건부 독립(conditional independent)이라는 가정을 사용한다. 이러한 가정을 나이브 가정(naive assumption)이라고 한다.

• 나이브 가정으로 사용하면 벡터 $x$의 결합확률분포함수는 개별 스칼라 원소 $x_{d}$의 확률분포함수의 곱이 된다.

$$P(x_1, \ldots, x_D \mid y = k) = \prod_{d=1}^D P(x_d \mid y = k)$$

• 스칼라 원소 $x_{d}$의 확률분포함수는 결합확률분포함수보다 추정하기 훨씬 쉽다. 가능도함수를 추정한 후에는 베이즈 정리를 사용하며 조건부확률을 계산할 수 있다.

$$\begin{align} P(y = k \mid x) &= \dfrac{ P(x_1, \ldots, x_D \mid y = k) P(y = k) }{P(x)} \\ &= \dfrac{ \left(\prod_{d=1}^D P(x_{d} \mid y = k) \right) P(y = k) }{P(x)} \end{align}$$

가우시안 분포(정규분포) 가능도 모형

• $x$벡터의 원소가 모두 실수이고 클래스마다 특정한 값 주변에서 발생한다고 하면 가능도 분포로 정규분포를 사용한다. 각 독립변수 $x_{d}$마다, 그리고 클래스 $k$마다 정규 분포의 기대값 $\mu_{d,k}$, 표준편차 $\sigma^{2}_{d,k}$가 달라진다. QDA모형과는 달리 모든 독립변수들이 서로 조건부독립이라고 가정한다.

$$P(x_d \mid y = k) = \dfrac{1}{\sqrt{2\pi\sigma_{d,k}^2}} \exp \left(-\dfrac{(x_d-\mu_{d,k})^2}{2\sigma_{d,k}^2}\right)$$

베르누이 분포 가능도 모형

• 베이누이분포 가능도 모형에서는 각각의 $x = (x_1,\ldots, x_D)$의 각 원소 $x_{d}$가 0 또는 1이라는 값만을 가질 수 있다. 독립변수는 $D$개의 독립적인 베르누이 확률변수, 동전으로 구성된 동전 세트로 표현할 수 있다. 이 동전들의 모수 $\mu_{d}$는 동전 $d$마다 다르다.

• 그런데 class $y=k (k = 1,\ldots, K)$마다도 $x_{d}$가 1이 될 확률이 다르다. 즉, 동전의 모수 $\mu_{d,k}$는 동전 $d$마다 다르고 class $k$마다도 다르다. 즉, 전체 $ D \times K $의 조합의 동전이 존재하며 같은 class에 속하는 D개의 동전이 하나의 동전 세트를 구성하고 이러한 동전 세트가 $K$개 있다고 생각할 수 있다.

$$P(x_d \mid y = k) = \mu_{d,k}^{x_d} (1-\mu_{d,k})^{(1-x_d)}$$ $$P(x_1, \ldots, x_D \mid y = k) = \prod_{d=1}^D \mu_{d,k}^{x_d} (1-\mu_{d,k})^{(1-x_d)}$$

• 이러한 동전 세트마다 확률 특성이 다르므로 베르누이분포 가능도 모형을 기반으로 하는 나이브베이즈 모형은 동전 세트를 $N$번 던진 결과로부터 $1, \ldots, K$ 중 어느 동전 세트를 던졌는지를 찾아내는 모형이라고 할 수 있다.

다항분포(Multinomial Distribution) 가능도 모형

• 다항분포 모형에서는 $x$ 벡터가 다항분포의 표본이라고 가정한다. 즉, $D$개의 면을 가지는 주사위를 $\sum_{d=1}^D x_d$ 번 던져서 나온 결과로 본다. 예를 들어 $x$가 다음과 같다면, $x=(1, 4, 0, 5)$ 4면체 주사위를 10번 던져서 1인 면이 1번, 2인 면이 4번, 4인 면이 5번 나온 결과로 해석한다. 각 class 마다 주사위가 다르다고 가정하므로 $K$개의 class를 구분하는 문제에서는 $D$개의 면을 가진 주사위가 $K$개 있다고 본다.

$$P(x_1, \ldots, x_D \mid y = k) \;\; \propto \;\; \prod_{d=1}^D \mu_{d,k}^{x_{d,k}}$$ $$\sum_{d=1}^{D} \mu_{d,k} = 1$$

• 따라서 다항분포 가능도 모형을 기반으로 하는 나이브 베이즈 모형은 주사위를 던진 결과로부터 $1, \ldots, K$ 중 어느 주사위를 던졌는지를 찾아내는 모형이라고 할 수 있다.

Naive Bayes 실습

1. Gaussian Naive Bayes

• 데이터, 모듈 불러오기
  • iris(붓꽃)데이터를 가지고 실습을 진행할 것이다.
  • 붓꽃의 종류는 3가지(‘setosa’, ‘versicolor’, ‘virginica’)이며, 각 feature는 sepal length (cm), sepal width (cm), petal length (cm), petal width (cm) 로 이루어져있다. 모든 피처가 실수의 값을 갖기 때문에 가우시안 나이브 베이즈 모형을 사용할 것이다.

import pandas as pd
import numpy as np
from sklearn import datasets
from sklearn.naive_bayes import GaussianNB
from sklearn.model_selection import train_test_split

iris = datasets.load_iris()
df_X=pd.DataFrame(iris.data, columns=iris.feature_names)
df_Y=pd.DataFrame(iris.target, columns=["target"])

• 예측의 성능은 항상 테스트 데이터를 통해 측정해야 하므로 분리해준다. 지금은 연습삼아 해보는 것이므로 필자는 따로 validation set을 구축하지 않겠다. 허나 실제로 데이터를 분석하고 그에 따른 성능을 측정하려면 꼭 validation set을 따로 두어 베이스라인 모델의 parameter들을 튜닝을 통해 모델의 성능을 높이도록 한 뒤 최종적으로 test set을 예측하여 성능을 검증해야 할 것이다.

train_X, test_X, train_Y, test_Y = train_test_split(df_X, df_Y, train_size=0.8, test_size=0.2, random_state=123)

• 모델 피팅

GNB = GaussianNB()
fitted = GNB.fit(train_X, train_Y)
y_pred = fitted.predict(test_X)
y_pred

• 확률 구하기
  • 먼저 예측한 클래스와 해당 예측 데이터의 클래스별 확률을 살펴 볼 것이다.

fitted.predict(test_X)[:1]

결과

1

fitted.predict_proba(test_X)[:1]

결과

array([[7.24143720e-126, 9.23061979e-001, 7.69380215e-002]])

• 위의 확률 직접 구해보기

1) 위의 확률 값이 나오게 된 중간과정을 살펴보자. 우선 추정된 독립변수의 모수와 정규 분포의 확률밀도 함수를 사용하여 가능도를 구할 수 있다.

$$p(x_1, x_2, x_3, x_4 | y) \propto p(x_1) p(x_2) p(x_3) p(x_4)$$

- 위에서 학습했을 때도 말했던 것 처럼, (class개수 * 변수의 개수)개의 조합의 모수를 갖고 있으므로 아래와 같이 12개의 모수를 갖는다. 아래에서는 각 class별로 정규분포의 모수인 평균과 분산을 보여준다.

predict_data = np.array(test_X.iloc[0])
predict_data

결과

array([6.3, 2.5, 4.9, 1.5])

• 추정한 모델의 클래스별 모수(평균과 분산)을 다음과 같이 알 수 있다.

fitted.theta_[0], fitted.sigma_[0]

결과

(array([5.01621622, 3.43243243, 1.46756757, 0.25945946]),
 array([0.10568298, 0.14975895, 0.02705625, 0.01214025]))

fitted.theta_[1], fitted.sigma_[1]

결과

(array([5.95      , 2.78409091, 4.24090909, 1.32272727]),
 array([0.27068182, 0.10042872, 0.22741736, 0.04221075]))

fitted.theta_[2], fitted.sigma_[2]

결과

(array([6.58717949, 2.95897436, 5.57948718, 2.02820513]),
 array([0.39752795, 0.11011177, 0.29188692, 0.0774096 ]))

• 위의 모수들을 통해 class별 가능도를 구하면 아래와 같을 것이다.

likelihood = [
(sp.stats.norm(fitted.theta_[0][0], np.sqrt(fitted.sigma_[0][0])).pdf(predict_data[0]) * \
sp.stats.norm(fitted.theta_[0][1], np.sqrt(fitted.sigma_[0][1])).pdf(predict_data[1]) * \
sp.stats.norm(fitted.theta_[0][2], np.sqrt(fitted.sigma_[0][2])).pdf(predict_data[2]) * \
sp.stats.norm(fitted.theta_[0][3], np.sqrt(fitted.sigma_[0][3])).pdf(predict_data[3])),\
(sp.stats.norm(fitted.theta_[1][0], np.sqrt(fitted.sigma_[1][0])).pdf(predict_data[0]) * \
sp.stats.norm(fitted.theta_[1][1], np.sqrt(fitted.sigma_[1][1])).pdf(predict_data[1]) * \
sp.stats.norm(fitted.theta_[1][2], np.sqrt(fitted.sigma_[1][2])).pdf(predict_data[2]) * \
sp.stats.norm(fitted.theta_[1][3], np.sqrt(fitted.sigma_[1][3])).pdf(predict_data[3])),\
(sp.stats.norm(fitted.theta_[2][0], np.sqrt(fitted.sigma_[0][0])).pdf(predict_data[0]) * \
sp.stats.norm(fitted.theta_[2][1], np.sqrt(fitted.sigma_[0][1])).pdf(predict_data[1]) * \
sp.stats.norm(fitted.theta_[2][2], np.sqrt(fitted.sigma_[0][2])).pdf(predict_data[2]) * \
sp.stats.norm(fitted.theta_[2][3], np.sqrt(fitted.sigma_[0][3])).pdf(predict_data[3]))    
]
likelihood

결과

[2.0700298536453225e-126, 0.2218869448618605, 7.497361843154609e-09]

• 여기에 사전확률을 곱하면 사후 확률에 비례하는 값이 나온다.

$$p(y|x_1, x_2) \propto p(x_1, x_2|y) p(y)$$

아직 정규화 상수 $p(x)$로 나누어주지 않았으므로 두 값의 합이 1이 아니다. 즉, 확률이라고 부를 수는 없다. 하지만 크기를 비교하면 이 데이터는 $y=1$일 확률이 가장 높다는 것을 알 수 있다.

fitted.class_prior_

결과

array([0.30833333, 0.36666667, 0.325     ])

posterior = likelihood * fitted.class_prior_
posterior

결과

array([6.38259205e-127, 8.13585464e-002, 2.43664260e-009])

• 이 값을 정규화하면 predict_proba 메서드로 구한 것과 같은 값이 나온다. 물론 완벽히 일치하진 않지만 그에 근사하는 값을 추정값으로 계산해내서 얻을 수 있다. 이는 계산시 소수점을 어느정도까지 사용하였는지에 따라 다르기 때문에 이 정도의 오차는 문제가 없다.

posterior / np.sum(posterior, axis=0)

결과

array([7.84501707e-126, 9.99999970e-001, 2.99494353e-008])

fitted.predict_proba(test_X)[:1]

결과

array([[7.24143720e-126, 9.23061979e-001, 7.69380215e-002]])

• Confusion matrix 구하기

from sklearn.metrics import confusion_matrix
confusion_matrix(y_pred, test_Y)

결과

array([[13,  0,  0],
       [ 0,  6,  1],
       [ 0,  0, 10]])

• Prior 설정하기
  • 이번에는 class가 발생되는 사전확률을 미리 알고 있었던 경우라고 가정하고 문제를 풀어볼 것이다.

GNB2 = GaussianNB(priors=[0.01, 0.01, 0.98])
set_prior_fitted_01 = GNB2.fit(train_X, train_Y)
set_prior_pred_01 = set_prior_fitted_01.predict(test_X)
confusion_matrix(set_prior_pred_01, test_Y)

결과

array([[13,  0,  0],
       [ 0,  4,  0],
       [ 0,  2, 11]])

GNB3 = GaussianNB(priors=[0.01, 0.98, 0.01])
set_prior_fitted_02 = GNB3.fit(train_X, train_Y)
set_prior_pred_02 = set_prior_fitted_02.predict(test_X)
confusion_matrix(set_prior_pred_02, test_Y)

array([[13,  0,  0],
       [ 0,  6,  4],
       [ 0,  0,  7]])

2. Bernoulli naive bayes

• e-mail과 같은 문서 내에 특정한 단어가 포함되어 있는지의 여부는 베르누이 확률변수로 모형화할 수 있다. 이렇게 독립변수가 0 또는 1의 값을 가지면 베르누이 나이브베이즈 모형을 사용한다.

• python의 sklearn의 베르누이분포 나이브베이즈 모형 클래스 BernoulliNB는 가능도 추정과 관련하여 다음 속성을 가진다.

  • feature_count_ : 각 class k에 대해 d번째 동전이 앞면이 나온 횟수 $N_{d,k}$
  • feature_log_prob_ : 베르누이분포 모수의 로그값

$$\log \mu_k = (\log \mu_{1,k}, \ldots, \log \mu_{D, k}) = \left( \log \dfrac{N_{1,k}}{N_k}, \ldots, \log \dfrac{N_{D,k}}{N_k} \right)$$

• 여기에서 $N_{k}$은 클래스 k에 대해 동전을 던진 총 횟수이다. 표본 데이터의 수가 적은 경우에는 모수에 대해 다음처럼 스무딩(smoothing)을 할 수도 있다.

스무딩(Smoothing)

• 표본 데이터의 수가 적은 경우에는 베르누이 모수가 0 또는 1이라는 극단적인 모수 추정값이 나올 수도 있다. 하지만 현실적으로는 실제 모수값이 이런 극단적인 값이 나올 가능성이 적다. 따라서 베르누이 모수가 0.5인 가장 일반적인 경우를 가정하여 0이 나오는 경우와 1이 나오는 경우, 두 개의 가상 표본 데이터를 추가한다. 그러면 0이나 1과 같은 극단적인 추정값이 0.5에 가까운 다음과 같은 값으로 변한다. 이를 라플라스 스무딩(Laplace smoothing) 또는 애드원(Add-One) 스무딩이라고 한다.

$$\hat{\mu}_{d,k} = \frac{ N_{d,k} + \alpha}{N_k + 2 \alpha}$$

• 가중치 $\alpha$를 사용하여 스무딩의 정도를 조절할 수도 있다. 가중치 $\alpha$는 정수가 아니라도 괜찮다. 가중치가 1인 경우는 무정보 사전확률을 사용한 베이즈 모수추정의 결과와 같다.

• 아래의 데이터는 4개의 key word를 사용하여 정상 메일 4개와 spam 메일 6개를 BOW 인코딩한 행렬이다. 예를 들어 첫번째 메일은 정상 메일이고 1번, 4번 key word는 포함하지 않지만 2번,3번 key word를 포함한다고 볼 수 있다.

from sklearn.naive_bayes import BernoulliNB
X = np.array([
    [0, 1, 1, 0],
    [1, 1, 1, 1],
    [1, 1, 1, 0],
    [0, 1, 0, 0],
    [0, 0, 0, 1],
    [0, 1, 1, 0],
    [0, 1, 1, 1],
    [1, 0, 1, 0],
    [1, 0, 1, 1],
    [0, 1, 1, 0]])
y = np.array([0, 0, 0, 0, 1, 1, 1, 1, 1, 1])

BNB = BernoulliNB()
fitted = BNB.fit(X, y)

• y 클래스의 종류와 각 클래스에 속하는 표본의 수, 그리고 그 값으로부터 구한 사전 확률의 값은 다음과 같다.

fitted.classes_

결과

array([0, 1])

fitted.class_count_

결과

array([4., 6.])

np.exp(fitted.class_log_prior_)

결과

array([0.4, 0.6])

• 각 클래스 k 별로, 그리고 각 독립변수 d 별로, 각각 다른 베르누이 확률변수라고 가정하여 모두 8개의 베르누이 확률변수의 모수를 구하면 다음과 같다.

# 행은 클래스의 개수, 열은 변수의 개수를 의미
count = fitted.feature_count_
count

결과

array([[2., 4., 3., 1.],
       [2., 3., 5., 3.]])

• 위의 각 변수의 클래스별로 몇번 나온지에 대한 행렬에 각 클래스의 전체 개수로 나누어 주어야 해당 변수들의 모수인 p값을 알 수 있으므로 class_count의 배열의 모양을 변형시켜주어야 한다. 현재는 1차원의 벡터(2,)이므로 2차원의 (2,1)의 모양을 갖도록 해주어야 나누어 줄 수 있기 때문이다.

count / fitted.class_count_[:,np.newaxis]

결과

array([[0.5       , 1.        , 0.75      , 0.25      ],
       [0.33333333, 0.5       , 0.83333333, 0.5       ]])

• 위에서는 필자는 Numpy의 brodcasting 연산을 사용하여 구한것인데, 혹시 count 행렬의 모양과 동일하게 만들어 확실하게 연산하고 싶다면, 아래와 같이 실행하면 동일한 결과를 얻을 수 있는 것을 볼 수 있다.

count / np.repeat(fitted.class_count_[:,np.newaxis], 4, axis=1)

결과

array([[0.5       , 1.        , 0.75      , 0.25      ],
       [0.33333333, 0.5       , 0.83333333, 0.5       ]])

• 그런데 이 값은 모형 내에서 구한 값과 다르다. 모형 내에서 스무딩(smoothing)이 이루어지기 때문이다. 스무딩은 동전의 각 면 즉, 0과 1이 나오는 가상의 데이터를 추가함으로서 추정한 모수의 값이 좀 더 0.5에 가까워지도록 하는 방법이다. 이 때 사용한 스무딩 가중치 값은 다음처럼 확인할 수 있다.

fitted.alpha

결과

1.0

theta = np.exp(fitted.feature_log_prob_)
theta

결과

array([[0.5       , 0.83333333, 0.66666667, 0.33333333],
       [0.375     , 0.5       , 0.75      , 0.5       ]])

• 이에 모형이 완성되었으니 테스트 데이터를 사용하여 예측을 해 본다. 예를 들어 1번, 2번 키워드를 포함한 메일이 정상 메일인지 스팸 메일인지 알아보자.

x_new = np.array([0, 1, 1, 1])
fitted.predict_proba([x_new])

결과

array([[0.34501348, 0.65498652]])

• 위 결과에서 정상 메일일 가능성이 약 3배임을 알 수 있다. 이 값은 다음처럼 구할 수도 있다.

# np.prod(axis=1)로 해준 이유는 아래 각 변수별로 독립이므로 가능도를 구하려면 곱을 해주어야 하기 때문
p = ((theta ** x_new) * (1 - theta) ** (1 - x_new)).prod(axis=1) * np.exp(fitted.class_log_prior_)
p / p.sum()

결과

array([0.34501348, 0.65498652])

• 반대로 3번, 4번 keyword가 포함된 메일은 스팸일 가능성이 약 90%이다.

x_new = np.array([0, 0, 1, 1])
fitted.predict_proba([x_new])

결과

array([[0.09530901, 0.90469099]])

p = ((theta ** x_new) * (1 - theta) ** (1 - x_new)).prod(axis=1) * np.exp(fitted.class_log_prior_)
p / p.sum()

결과

array([[0.09530901, 0.90469099]])

• MNIST 숫자 분류문제에서 sklearn.preprocessing.Binarizer로 x값을 0, 1로 바꾼다(값이 8 이상이면 1, 8 미만이면 0). 즉 흰색과 검은색 픽셀로만 구성된 이미지로 만든다(다음 코드 참조)

from sklearn.datasets import load_digits
digits = load_digits()
X = digits.data
y = digits.target
from sklearn.preprocessing import Binarizer
X = Binarizer(7).fit_transform(X)

• 이 이미지에 대해 베르누이 나이브베이즈 모형을 적용하자. 분류 결과를 분류보고서 형식으로 나타내라.

  • (1) BernoulliNB 클래스의 binarize 인수를 사용하여 같은 문제를 풀어본다.

  • (2) 계산된 모형의 모수 벡터 값을 각 클래스별로 8x8 이미지의 형태로 나타내라. 이 이미지는 무엇을 뜻하는가?

digits.images[0]

결과

array([[ 0.,  0.,  5., 13.,  9.,  1.,  0.,  0.],
       [ 0.,  0., 13., 15., 10., 15.,  5.,  0.],
       [ 0.,  3., 15.,  2.,  0., 11.,  8.,  0.],
       [ 0.,  4., 12.,  0.,  0.,  8.,  8.,  0.],
       [ 0.,  5.,  8.,  0.,  0.,  9.,  8.,  0.],
       [ 0.,  4., 11.,  0.,  1., 12.,  7.,  0.],
       [ 0.,  2., 14.,  5., 10., 12.,  0.,  0.],
       [ 0.,  0.,  6., 13., 10.,  0.,  0.,  0.]])

X.shape

결과

(1797, 64)

BNB =  BernoulliNB()
fitted = BNB.fit(X, y)
theta = np.exp(fitted.feature_log_prob_)
theta = theta.reshape((10, 8, 8))

fig, axes = plt.subplots(2, 5, figsize=(12, 3),
                         subplot_kw={'xticks': [], 'yticks': []})
for i in range(5):
    axes[0][i].set_title("class {}".format(i))
    axes[0][i].imshow(theta[i], interpolation='nearest', cmap=plt.cm.Blues)
    axes[1][i].set_title("class {}".format(i+5))
    axes[1][i].imshow(theta[i+5], interpolation='nearest', cmap=plt.cm.Blues)
plt.show()

• 위의 이미지에서는 모수값이 높은 변수가 진한 파란색을 띄게 된다. sklearn.preprocessing.Binarizer를 통해 x값을 값이 8 이상이면 1, 8 미만이면 0으로 바꾸어주었으므로 8 미만인 데이터보다는 8이상인 데이터가 각 클래스를 구분하는데 좀 더 영향을 주는 공간을 알 수 있게 해준다.

3. Multinomial naive bayes

• 다항분포 나이브베이즈 모형 클래스 MultinomialNB는 가능도 추정과 관련하여 다음 속성을 가진다.
  • feature_count_ : 각 클래스 $k$에서 $d$번째 면이 나온 횟수 $N_{d,k}$
  • feature_log_prob_ : 다항분포의 모수의 로그

$$\log \mu_k = (\log \mu_{1,k}, \ldots, \log \mu_{D, k}) = \left( \log \dfrac{N_{1,k}}{N_k}, \ldots, \log \dfrac{N_{D,k}}{N_k} \right)$$

• 여기에서 $N_{k}$은 클래스 $k$에 대해 주사위를 던진 총 회수를 뜻한다.

• 스무딩공식은 아래와 같다.

$$\hat{\mu}_{d,k} = \frac{ N_{d,k} + \alpha}{N_k + D \alpha} , \; (D=변수의 개수)$$

• 이번에도 스팸 메일 필터링을 예로 들어보다. 다만 BOW 인코딩을 할 때, 각 키워드가 출현한 빈도를 직접 입력 변수로 사용한다.

X = np.array([
    [3, 4, 1, 2],
    [3, 5, 1, 1],
    [3, 3, 0, 4],
    [3, 4, 1, 2],
    [1, 2, 1, 4],
    [0, 0, 5, 3],
    [1, 2, 4, 1],
    [1, 1, 4, 2],
    [0, 1, 2, 5],
    [2, 1, 2, 3]])
y = np.array([0, 0, 0, 0, 1, 1, 1, 1, 1, 1])

from sklearn.naive_bayes import MultinomialNB
MNB = MultinomialNB()
fitted = MNB.fit(X, y)

fitted.classes_

결과

array([0, 1])

fitted.class_count_

결과

array([4., 6.])

np.exp(fitted.class_log_prior_)

결과

array([0.4, 0.6])

• 다음으로 각 클래스에 대한 가능도 확률분포를 구한다. 다항분포 모형을 사용하므로 각 클래스를 4개의 면을 가진 주사위로 생각할 수 있다. 그리고 각 면이 나올 확률은 각 면이 나온 횟수를 주사위를 던진 전체 횟수로 나누면 된다. 우선 각 클래스 별로 각각의 면이 나온 횟수는 다음과 같다.

count = fitted.feature_count_
count

결과

array([[12., 16.,  3.,  9.],
       [ 5.,  7., 18., 18.]])

• 이 데이터에서 클래스 $Y=0$인 주사위를 던진 횟수는 첫번째 행의 값의 합인 40이므로 클래스 $Y=0$인 주사위를 던져 1이라는 면이 나올 확률은 다음처럼 계산할 수 있다.

$$\mu_{1,Y=0} = \dfrac{12}{40} = 0.3$$

count / np.repeat(count.sum(axis=1)[:,np.newaxis], 4, axis=1)

결과

array([[0.3       , 0.4       , 0.075     , 0.225     ],
       [0.10416667, 0.14583333, 0.375     , 0.375     ]])

• 실제로는 극단적인 추정을 피하기 위해 이 값을 가중치 1인 스무딩을 한 추정값을 사용한다.

fitted.alpha

결과

1.0

(count + fitted.alpha) / \
(np.repeat(count.sum(axis=1)[:,np.newaxis], 4, axis=1) + fitted.alpha * X.shape[1])

결과

array([[0.29545455, 0.38636364, 0.09090909, 0.22727273],
       [0.11538462, 0.15384615, 0.36538462, 0.36538462]])

theta = np.exp(fitted.feature_log_prob_)
theta

결과

array([[0.29545455, 0.38636364, 0.09090909, 0.22727273],
       [0.11538462, 0.15384615, 0.36538462, 0.36538462]])

• 이제 이 값을 사용하여 예측을 해 보자. 만약 어떤 메일에 1번부터 4번까지의 키워드가 각각 10번씩 나왔다면 다음처럼 확률을 구할 수 있다. 구해진 확률로부터 이 메일이 스팸임을 알 수 있다.

x_new = np.array([10, 10, 10, 10])
fitted.predict_proba([x_new])

결과

array([[0.38848858, 0.61151142]])

p = (theta ** x_new).prod(axis=1) * np.exp(fitted.class_log_prior_)
p / p.sum(axis=0)

결과

array([0.38848858, 0.61151142])

• MNIST 숫자 분류문제를 다항분포 나이브베이즈 모형을 사용하여 풀고 이진화(Binarizing)를 하여 베르누이 나이브베이즈 모형을 적용했을 경우와 성능을 비교하라.

from sklearn.datasets import load_digits
from sklearn.metrics import classification_report
digits = load_digits()
X = digits.data
y = digits.target
train_X, test_X, train_Y, test_Y = train_test_split(X, y, test_size=0.3, random_state=123)
from sklearn.preprocessing import Binarizer
binary_train_X = Binarizer(7).fit_transform(train_X)
binary_test_X = Binarizer(7).fit_transform(test_X)
BNB = BernoulliNB().fit(binary_train_X, train_Y)
bnb_pred = BNB.predict(binary_test_X)
MNB = MultinomialNB().fit(train_X, train_Y)
mnb_pred = MNB.predict(test_X)

이진화 한 베르누이 나이브베이즈 모형 성능

print(classification_report(bnb_pred, test_Y))

결과

            precision    recall  f1-score   support

           0       0.98      0.98      0.98        59
           1       0.82      0.84      0.83        55
           2       0.87      0.90      0.88        51
           3       0.80      0.93      0.86        40
           4       0.95      0.97      0.96        60
           5       0.84      0.94      0.89        51
           6       0.96      1.00      0.98        55
           7       1.00      0.89      0.94        56
           8       0.85      0.80      0.83        51
           9       0.87      0.74      0.80        62

    accuracy                           0.90       540
   macro avg       0.90      0.90      0.90       540
weighted avg       0.90      0.90      0.90       540

다항 분포 나이브 베이즈 모형 성능

print(classification_report(mnb_pred, test_Y))

결과

              precision    recall  f1-score   support

           0       0.98      1.00      0.99        58
           1       0.77      0.88      0.82        49
           2       0.83      0.88      0.85        50
           3       0.89      1.00      0.94        41
           4       0.97      0.92      0.94        64
           5       0.74      1.00      0.85        42
           6       0.98      0.98      0.98        57
           7       1.00      0.89      0.94        56
           8       0.81      0.74      0.77        53
           9       0.89      0.67      0.76        70

    accuracy                           0.89       540
   macro avg       0.89      0.90      0.89       540
weighted avg       0.89      0.89      0.89       540

• 텍스트 분석에서 TF-IDF 인코딩을 하면 단어의 빈도수가 정수가 아닌 실수값이 된다. 이런 경우에도 다항분포 모형을 적용할 수 있는가?

  • 정수가 아니더라도 해당 적용 가능하다! 다항분포의 모수를 추정할 때 해당 관측치의 변수가 갖는 값의 합으로 나누어주어 모수값을 구했는데, TF-idf 행렬은 row가 문서를 의미하고 전체 문서에서의 토큰들이 열을 이루게 되므로 해당 문서에서 어떠한 단어가 몇번 나온것인지에 대해 다항분포를 통해 계산할 수 있기 때문이다.

• MultinomialNB를 사용할 경우 범주형으로 정수이면 사용하는 것이라고 생각하지말고 위의 예시 처럼 데이터 당 각 피처가 유기적으로 하나의 사건에서 파생되어 이루어질 수 있는지에 대해서 먼저 생각해보자. 통계적인 분포를 다항분포로 생각할 수 있는지를 확인해보자는 이야기이다.

• 아래의 뉴스그룹 분류 문제를 통해 검증해보자.

뉴스그룹 분류

• 다음은 뉴스그룹 데이터에 대해 나이브베이즈 분류모형을 적용한 결과이다.
  • 문서는18846건

from sklearn.datasets import fetch_20newsgroups

news = fetch_20newsgroups(subset="all")
X = news.data
y = news.target

from sklearn.feature_extraction.text import TfidfVectorizer, HashingVectorizer, CountVectorizer
from sklearn.naive_bayes import MultinomialNB
from sklearn.pipeline import Pipeline

model1 = Pipeline([
    ('vect', CountVectorizer()),
    ('model', MultinomialNB()),
])
model2 = Pipeline([
    ('vect', TfidfVectorizer()),
    ('model', MultinomialNB()),
])
model3 = Pipeline([
    ('vect', TfidfVectorizer(stop_words="english")),
    ('model', MultinomialNB()),
])
model4 = Pipeline([
    ('vect', TfidfVectorizer(stop_words="english",
                             token_pattern=r"\b[a-z0-9_\-\.]+[a-z][a-z0-9_\-\.]+\b")),
    ('model', MultinomialNB()),
])

%%time
from sklearn.model_selection import cross_val_score, KFold

for i, model in enumerate([model1, model2, model3, model4]):
    scores = cross_val_score(model, X, y, cv=5)
    print(("Model{0:d}: Mean score: {1:.3f}").format(i + 1, np.mean(scores)))

결과

Model1: Mean score: 0.855
Model2: Mean score: 0.856
Model3: Mean score: 0.883
Model4: Mean score: 0.888
CPU times: user 1min 35s, sys: 4.54 s, total: 1min 40s
Wall time: 1min 53s

• (1) 만약 독립변수로 실수 변수, 0 또는 1 값을 가지는 변수, 자연수 값을 가지는 변수가 섞여있다면 사이킷런에서 제공하는 나이브베이즈 클래스를 사용하여 풀 수 있는가?

- 위에서 likelihood를 직접 계산했던 것과 같이 `likelihood만을 각각 계산하여 각 변수들은 독립이라는 가정을 전제로하기 때문에 서로 곱한뒤에 베이즈정리식에 따라 최종적으로 확률값을 구해 클래스를 구분`할 수 있다.

• (2) 사이킷런에서 제공하는 분류문제 예제 중 숲의 수종을 예측하는 covtype 분류문제는 연속확률분포 특징과 베르누이확률분포 특징이 섞여있다. 이 문제를 사이킷런에서 제공하는 나이브베이즈 클래스를 사용하여 풀어라.

대표 수종 데이터(covtype)

• 대표 수종 데이터는 미국 삼림을 30×30m 영역으로 나누어 각 영역의 특징으로부터 대표적인 나무의 종류(species of tree)을 예측하기위한 데이터이다. 수종은 7종류이지만 특징 데이터가 54종류, 표본 데이터의 갯수가 581,012개에 달하는 대규모 데이터이다.

from sklearn.datasets import fetch_covtype
covtype = fetch_covtype()
# print(covtype.DESCR)

df = pd.DataFrame(covtype.data,
                  columns=["x{:02d}".format(i + 1) for i in range(covtype.data.shape[1])],
                  dtype=int)
sy = pd.Series(covtype.target, dtype="category")
df['covtype'] = sy

• 각 특징 데이터가 가지는 값의 종류를 보면 1번부터 10번 특징은 실수값이고 11번부터 54번 특징은 이진 카테고리값이라는 것을 알 수 있다.

df.iloc[:, 10:54] = df.iloc[:, 10:54].astype('category')

• 다음 플롯은 카테고리값에 따라 “x14” 특징의 값이 어떻게 변하는지 나타낸 것이다. “x14” 특징이 0인가 1인가를 사용하면 1, 5, 7번 클래스와 4번 클래스는 완벽하게 분류할 수 있다는 것을 알 수 있다.

plt.figure(figsize=(10,12))
df_count = df.pivot_table(index="covtype", columns="x14", aggfunc="size")
sns.heatmap(df_count, cmap=sns.light_palette("gray", as_cmap=True), annot=True, fmt="0")
plt.show()

Gaussian_df_X = df.iloc[:, :10]
Bern_df_X = df.iloc[:, 10:-1]
df_Y = df.iloc[:, -1]

GNB = GaussianNB()
fitted_GNB = GNB.fit(Gaussian_df_X, df_Y)

theta = fitted_GNB.theta_
sigma = fitted_GNB.sigma_

• 가능도를 계산하기 위한 함수를 작성하였다.
  • 아래 함수는 반복문을 통해 실행하는 방식인데 데이터 수가 많다면 너무 비효율적이다. 그러므로, 데이터(관측치)의 수가 적은 경우에만 이용하는 것을 권장한다.
  • 그래서 또 직접적으로 for문을 돌리지 않고 사용할 수 있는 방식의 함수를 다시 구현하였다. for문으로 반복문을 작성한것과 비교했을 때는 직관적으로 어떻게 가능도를 구하는 지 알 수 있다.

def Gaussian_likelihood_cal(predict_data, theta, sigma, class_count, feature_count):
    likelihood = []
    for c in np.arange(class_count):
        prod = 1
        for f in np.arange(feature_count):
            prod = prod * sp.stats.norm(theta[c][f], np.sqrt(sigma[c][f])).pdf(predict_data[f])
        likelihood.append(prod)
    return likelihood

def Gaussian_likelihood_cal(predict_data, theta, sigma):
    likelihood = [(sp.stats.norm(theta[0][0], np.sqrt(sigma[0][0])).pdf(predict_data[0]) *\
                   sp.stats.norm(theta[0][1], np.sqrt(sigma[0][1])).pdf(predict_data[1]) *\
                   sp.stats.norm(theta[0][2], np.sqrt(sigma[0][2])).pdf(predict_data[2]) *\
                   sp.stats.norm(theta[0][3], np.sqrt(sigma[0][3])).pdf(predict_data[3]) *\
                   sp.stats.norm(theta[0][4], np.sqrt(sigma[0][4])).pdf(predict_data[4]) *\
                   sp.stats.norm(theta[0][5], np.sqrt(sigma[0][5])).pdf(predict_data[5]) *\
                   sp.stats.norm(theta[0][6], np.sqrt(sigma[0][6])).pdf(predict_data[6]) *\
                   sp.stats.norm(theta[0][7], np.sqrt(sigma[0][7])).pdf(predict_data[7]) *\
                   sp.stats.norm(theta[0][8], np.sqrt(sigma[0][8])).pdf(predict_data[8]) *\
                   sp.stats.norm(theta[0][9], np.sqrt(sigma[0][9])).pdf(predict_data[9])),\
                  (sp.stats.norm(theta[1][0], np.sqrt(sigma[1][0])).pdf(predict_data[0]) *\
                   sp.stats.norm(theta[1][1], np.sqrt(sigma[1][1])).pdf(predict_data[1]) *\
                   sp.stats.norm(theta[1][2], np.sqrt(sigma[1][2])).pdf(predict_data[2]) *\
                   sp.stats.norm(theta[1][3], np.sqrt(sigma[1][3])).pdf(predict_data[3]) *\
                   sp.stats.norm(theta[1][4], np.sqrt(sigma[1][4])).pdf(predict_data[4]) *\
                   sp.stats.norm(theta[1][5], np.sqrt(sigma[1][5])).pdf(predict_data[5]) *\
                   sp.stats.norm(theta[1][6], np.sqrt(sigma[1][6])).pdf(predict_data[6]) *\
                   sp.stats.norm(theta[1][7], np.sqrt(sigma[1][7])).pdf(predict_data[7]) *\
                   sp.stats.norm(theta[1][8], np.sqrt(sigma[1][8])).pdf(predict_data[8]) *\
                   sp.stats.norm(theta[1][9], np.sqrt(sigma[1][9])).pdf(predict_data[9])),\
                  (sp.stats.norm(theta[2][0], np.sqrt(sigma[2][0])).pdf(predict_data[0]) *\
                   sp.stats.norm(theta[2][1], np.sqrt(sigma[2][1])).pdf(predict_data[1]) *\
                   sp.stats.norm(theta[2][2], np.sqrt(sigma[2][2])).pdf(predict_data[2]) *\
                   sp.stats.norm(theta[2][3], np.sqrt(sigma[2][3])).pdf(predict_data[3]) *\
                   sp.stats.norm(theta[2][4], np.sqrt(sigma[2][4])).pdf(predict_data[4]) *\
                   sp.stats.norm(theta[2][5], np.sqrt(sigma[2][5])).pdf(predict_data[5]) *\
                   sp.stats.norm(theta[2][6], np.sqrt(sigma[2][6])).pdf(predict_data[6]) *\
                   sp.stats.norm(theta[2][7], np.sqrt(sigma[2][7])).pdf(predict_data[7]) *\
                   sp.stats.norm(theta[2][8], np.sqrt(sigma[2][8])).pdf(predict_data[8]) *\
                   sp.stats.norm(theta[2][9], np.sqrt(sigma[2][9])).pdf(predict_data[9])),\
                  (sp.stats.norm(theta[3][0], np.sqrt(sigma[3][0])).pdf(predict_data[0]) *\
                   sp.stats.norm(theta[3][1], np.sqrt(sigma[3][1])).pdf(predict_data[1]) *\
                   sp.stats.norm(theta[3][2], np.sqrt(sigma[3][2])).pdf(predict_data[2]) *\
                   sp.stats.norm(theta[3][3], np.sqrt(sigma[3][3])).pdf(predict_data[3]) *\
                   sp.stats.norm(theta[3][4], np.sqrt(sigma[3][4])).pdf(predict_data[4]) *\
                   sp.stats.norm(theta[3][5], np.sqrt(sigma[3][5])).pdf(predict_data[5]) *\
                   sp.stats.norm(theta[3][6], np.sqrt(sigma[3][6])).pdf(predict_data[6]) *\
                   sp.stats.norm(theta[3][7], np.sqrt(sigma[3][7])).pdf(predict_data[7]) *\
                   sp.stats.norm(theta[3][8], np.sqrt(sigma[3][8])).pdf(predict_data[8]) *\
                   sp.stats.norm(theta[3][9], np.sqrt(sigma[3][9])).pdf(predict_data[9])),\
                  (sp.stats.norm(theta[4][0], np.sqrt(sigma[4][0])).pdf(predict_data[0]) *\
                   sp.stats.norm(theta[4][1], np.sqrt(sigma[4][1])).pdf(predict_data[1]) *\
                   sp.stats.norm(theta[4][2], np.sqrt(sigma[4][2])).pdf(predict_data[2]) *\
                   sp.stats.norm(theta[4][3], np.sqrt(sigma[4][3])).pdf(predict_data[3]) *\
                   sp.stats.norm(theta[4][4], np.sqrt(sigma[4][4])).pdf(predict_data[4]) *\
                   sp.stats.norm(theta[4][5], np.sqrt(sigma[4][5])).pdf(predict_data[5]) *\
                   sp.stats.norm(theta[4][6], np.sqrt(sigma[4][6])).pdf(predict_data[6]) *\
                   sp.stats.norm(theta[4][7], np.sqrt(sigma[4][7])).pdf(predict_data[7]) *\
                   sp.stats.norm(theta[4][8], np.sqrt(sigma[4][8])).pdf(predict_data[8]) *\
                   sp.stats.norm(theta[4][9], np.sqrt(sigma[4][9])).pdf(predict_data[9])),\
                  (sp.stats.norm(theta[5][0], np.sqrt(sigma[5][0])).pdf(predict_data[0]) *\
                   sp.stats.norm(theta[5][1], np.sqrt(sigma[5][1])).pdf(predict_data[1]) *\
                   sp.stats.norm(theta[5][2], np.sqrt(sigma[5][2])).pdf(predict_data[2]) *\
                   sp.stats.norm(theta[5][3], np.sqrt(sigma[5][3])).pdf(predict_data[3]) *\
                   sp.stats.norm(theta[5][4], np.sqrt(sigma[5][4])).pdf(predict_data[4]) *\
                   sp.stats.norm(theta[5][5], np.sqrt(sigma[5][5])).pdf(predict_data[5]) *\
                   sp.stats.norm(theta[5][6], np.sqrt(sigma[5][6])).pdf(predict_data[6]) *\
                   sp.stats.norm(theta[5][7], np.sqrt(sigma[5][7])).pdf(predict_data[7]) *\
                   sp.stats.norm(theta[5][8], np.sqrt(sigma[5][8])).pdf(predict_data[8]) *\
                   sp.stats.norm(theta[5][9], np.sqrt(sigma[5][9])).pdf(predict_data[9])),\
                  (sp.stats.norm(theta[6][0], np.sqrt(sigma[6][0])).pdf(predict_data[0]) *\
                   sp.stats.norm(theta[6][1], np.sqrt(sigma[6][1])).pdf(predict_data[1]) *\
                   sp.stats.norm(theta[6][2], np.sqrt(sigma[6][2])).pdf(predict_data[2]) *\
                   sp.stats.norm(theta[6][3], np.sqrt(sigma[6][3])).pdf(predict_data[3]) *\
                   sp.stats.norm(theta[6][4], np.sqrt(sigma[6][4])).pdf(predict_data[4]) *\
                   sp.stats.norm(theta[6][5], np.sqrt(sigma[6][5])).pdf(predict_data[5]) *\
                   sp.stats.norm(theta[6][6], np.sqrt(sigma[6][6])).pdf(predict_data[6]) *\
                   sp.stats.norm(theta[6][7], np.sqrt(sigma[6][7])).pdf(predict_data[7]) *\
                   sp.stats.norm(theta[6][8], np.sqrt(sigma[6][8])).pdf(predict_data[8]) *\
                   sp.stats.norm(theta[6][9], np.sqrt(sigma[6][9])).pdf(predict_data[9]))]
    return likelihood

• 위의 함수를 사용하여 10개의 실수 변수들에 대한 모수를 계산하여 가능도를 구하는 반복문을 작성하였다.

%%time
total_num = np.array(Gaussian_df_X).shape[0]
gaussian_likelihood_matrix = []
percentage = 0
for num, predict_data in enumerate(np.array(Gaussian_df_X)):
    if (percentage != int(num / total_num * 100)) and (int(num / total_num * 100) in list(np.arange(10,101,10))):
        percentage = int(num / total_num * 100)
        print("완성도 {} %".format(percentage))
    likelihood = Gaussian_likelihood_cal(predict_data, theta, sigma)
    gaussian_likelihood_matrix.append(likelihood)

• 위에서 가우시안 나이브 베이즈 모형의 가능도를 구했으므로 이젠 베르누이 나이브 베이즈 모형의 가능도를 구할 것이다.

BNB = BernoulliNB()
fitted_BNB = BNB.fit(Bern_df_X, df_Y)

theta = np.exp(fitted_BNB.feature_log_prob_)
theta.shape

• 상대적으로 가우시안 나이브 베이즈 모형의 가능도를 계산하는 것보단 단순 연산으로 이루어져 있어 속도가 훨씬 빠르다.

%%time
bern_likelihood_matrix = []
for data in np.array(Bern_df_X):
    bern_likelihood_matrix.append(list(((theta ** data) * (1 - theta) ** (1 - data)).prod(axis=1)))

결과

CPU times: user 49.2 s, sys: 1.28 s, total: 50.5 s
Wall time: 50.4 s

likelihood = np.array(gaussian_likelihood_matrix) * np.array(bern_likelihood_matrix)
posterior = likelihood * np.exp(BNB.class_log_prior_)
prob = posterior / np.repeat(posterior.sum(axis=1)[:, np.newaxis], 7, axis=1)
result = np.argmax(prob, axis=1)

result_bern = fitted_BNB.predict(Bern_df_X)
result_gaussian = fitted_GNB.predict(Gaussian_df_X)

가우시안 나이브 베이즈 모형의 성능

print(classification_report(result_gaussian, df_Y))

결과

              precision    recall  f1-score   support

           1       0.67      0.63      0.65    225973
           2       0.66      0.73      0.69    255934
           3       0.65      0.50      0.56     46785
           4       0.47      0.41      0.44      3099
           5       0.22      0.18      0.20     11425
           6       0.31      0.33      0.32     16424
           7       0.28      0.27      0.28     21372

    accuracy                           0.63    581012
   macro avg       0.47      0.44      0.45    581012
weighted avg       0.63      0.63      0.63    581012

베르누이 나이브 베이즈 모형의 성능

print(classification_report(result_bern, df_Y))

결과

              precision    recall  f1-score   support

           1       0.67      0.63      0.65    225973
           2       0.66      0.73      0.69    255934
           3       0.65      0.50      0.56     46785
           4       0.47      0.41      0.44      3099
           5       0.22      0.18      0.20     11425
           6       0.31      0.33      0.32     16424
           7       0.28      0.27      0.28     21372

    accuracy                           0.63    581012
   macro avg       0.47      0.44      0.45    581012
weighted avg       0.63      0.63      0.63    581012

가우시안 나이브 베이즈 모형과 베르누이 나이브 베이즈 모형의 가능도를 곱해 확률을 계산한 모형의 성능

print(classification_report(result, df_Y))

결과

              precision    recall  f1-score   support

           1       0.03      0.25      0.06     26481
           2       0.93      0.48      0.63    554531
           3       0.00      0.00      0.00         0
           4       0.00      0.00      0.00         0
           5       0.00      0.00      0.00         0
           6       0.00      0.00      0.00         0
           7       0.00      0.00      0.00         0

    accuracy                           0.47    581012
   macro avg       0.14      0.10      0.10    581012
weighted avg       0.89      0.47      0.60    581012

• 위의 성능 보면 3가지 모형 다 성능이 좋지 않다는 것을 확인 할 수 있다. 이는 적절한 피처의 선택이 이루어지지 않은 모형이기 때문일 것이며, 또한 아래 그림에서와 같이 클래스간의 비율차이가 극심하게 차이가 나는데, 특히 1,2 클래스가 대다수를 이루고 있기 때문에 1, 2클래스에 대한 학습이 많이 된 결과라고 해석 할 수 있을 것이다. 이는 마지막 두 나이브 베이즈 모형의 성능을 보아도 확인 할 수 있다. 마지막 모형의 성능은 다른 클래스로 예측한 데이터는 존재하지 않고 오로지 1과 2로 예측을 했다.