1、貝葉斯介紹

我個人一直很喜歡演算法一類的東西，在我看來演算法是人類智慧的精華，其中蘊含著無與倫比的美感。而每次將學過的演算法應用到實際中，並解決了實際問題後，那種快感更是我在其它地方體會不到的。

      一直想寫關於演算法的博文，也曾寫過零散的兩篇，但也許是相比於工程性文章來說太小眾，並沒有引起大家的興趣。最近面臨畢業找工作，為了能給自己增加籌碼，決定再次複習演算法方面的知識，我決定趁這個機會，寫一系列關於演算法的文章。這樣做，主要是為了加強自己複習的效果，我想，如果能將複習的東西用自己的理解寫成文章，勢必比單純的讀書做題掌握的更牢固，也更能觸發自己的思考。如果能有感興趣的朋友從中有所收穫，那自然更好。

      這個系列我將其命名為“演算法雜貨鋪”，其原因就是這些文章一大特徵就是“雜”，我不會專門討論堆疊、連結串列、二叉樹、查詢、排序等任何一本資料結構教科書都會講的基礎內容，我會從一個“專題”出發，如概率演算法、分類演算法、NP問題、遺傳演算法等，然後做一個引申，可能會涉及到演算法與資料結構、離散數學、概率論、統計學、運籌學、資料探勘、形式語言與自動機等諸多方面，因此其內容結構就像一個雜貨鋪。當然，我會竭盡所能，儘量使內容“雜而不亂”。

1.1、摘要

      貝葉斯分類是一類分類演算法的總稱，這類演算法均以貝葉斯定理為基礎，故統稱為貝葉斯分類。本文作為分類演算法的第一篇，將首先介紹分類問題，對分類問題進行一個正式的定義。然後，介紹貝葉斯分類演算法的基礎——貝葉斯定理。最後，通過例項討論貝葉斯分類中最簡單的一種：樸素貝葉斯分類。

1.2、分類問題綜述

      對於分類問題，其實誰都不會陌生，說我們每個人每天都在執行分類操作一點都不誇張，只是我們沒有意識到罷了。例如，當你看到一個陌生人，你的腦子下意識判斷TA是男是女；你可能經常會走在路上對身旁的朋友說“這個人一看就很有錢、那邊有個非主流”之類的話，其實這就是一種分類操作。

      從數學角度來說，分類問題可做如下定義：

      已知集合：和，確定對映規則，使得任意有且僅有一個使得成立。（不考慮模糊數學裡的模糊集情況）

      其中C叫做類別集合，其中每一個元素是一個類別，而I叫做項集合，其中每一個元素是一個待分類項，f叫做分類器。分類演算法的任務就是構造分類器f。

      這裡要著重強調，分類問題往往採用經驗性方法構造對映規則，即一般情況下的分類問題缺少足夠的資訊來構造100%正確的對映規則，而是通過對經驗資料的學習從而實現一定概率意義上正確的分類，因此所訓練出的分類器並不是一定能將每個待分類項準確對映到其分類，分類器的質量與分類器構造方法、待分類資料的特性以及訓練樣本數量等諸多因素有關。

      例如，醫生對病人進行診斷就是一個典型的分類過程，任何一個醫生都無法直接看到病人的病情，只能觀察病人表現出的症狀和各種化驗檢測資料來推斷病情，這時醫生就好比一個分類器，而這個醫生診斷的準確率，與他當初受到的教育方式（構造方法）、病人的症狀是否突出（待分類資料的特性）以及醫生的經驗多少（訓練樣本數量）都有密切關係。

1.3、貝葉斯分類的基礎——貝葉斯定理

      每次提到貝葉斯定理，我心中的崇敬之情都油然而生，倒不是因為這個定理多高深，而是因為它特別有用。這個定理解決了現實生活裡經常遇到的問題：已知某條件概率，如何得到兩個事件交換後的概率，也就是在已知P(A|B)的情況下如何求得P(B|A)。這裡先解釋什麼是條件概率：

      表示事件B已經發生的前提下，事件A發生的概率，叫做事件B發生下事件A的條件概率。其基本求解公式為：。

      貝葉斯定理之所以有用，是因為我們在生活中經常遇到這種情況：我們可以很容易直接得出P(A|B)，P(B|A)則很難直接得出，但我們更關心P(B|A)，貝葉斯定理就為我們打通從P(A|B)獲得P(B|A)的道路。

      下面不加證明地直接給出貝葉斯定理：

1.4、樸素貝葉斯分類

1.4.1、樸素貝葉斯分類的原理與流程

      樸素貝葉斯分類是一種十分簡單的分類演算法，叫它樸素貝葉斯分類是因為這種方法的思想真的很樸素，樸素貝葉斯的思想基礎是這樣的：對於給出的待分類項，求解在此項出現的條件下各個類別出現的概率，哪個最大，就認為此待分類項屬於哪個類別。通俗來說，就好比這麼個道理，你在街上看到一個黑人，我問你你猜這哥們哪裡來的，你十有八九猜非洲。為什麼呢？因為黑人中非洲人的比率最高，當然人家也可能是美洲人或亞洲人，但在沒有其它可用資訊下，我們會選擇條件概率最大的類別，這就是樸素貝葉斯的思想基礎。

      樸素貝葉斯分類的正式定義如下：

      1、設為一個待分類項，而每個a為x的一個特徵屬性。

      2、有類別集合。

      3、計算。

      4、如果，則。

      那麼現在的關鍵就是如何計算第3步中的各個條件概率。我們可以這麼做：

      1、找到一個已知分類的待分類項集合，這個集合叫做訓練樣本集。

      2、統計得到在各類別下各個特徵屬性的條件概率估計。即。

      3、如果各個特徵屬性是條件獨立的，則根據貝葉斯定理有如下推導：

      因為分母對於所有類別為常數，因為我們只要將分子最大化皆可。又因為各特徵屬性是條件獨立的，所以有：

      根據上述分析，樸素貝葉斯分類的流程可以由下圖表示（暫時不考慮驗證）：

      可以看到，整個樸素貝葉斯分類分為三個階段：

      第一階段——準備工作階段，這個階段的任務是為樸素貝葉斯分類做必要的準備，主要工作是根據具體情況確定特徵屬性，並對每個特徵屬性進行適當劃分，然後由人工對一部分待分類項進行分類，形成訓練樣本集合。這一階段的輸入是所有待分類資料，輸出是特徵屬性和訓練樣本。這一階段是整個樸素貝葉斯分類中唯一需要人工完成的階段，其質量對整個過程將有重要影響，分類器的質量很大程度上由特徵屬性、特徵屬性劃分及訓練樣本質量決定。

      第二階段——分類器訓練階段，這個階段的任務就是生成分類器，主要工作是計算每個類別在訓練樣本中的出現頻率及每個特徵屬性劃分對每個類別的條件概率估計，並將結果記錄。其輸入是特徵屬性和訓練樣本，輸出是分類器。這一階段是機械性階段，根據前面討論的公式可以由程式自動計算完成。

      第三階段——應用階段。這個階段的任務是使用分類器對待分類項進行分類，其輸入是分類器和待分類項，輸出是待分類項與類別的對映關係。這一階段也是機械性階段，由程式完成。

1.4.2、估計類別下特徵屬性劃分的條件概率及Laplace校準

這一節討論P(a|y)的估計。

      由上文看出，計算各個劃分的條件概率P(a|y)是樸素貝葉斯分類的關鍵性步驟，當特徵屬性為離散值時，只要很方便的統計訓練樣本中各個劃分在每個類別中出現的頻率即可用來估計P(a|y)，下面重點討論特徵屬性是連續值的情況。

      當特徵屬性為連續值時，通常假定其值服從高斯分佈（也稱正態分佈）。即：

      而

      因此只要計算出訓練樣本中各個類別中此特徵項劃分的各均值和標準差，代入上述公式即可得到需要的估計值。均值與標準差的計算在此不再贅述。

      另一個需要討論的問題就是當P(a|y)=0怎麼辦，當某個類別下某個特徵項劃分沒有出現時，就是產生這種現象，這會令分類器質量大大降低。為了解決這個問題，我們引入Laplace校準，它的思想非常簡單，就是對沒類別下所有劃分的計數加1，這樣如果訓練樣本集數量充分大時，並不會對結果產生影響，並且解決了上述頻率為0的尷尬局面。

1.4.3、樸素貝葉斯分類例項：檢測SNS社群中不真實賬號

      下面討論一個使用樸素貝葉斯分類解決實際問題的例子，為了簡單起見，對例子中的資料做了適當的簡化。

      這個問題是這樣的，對於SNS社群來說，不真實賬號（使用虛假身份或使用者的小號）是一個普遍存在的問題，作為SNS社群的運營商，希望可以檢測出這些不真實賬號，從而在一些運營分析報告中避免這些賬號的干擾，亦可以加強對SNS社群的瞭解與監管。

      如果通過純人工檢測，需要耗費大量的人力，效率也十分低下，如能引入自動檢測機制，必將大大提升工作效率。這個問題說白了，就是要將社群中所有賬號在真實賬號和不真實賬號兩個類別上進行分類，下面我們一步一步實現這個過程。

      首先設C=0表示真實賬號，C=1表示不真實賬號。

      1、確定特徵屬性及劃分

      這一步要找出可以幫助我們區分真實賬號與不真實賬號的特徵屬性，在實際應用中，特徵屬性的數量是很多的，劃分也會比較細緻，但這裡為了簡單起見，我們用少量的特徵屬性以及較粗的劃分，並對資料做了修改。

      我們選擇三個特徵屬性：a1：日誌數量/註冊天數，a2：好友數量/註冊天數，a3：是否使用真實頭像。在SNS社群中這三項都是可以直接從資料庫裡得到或計算出來的。

      下面給出劃分：a1：{a<=0.05, 0.05<a<0.2, a>=0.2}，a1：{a<=0.1, 0.1<a<0.8, a>=0.8}，a3：{a=0（不是）,a=1（是）}。

      2、獲取訓練樣本

      這裡使用運維人員曾經人工檢測過的1萬個賬號作為訓練樣本。

      3、計算訓練樣本中每個類別的頻率

      用訓練樣本中真實賬號和不真實賬號數量分別除以一萬，得到：

      4、計算每個類別條件下各個特徵屬性劃分的頻率

      5、使用分類器進行鑑別

      下面我們使用上面訓練得到的分類器鑑別一個賬號，這個賬號使用非真實頭像，日誌數量與註冊天數的比率為0.1，好友數與註冊天數的比率為0.2。

      可以看到，雖然這個使用者沒有使用真實頭像，但是通過分類器的鑑別，更傾向於將此賬號歸入真實賬號類別。這個例子也展示了當特徵屬性充分多時，樸素貝葉斯分類對個別屬性的抗干擾性。

1.5、分類器的評價

      雖然後續還會提到其它分類演算法，不過這裡我想先提一下如何評價分類器的質量。

      首先要定義，分類器的正確率指分類器正確分類的專案佔所有被分類專案的比率。

      通常使用迴歸測試來評估分類器的準確率，最簡單的方法是用構造完成的分類器對訓練資料進行分類，然後根據結果給出正確率評估。但這不是一個好方法，因為使用訓練資料作為檢測資料有可能因為過分擬合而導致結果過於樂觀，所以一種更好的方法是在構造初期將訓練資料一分為二，用一部分構造分類器，然後用另一部分檢測分類器的準確率。

2、MATLAB的程式程式碼：

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%By Shelley from NCUT，April 14th 2011
%Email:[email protected]
%此程式利用貝葉斯分類演算法，首先對兩類樣本進行訓練，
%進而可在螢幕上任意取點，程式可輸出屬於第一類，還是第二類
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
close all
 
%讀入兩類訓練樣本資料
load data
%求兩類訓練樣本的均值和方差
u1=mean(Sample1);
u2=mean(Sample2);
sigm1=cov(Sample1); 
sigm2=cov(Sample2);
%計算兩個樣本的密度函式並顯示
x=-20:0.5:40;
y= -20:0.5:20;
[X,Y] = meshgrid(x,y);
F1 = mvnpdf([X(:),Y(:)],u1,sigm1);
F2 = mvnpdf([X(:),Y(:)],u2,sigm2);
P1=reshape(F1,size(X));
P2=reshape(F2,size(X));
figure(2)
surf(X,Y,P1)
hold on
surf(X,Y,P2)
shading interp
colorbar
title('條件概率密度函式曲線');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%以下為測試部分
%利用ginput隨機選取螢幕上的點（可連續取10個點）
%程式可根據點的位置自動地顯示出屬於那個類
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pw1=0.4;pw2=0.6;
 
figure(1)
plot(Sample1(:,1),Sample1(:,2),'r.')
hold on
plot(Sample2(:,1),Sample2(:,2),'b.')
 
for i=1:10
    [u,v]=ginput(1);
    plot(u,v,'m*');
    P1=pw1*mvnpdf([u,v],u1,sigm1);
    P2=pw2*mvnpdf([u,v],u2,sigm2);
    hold all
    if(P1>P2)
         disp('it belong to the first class');
    else
         disp('it belong to the second class');
    end;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

3、MATLAB進行二維正太分佈圖

   為什麼要介紹這個原因，是因為最小風險貝葉斯里面的那個P(X|Wi)是二維正態分佈概率密度，其x軸是每個樣本的屬性，y軸是每個樣本中屬性對於的個數。其是程正態分佈的，其中的期望值是所有屬性對應個數加起來的平均值。則其二維正態分佈概率密度公式是：

clc
clear
close all

mu=[0,0];% 均值向量
Sigma=[1 0.8;0.8 1];% 協方差矩陣
[X,Y]=meshgrid(-3:0.1:3,-3:0.1:3);%在XOY面上，產生網格資料
p=mvnpdf([X(:) Y(:)],mu,Sigma);%求取聯合概率密度，相當於Z軸
p=reshape(p,size(X));%將Z值對應到相應的座標上

figure
set(gcf,'Position',get(gcf,'Position').*[1 1 1.3 1])

subplot(2,3,[1 2 4 5])
surf(X,Y,p),axis tight,title('二維正態分佈圖')
subplot(2,3,3)
surf(X,Y,p),view(2),axis tight,title('在XOY面上的投影'

3、最小風險貝葉斯分類器的構建

其程式碼是：

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%功能：演示貝葉斯學習演算法在計算機視覺中的應用
%基於貝葉斯學習實現目標分類；

%其是真實訓練跟測試資料都是k_d-k_s=120-11=109
%環境：Win7，Matlab2012b
%Modi: NUDT-VAP
%時間：2014-02-04
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all
clear all
%載入背景和行人樣本
load human;
load background;
%隨機劃分訓練樣本和測試樣本
for j=1:10
k_s=11;
k_d=120;
hm_tr(250,252)=0;    %這隻初始化了一個元素
hm_test(250,252)=0;
% ceshi=ones(size(hm_test));
[hm_tr,hm_test]=randQ(human);
bg_tr(250,252)=0;
bg_test(250,252)=0;
[bg_tr,bg_test]=randQ(background);   %這裡的測試資料值跟訓練的資料值是不一樣的，所以這裡要使用區域範圍來把測試資料歸為訓練資料
%設計最小風險貝葉斯分類器
M_hm=hm_tr(1,k_s:k_d);
for i=2:250
M_hm=M_hm+(hm_tr(i,k_s:k_d));  %每一列的資料進行累加
end
M_hm=M_hm/250;                 %累加求平均值
M_bg=bg_tr(1,k_s:k_d);
for i=2:250
M_bg=M_bg+(bg_tr(i,k_s:k_d));
end
M_bg=M_bg/250; 
E_hm(k_d-k_s+1,k_d-k_s+1)=0;
E_bg(k_d-k_s+1,k_d-k_s+1)=0;
for i=1:250
    E_hm=E_hm+(hm_tr(i,k_s:k_d)-M_hm)'*(hm_tr(i,k_s:k_d)-M_hm);
    E_bg=E_bg+(bg_tr(i,k_s:k_d)-M_bg)'*(bg_tr(i,k_s:k_d)-M_bg);
end
E_hm=E_hm/250;  %這裡不是歸一化
E_bg=E_bg/250;

%分類測試
%設計判決函式並分類
flag_hm(250)=0;
flag_bg(250)=0;
for i=1:250
    dk_hm=log(0.5)+log(0.95)-log(abs(det(E_bg)))/2-((hm_test(i,k_s:k_d)-M_bg)*inv(E_bg)*(hm_test(i,k_s:k_d)-M_bg)')/2;   %其中的0.5是背景判為行人的風險；1.5為行人判為背景的風險；出現背景的先驗概率為0.95；出現行人的先驗概率為0.05；
    dk_bg=log(1.5)+log(0.05)-log(abs(det(E_hm)))/2-((hm_test(i,k_s:k_d)-M_hm)*inv(E_hm)*(hm_test(i,k_s:k_d)-M_hm)')/2;
    if(dk_hm<=dk_bg)
        flag_hm(i)=1;
    else flag_hm(i)=0;
    end
end
s2_hm=sum(flag_hm);
for i=1:250
    %因為只有方陣才有逆矩陣
    h=inv(E_bg);
    dk_hm=log(0.5)+log(0.95)-log(abs(det(E_bg)))/2-((bg_test(i,k_s:k_d)-M_bg)*inv(E_bg)*(bg_test(i,k_s:k_d)-M_bg)')/2;
    dk_bg=log(1.5)+log(0.05)-log(abs(det(E_hm)))/2-((bg_test(i,k_s:k_d)-M_hm)*inv(E_hm)*(bg_test(i,k_s:k_d)-M_hm)')/2;
    if(dk_hm>dk_bg)
        flag_bg(i)=1;
    else flag_bg(i)=0;
    end
end
%統計分類結果
s2_bg=sum(flag_bg);
err_hg(j)=(250-s2_hm)/250;
err_gh(j)=(250-s2_bg)/250;
end
ERR_HG=0;
ERR_GH=0;
for j=1:10
    ERR_HG=ERR_HG+err_hg(j);
    ERR_GH=ERR_GH+err_gh(j);
end
ERR_HG=ERR_HG/10
ERR_GH=ERR_GH/10