NB樸素貝葉斯理論推導與三種常見模型
轉自:http://www.tuicool.com/articles/zEJzIbR
樸素貝葉斯(Naive Bayes)是一種簡單的分類演算法,它的經典應用案例為人所熟知:文字分類(如垃圾郵件過濾)。很多教材都從這些案例出發,本文就不重複這些內容了,而把重點放在理論推導(其實很淺顯,別被“理論”嚇到),三種常用模型及其編碼實現(Python)。
如果你對理論推導過程不感興趣,可以直接逃到三種常用模型及編碼實現部分,但我建議你還是看看理論基礎部分。
另外,本文的所有程式碼都可以在 這裡獲取
1. 樸素貝葉斯的理論基礎
樸素貝葉斯演算法是基於貝葉斯定理與特徵條件獨立假設的分類方法。
這裡提到的 貝葉斯定理
1.1 貝葉斯定理
先看什麼是 條件概率 。
P(A|B)表示事件B已經發生的前提下,事件A發生的概率,叫做事件B發生下事件A的條件概率。其基本求解公式為:
$$P(A|B)=\frac{P(AB)}{P(B)}$$
貝葉斯定理便是基於條件概率,通過P(A|B)來求P(B|A):
$$P(B|A)=\frac{P(A|B)P(B)}{P(A)}$$
順便提一下,上式中的分母P(A),可以根據全概率公式分解為:
$$ P(A)= \sum_{i=1}^{n}P(B_i)P(A|B_i)$$
1.2 特徵條件獨立假設
這一部分開始樸素貝葉斯的理論推導,從中你會深刻地理解什麼是特徵條件獨立假設。
給定訓練資料集(X,Y),其中每個樣本x都包括n維特徵,即$x=({x_1,x_2,x_3,…,x_n})$,類標記集合含有k種類別,即$y=({y_1,y_2,…,y_k})$。
如果現在來了一個新樣本x,我們要怎麼判斷它的類別?從概率的角度來看,這個問題就是給定x,它屬於哪個類別的概率最大。那麼問題就轉化為求解$P(y_1|x)$,$P(y_2|x)$…$P(y_k|x)$中最大的那個,即求後驗概率最大的輸出:
$$argmax_{y_k}P(y_k|x)$$
那$P(y_k|x)$怎麼求解?答案就是貝葉斯定理:
$$P(y_k|x)=\frac{P(x|y_k)P(y_k)}{P(x)}$$
根據全概率公式,可以進一步地分解上式中的分母:
$$P(y_k|x)=\frac{P(x|y_k)P(y_k)}{\sum_kP(x|y_k)P(y_k)} 【公式1】$$
這裡休息兩分鐘
先不管分母,分子中的$P(y_k)$是先驗概率,根據訓練集就可以簡單地計算出來。
而條件概率$P(x|y_k)=P(x_1,x_2,…,x_n|y_k)$,它的引數規模是指數數量級別的,假設第i維特徵$x_i$可取值的個數有 $S_i$ 個,類別取值個數為k個,那麼引數個數為:
$$k\prod_{i=1}^{n}S_i$$
這顯然不可行。 針對這個問題,樸素貝葉斯演算法對條件概率分佈作出了獨立性的假設 ,通俗地講就是說假設各個維度的特徵$x_1,x_2,…,x_n$互相獨立,在這個假設的前提上,條件概率可以轉化為:
$$\prod_{i=1}^{n}P(x_i|y_k) 【公式2】$$
這樣,引數規模就降到$\sum_{i=1}^{n}S_i k$
以上就是針對條件概率所作出的特徵條件獨立性假設,至此,先驗概率$P(y_k)$和條件概率$P(x|y_k)$的求解問題就都解決了,那麼我們是不是可以求解我們所要的後驗概率$P(y_k|x)$了?
這裡思考兩分鐘
答案是肯定的。我們繼續上面關於$P(y_k|x)$的推導,將【公式2】代入【公式1】得到:
$P(y_k|x)=\frac{P(y k)\prod {i=1}^{n}P(x_i|y k)}{\sum {k}P(y k)\prod {i=1}^{n}P(x_i|y_k)}$
於是樸素貝葉斯分類器可表示為:
$f(x)=argmax_{y_k} P(y k|x)=argmax {y_k} \frac{P(y k)\prod {i=1}^{n}P(x_i|y k)}{\sum {k}P(y k)\prod {i=1}^{n}P(x_i|y_k)}$
因為對所有的$y_k$,上式中的分母的值都是一樣的(為什麼?注意到全加符號就容易理解了),所以可以忽略分母部分,樸素貝葉斯分類器最終表示為:
$f(x)=argmax P(y k)\prod {i=1}^{n}P(x_i|y_k)$
關於$P(y_k)$,$P(x_i|y_k)$的求解,有以下三種常見的模型.
2. 三種常見的模型及程式設計實現
2.1 多項式模型
當特徵是離散的時候,使用多項式模型。多項式模型在計算先驗概率$P(y_k)$和條件概率$P(x_i|y_k)$時,會做一些 平滑處理 ,具體公式為:
$P(y k)=\frac{N {y_k}+\alpha}{N+k\alpha}$
N是總的樣本個數,k是總的類別個數,$N_{y_k}$是類別為$y_k$的樣本個數,$\alpha$是平滑值。
$P(x_i|y k)=\frac{N {y_k,x i}+\alpha}{N {y_k}+n\alpha}$
$N_{y_k}$是類別為$y k$的樣本個數,n是特徵的維數,$N {y_k,x_i}$是類別為$y_k$的樣本中,第i維特徵的值是$x_i$的樣本個數,$\alpha$是平滑值。
當$\alpha=1$時,稱作Laplace平滑,當$0<\alpha<1$時,稱作Lidstone平滑,$\alpha=0$時不做平滑。
如果不做平滑,當某一維特徵的值$x_i$沒在訓練樣本中出現過時,會導致$P(x_i|y_k)=0$,從而導致後驗概率為0。加上平滑就可以克服這個問題。
2.1.1 舉例
有如下訓練資料,15個樣本,2維特徵$X^1,X^2$,2種類別-1,1。給定測試樣本$x=(2,S)^{T}$,判斷其類別。
解答如下:
運用多項式模型,令$\alpha=1$
- 計算先驗概率
- 計算各種條件概率
- 對於給定的$x=(2,S)^{T}$,計算:
由此可以判定y=-1。
2.1.2 程式設計實現(基於Python,Numpy)
從上面的例項可以看到,當給定訓練集時,我們無非就是先計算出所有的先驗概率和條件概率,然後把它們存起來(當成一個查詢表)。當來一個測試樣本時,我們就計算它所有可能的後驗概率,最大的那個對應的就是測試樣本的類別,而後驗概率的計算無非就是在查詢表裡查詢需要的值。
我的程式碼就是根據這個思想來寫的。定義一個MultinomialNB類,它有兩個主要的方法:fit(X,y)和predict(X)。fit方法其實就是訓練,呼叫fit方法時,做的工作就是構建查詢表。predict方法就是預測,呼叫predict方法時,做的工作就是求解所有後驗概率並找出最大的那個。此外,類的建構函式__init__()中,允許設定$\alpha$的值,以及設定先驗概率的值。具體程式碼及如下:
"""
Created on 2015/09/06
@author: wepon (http://2hwp.com)
API Reference: http://scikit-learn.org/stable/modules/naive_bayes.html#naive-bayes
"""
import numpy as np
class MultinomialNB(object):
"""
Naive Bayes classifier for multinomial models
The multinomial Naive Bayes classifier is suitable for classification with
discrete features
Parameters
----------
alpha : float, optional (default=1.0)
Setting alpha = 0 for no smoothing
Setting 0 < alpha < 1 is called Lidstone smoothing
Setting alpha = 1 is called Laplace smoothing
fit_prior : boolean
Whether to learn class prior probabilities or not.
If false, a uniform prior will be used.
class_prior : array-like, size (n_classes,)
Prior probabilities of the classes. If specified the priors are not
adjusted according to the data.
Attributes
----------
fit(X,y):
X and y are array-like, represent features and labels.
call fit() method to train Naive Bayes classifier.
predict(X):
"""
def __init__(self,alpha=1.0,fit_prior=True,class_prior=None):
self.alpha = alpha
self.fit_prior = fit_prior
self.class_prior = class_prior
self.classes = None
self.conditional_prob = None
def _calculate_feature_prob(self,feature):
values = np.unique(feature)
total_num = float(len(feature))
value_prob = {}
for v in values:
value_prob[v] = (( np.sum(np.equal(feature,v)) + self.alpha ) /( total_num + len(values)*self.alpha))
return value_prob
def fit(self,X,y):
#TODO: check X,y
self.classes = np.unique(y)
#calculate class prior probabilities: P(y=ck)
if self.class_prior == None:
class_num = len(self.classes)
if not self.fit_prior:
self.class_prior = [1.0/class_num for _ in range(class_num)] #uniform prior
else:
self.class_prior = []
sample_num = float(len(y))
for c in self.classes:
c_num = np.sum(np.equal(y,c))
self.class_prior.append((c_num+self.alpha)/(sample_num+class_num*self.alpha))
#calculate Conditional Probability: P( xj | y=ck )
self.conditional_prob = {} # like { c0:{ x0:{ value0:0.2, value1:0.8 }, x1:{} }, c1:{...} }
for c in self.classes:
self.conditional_prob[c] = {}
for i in range(len(X[0])): #for each feature
feature = X[np.equal(y,c)][:,i]
self.conditional_prob[c][i] = self._calculate_feature_prob(feature)
return self
#given values_prob {value0:0.2,value1:0.1,value3:0.3,.. } and target_value
#return the probability of target_value
def _get_xj_prob(self,values_prob,target_value):
return values_prob[target_value]
#predict a single sample based on (class_prior,conditional_prob)
def _predict_single_sample(self,x):
label = -1
max_posterior_prob = 0
#for each category, calculate its posterior probability: class_prior * conditional_prob
for c_index in range(len(self.classes)):
current_class_prior = self.class_prior[c_index]
current_conditional_prob = 1.0
feature_prob = self.conditional_prob[self.classes[c_index]]
j = 0
for feature_i in feature_prob.keys():
current_conditional_prob *= self._get_xj_prob(feature_prob[feature_i],x[j])
j += 1
#compare posterior probability and update max_posterior_prob, label
if current_class_prior * current_conditional_prob > max_posterior_prob:
max_posterior_prob = current_class_prior * current_conditional_prob
label = self.classes[c_index]
return label
#predict samples (also single sample)
def predict(self,X):
#TODO1:check and raise NoFitError
#ToDO2:check X
if X.ndim == 1:
return self._predict_single_sample(X)
else:
#classify each sample
labels = []
for i in range(X.shape[0]):
label = self._predict_single_sample(X[i])
labels.append(label)
return labels
我們用上面舉的例子來檢驗一下,注意S,M,L我這裡用4,5,6替換:
import numpy as np
X = np.array([
[1,1,1,1,1,2,2,2,2,2,3,3,3,3,3],
[4,5,5,4,4,4,5,5,6,6,6,5,5,6,6]
])
X = X.T
y = np.array([-1,-1,1,1,-1,-1,-1,1,1,1,1,1,1,1,-1])
nb = MultinomialNB(alpha=1.0,fit_prior=True)
nb.fit(X,y)
print nb.predict(np.array([2,4]))#輸出-1
2.2 高斯模型
當特徵是連續變數的時候,運用多項式模型就會導致很多$P(x_i|y_k)=0$(不做平滑的情況下),此時即使做平滑,所得到的條件概率也難以描述真實情況。所以處理連續的特徵變數,應該採用高斯模型。
2.2.1 通過一個例子來說明:
下面是一組人類身體特徵的統計資料。
| 性別 | 身高(英尺) | 體重(磅) |腳掌(英寸)|
| :————-: |:————-:| :—–:|:—–:|
| 男 | 6 | 180 | 12 |
|男 | 5.92 | 190 | 11 |
|男 | 5.58 | 170 | 12 |
| 男 | 5.92 | 165 | 10 |
| 女 | 5 | 100 | 6 |
| 女 | 5.5 | 150 | 8 |
| 女 | 5.42 | 130 | 7 |
| 女 | 5.75 | 150 | 9|
已知某人身高6英尺、體重130磅,腳掌8英寸,請問該人是男是女?
根據樸素貝葉斯分類器,計算下面這個式子的值。
P(身高|性別) x P(體重|性別) x P(腳掌|性別) x P(性別)
這裡的困難在於,由於身高、體重、腳掌都是連續變數,不能採用離散變數的方法計算概率。而且由於樣本太少,所以也無法分成區間計算。怎麼辦?
這時,可以假設男性和女性的身高、體重、腳掌都是正態分佈,通過樣本計算出均值和方差,也就是得到正態分佈的密度函式。有了密度函式,就可以把值代入,算出某一點的密度函式的值。
比如,男性的身高是均值5.855、方差0.035的正態分佈。所以,男性的身高為6英尺的概率的相對值等於1.5789(大於1並沒有關係,因為這裡是密度函式的值,只用來反映各個值的相對可能性)。
對於腳掌和體重同樣可以計算其均值與方差。有了這些資料以後,就可以計算性別的分類了。
P(身高=6|男) x P(體重=130|男) x P(腳掌=8|男) x P(男)
= 6.1984 x e-9
P(身高=6|女) x P(體重=130|女) x P(腳掌=8|女) x P(女)
= 5.3778 x e-4
可以看到,女性的概率比男性要高出將近10000倍,所以判斷該人為女性。
- 總結
高斯模型假設每一維特徵都服從高斯分佈(正態分佈):
$P(x_i|y k)=\frac{1}{\sqrt{2\pi\sigma {y_k,i}^{2}}}e^{-\frac{(x i-\mu {y k,i})^{2}}{2 \sigma {y_k,i}^{2}}}$
$\mu_{y_k,i}$表示類別為$yk$的樣本中,第i維特徵的均值。
$\sigma
{y_k,i}^{2}$表示類別為$y_k$的樣本中,第i維特徵的方差。2.2.2 程式設計實現
高斯模型與多項式模型唯一不同的地方就在於計算 $ P( x_i | y_k) $,高斯模型假設各維特徵服從正態分佈,需要計算的是各維特徵的均值與方差。所以我們定義GaussianNB類,繼承自MultinomialNB並且過載相應的方法即可。程式碼如下:
#GaussianNB differ from MultinomialNB in these two method:
# _calculate_feature_prob, _get_xj_prob
class GaussianNB(MultinomialNB):
"""
GaussianNB inherit from MultinomialNB,so it has self.alpha
and self.fit() use alpha to calculate class_prior
However,GaussianNB should calculate class_prior without alpha.
Anyway,it make no big different
"""
#calculate mean(mu) and standard deviation(sigma) of the given feature
def_calculate_feature_prob(self,feature):
mu = np.mean(feature)
sigma = np.std(feature)
return (mu,sigma)
#the probability density for the Gaussian distribution
def_prob_gaussian(self,mu,sigma,x):
return ( 1.0/(sigma * np.sqrt(2 * np.pi)) *
np.exp( - (x - mu)**2 / (2 * sigma**2)) )
#given mu and sigma , return Gaussian distribution probability for target_value
def_get_xj_prob(self,mu_sigma,target_value):
return self._prob_gaussian(mu_sigma[0],mu_sigma[1],target_value)
2.3 伯努利模型
與多項式模型一樣,伯努利模型適用於離散特徵的情況,所不同的是,伯努利模型中每個特徵的取值只能是1和0(以文字分類為例,某個單詞在文件中出現過,則其特徵值為1,否則為0).
伯努利模型中,條件概率$P(x_i|y_k)$的計算方式是:
當特徵值$x_i$為1時,$P(x_i|y_k)=P(x_i=1|y_k)$;
當特徵值$x_i$為0時,$P(x_i|y_k)=1-P(x_i=1|y_k)$;
2.3.1 程式設計實現
伯努利模型和多項式模型是一致的,BernoulliNB需要比MultinomialNB多定義一個二值化的方法,用於將輸入的特徵二值化(1,0)。當然也可以直接採用MultinomialNB,但需要將輸入的特徵預先二值化。寫到這裡不想寫了,程式設計實現留給讀者。