模式識別導論大作業

一、      K均值聚類

1.      功能描述：

利用K-均值演算法將150個模式樣本分成3類別。分別計算最後演算法所用的迭代次數，最終聚類中心以及每個類別中對應模式樣本的序號。

2.      帶註釋的原始碼

#include "stdio.h"

#include "math.h"

typedef struct mydata{

    float data[4];   //儲存原始資料

    int index;

}MyData;

MyData myData[150];

float classdata[3][150][4];

int saveindex[3][150];  //儲存下表

float center[3][4];  //儲存三個聚類中心

float newcenter[3][4];//儲存三個新的聚類中心

int findMinCenter(float*temp){//找到當前資料的類別

    float dis[3]={0,0,0};

    float min = 100000;

    int returnk=0;

    for(int k=0;k<3;k++){

        float mytemp = 0;

        for(int m=0;m<4;m++)

            mytemp += (temp[m]-center[k][m])*(temp[m]-center[k][m]);

        dis[k] = sqrt(mytemp);

        //printf("%f\n",dis[k]);

        if(dis[k]<min){

            min = dis[k];

            returnk = k;

        }}  return returnk;

}

void counterNewCenter(floattempdata[][4],int len,int num){

    float temp[4] = {0,0,0,0}; //計算新的聚類中心方法

    for(int j=0;j<len;j++){

        for(int i=0;i<4;i++){

            temp[i] += tempdata[j][i];

        }   }

    for(int k=0;k<4;k++){

        newcenter[num][k] = temp[k]/len;

        printf("k = %f\n",newcenter[num][k]);

}}

int main(){

    int k=0;

    int kt[3] = {0,0,0};

    int i=0,u=0;

    FILE *fp;

    if ((fp = fopen("D:\\Iris.txt", "r")) ==NULL){  //開啟檔案

        printf("開啟檔案失敗\n");

        return 0;

    }   float numtemp;

    while (!feof(fp)){     //讀取資料

        fscanf(fp, "%f", &numtemp);

        myData[u].data[i++] = numtemp;

        if(i==4){

            u++;

            i=0;}   }

    fclose(fp);//關閉檔案 

    for(k=0;k<4;k++){

        center[0][k] = myData[0].data[k];

        center[1][k] = myData[1].data[k];

        center[2][k] = myData[2].data[k];

    }//設定初始聚類中心

    int s=0;int num = 0;i=0;//聚類所用的迭代次數

    while(i<150){

        i++;

        k = findMinCenter(myData[i].data);

        //printf("k = %d\n",k);

        saveindex[k][kt[k]] = i;

        for(s=0;s<4;s++){

            classdata[k][kt[k]][s] = myData[i].data[s];}//每次將資料分到不同的類中

        kt[k]++;if (i == 150){

            printf("******************");

            for (int m = 0; m<3; m++)//計算新的聚類中心

                counterNewCenter(classdata[m],kt[m],m);

            bool flag = true;//標示前後聚類中心是否相等的布林型別

            for(int n=0;n<3;n++){

                for(int t=0;t<4;t++){

                   if(newcenter[n][t] != center[n][t]){

                       flag = false;   break;}//if

                   if(!flag)   break;}//for   

            }//for

            if(!flag){  i = 0;num++; //如果前後兩次的聚類中心不相等，則繼續分類

                for(int y=0;y<3;y++){

                   kt[y] = 0;

                   for(int w=0;w<4;w++)

                       center[y][w] = newcenter[y][w];//構建新的聚類中心

                   }//for

            }//if(!flag)

        }//if (i == 150)

            }//while

    printf("%d\n",num); //輸出聚類所用的迭代次數

    if ((fp = fopen("D:\\result.txt", "wr")) ==NULL){  //建立result.txt檔案

        printf("建立檔案失敗\n");

        return 0;   }

    for(i=0;i<3;i++){

        for(u=0;u<kt[i];u++){

            printf("%d  ",saveindex[i][u]);//將3組資料的下標寫入檔案儲存

            fprintf(fp,"%d ",saveindex[i][u]);    }

        fprintf(fp,"\n\n\n");

    }fclose(fp);

    return 0;

}

 

3.      分類結果

聚類所用的迭代次數: 11 

第一組：

聚類中心 ( 5.004082  , 3.426531 , 1.463265  , 0.246939 )

模式樣本序號：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

第二組：

聚類中心 (5.883606  ， 2.740984 ， 4.409837 ， 1.434426 )

模式樣本序號：52  54 55  56  57 58  59  60 61  62  63 64  65  66 67  68  69 70  71  72 73  74  75 76  77  79 80  81  82 83  84  85 86  87  88 89  90  91 92  93  94 95  96  97 98  99  100 102  107  114 115  120  122 124  127  128 134  139  143 147  150

第三組：

聚類中心 (6.853845    3.076923   5.715384    2.053846 )

模式樣本序號：51  53 78  101  103 104  105  106 108  109  110 111  112  113 116  117  118 119  121  123 125  126  129 130  131  132 133  135  136 137  138  140 141  142  144 145  146  148 149

4.      演算法分析

檔案中共有150個模式樣本，要求分成3個模式類別，具體步驟如下：

1.  讀取檔案中資料，由於樣本維數為4。所以建立一個二維組用fscanf格式化讀取檔案中的樣本資料。

2.  選取3個初始樣本中心，這裡我選擇的是樣本前三個資料作為初始樣本中。用center資料儲存。

3.  建立while迴圈，迴圈內依次將樣本資料按照最短距離原則分類到3個聚類中心中的某一個。這裡我用findMinCenter()方法返回當前資料到哪個聚類中心最近的序號。除此之外，用saveindex陣列分別記錄每個聚類中心中所包含模式樣本的序號，方便最後將結果輸出。

4.  用counterNewCenter()計算當前類別的新的聚類中心。設定布林型別的標識flag，初始為true。將3個新的聚類中心，分別和3個原始聚類中心相比較，如有有一個不相同，則設定flag為false，置樣本索引i=0，繼續迴圈。否則跳出迴圈。得到分類結果。

二、      線性分類器設計（20分）

以第一題的分類結果作為樣本集，首先選取訓練集與測試集（訓練集大概是整體樣本的2/3），請分別給出三個類別的訓練集與測試集包含的樣本編號：

第一類訓練集樣本編號（樣本個數26個）: 51  53  78 101  103  104 105  106  108 109  110  111 112  113  116 117  118  119 121  123  125 126  129  130 131  132

第二類訓練集樣本編號（樣本個數40個）: 52 54  55  56 57  58  59 60  61  62 63  64  65 66  67  68 69  70  71 72  73  74 75  76  77 79  80  81 82  83  84 85  86  87 88  89  90 91  92  93

第三類訓練集樣本編號（樣本個數33個）: 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

第一類測試集樣本編號: 133  135 136  137  138 140  141  142 144  145  146 148  149

第二類測試集樣本編號: 94  95 96  97  98 99  100  102 107  114  115 120  122  124 127  128  134 139  143  147 150

第三類測試集樣本編號: 34  35  36 37  38  39 40  41  42 43  44  45 46  47  48 49  50

1.  設定分類規則（一對多還是一對一），並利用迭代法在訓練集上設計線性分類器（給出懲罰係數），給出判別函式。同時將測試集中的資料代入判別函式，給出測試的準確率；

懲罰係數: factor = 0.1

最終計算的得出的權向量為:

W1  =  (-15.1156  -13.7006   -1.32033   0.448794  -112.199)T

W2  =  (-34.1219  -19.1983   -14.6433   -3.18996  101.799)T

W3  =  (-18.77        1.07      -28.8         -10.58         0.5)T

判別函式:

D1(X)= -15.1156x1-13.7006x2-1.32033x3+0.448794x4-112.199

D2(X)= -34.1219x1-19.1983x2-14.6433x3+-3.18996x4+101.799

D3(X)= -18.77x1+1.07x2 -28.8x3 -10.58x4+0.5

測試的準確率:  96.0784%

2.  設定分類規則（一對多還是一對一），並利用Fisher法在訓練集上設計線性分類器，給出判別函式（給出各類的類間離散度矩陣與類內離散度矩陣等引數）。同時將測試集中的資料代入判別函式，給出測試的準確率；

類的類間離散度矩陣:

Sb =

類內離散度矩陣:

第一類類內離散矩陣 :

第一類類內離散矩陣:

第一類類內離散矩陣:

測試的準確率:  86.28%

3.  對上述的兩種演算法進行分析

感知器演算法通過賞罰原則依據每次對訓練集的訓練不斷修正判別函式的權向量，當分類器發生錯誤分類的時候對分類器進行“罰”,即對權向量進行修改，當感知器正確分類的時候對分類器進行“賞”，對全向量不進行修改。這樣經過迭代計算後，通過訓練集的訓練得到最優的判別函式的權向量。可以從程式上看出感知器演算法的實現是十分簡單清晰的，而且對測試集的判斷正確率也十分理想，高達96%以上。但是感知器一般只用於小樣本資料的學習和識別，對於大樣本資料感知器演算法就顯出劣勢了。相對於感知器演算法的fisher分類演算法在樣本的判別正確率上要遠低於感知器演算法，而且fisher演算法的內部實現相對於感知器演算法要複雜得多，這大多是由於矩陣計算相當消耗計算機資源。Fisher演算法在判別上劣於感知器演算法主要是樣本的訓練集並不能達到fisher演算法訓練的標準。另一個在樣本的不平衡條件下，會導致fisher演算法的協方差矩陣不好估計，從而導致矩陣計算的誤差。從綜合性能和判別結果上分析，感知器演算法對於小樣本的判別要更適合fisher演算法，而且在時間和演算法複雜度上都要更好些。

三、      貝葉斯決策（20分）

保持第二題中已選擇的訓練集與測試集不變，將所有的樣本從4維降低為兩維（隨機選取兩維），將訓練集的所有樣本點在二維座標系下標註出來，注意不同的類別用不同形狀的點加以區分。

選取樣本的第一維和第二維。

1.      利用貝葉斯決策與訓練集設計分類器，寫出判別函式與判別規則。

判別函式：d(1)= 42.25*x1 - 4.126*x2+13.19x1*x2-9.01*x2*x2-8.73x1*x1-97.21

                d(2)=19.1*x1 + 10.79*x2-2.489x1*x1-6.145x2*x2+3.888x1*x2-70.71

                d(3)= 9.813*x1 + 14.42*x2 -1.301x1*x1+5.26x2*x2+2.502x1*x2-5207

判別規則：

如果d(1)>d(2) 且d(1)>d(3)則樣本屬於第一類

              如果d(2)>d(1)且d(2)>d(3)則樣本屬於第二類

如果d(3)>d(1) 且d(3)>d(2)則樣本屬於第三類

2.      將判別函式作用於測試集，在二維座標系下將測試集的資料標註出來，注意不同的類別用不同形狀的點加以區分，用深淺不同的灰度表示錯分與正確分類的樣本點，給出分類的準確率。

測試樣本個數:  51

錯誤分類樣本個數:9

分類準確率:82.35%

3.      實驗結果分析

貝葉斯分類決策是通過統計學理論知識中樣本的概率問題對樣本進行分類判別的。通過樣本資料計算樣本的後驗概率的大小決定樣本最終屬於的類別。對於正態分佈的貝葉斯判別決策更是需要特徵空間中某一類樣本較多的分佈在其均值附近，遠離均值樣點的個數較少。這樣決策才算合理。這使得貝葉斯分類決策十分依賴樣本資料。所以對於本實驗中4維的資料降到2維，如果選擇的維數不同會使得實驗結果有較大差距。根據實驗結果可以看出某些樣本中的兩個資料並不能準確的反應這個樣本的類別，這導致貝葉斯判別失誤。所以在選擇樣本維數的時候要儘量選擇容易區分樣本。

四、      特徵選擇（20分）

1. 保持第二題中已選擇的訓練集與測試集不變，利用特徵選擇演算法，將樣本從4維特徵空間降至2維特徵空間。給出特徵選擇的依據與演算法過程。並將訓練集的所有樣本點在二維座標系下標註出來，注意不同的類別用不同形狀的點加以區分。

2.  特徵選擇的依據：K-L變換

3.  演算法過程：

4.      給定樣本X後，首先統計並計算其協方差矩陣C_X ；

5.      求C_X的特徵根，選取前M個最大的特徵根對應的特徵向量，這些特徵向量構成K-L變換的正交陣 ；

6.      對X做K-L變換後得到M維向量Y，作為表示原始樣本X的特徵

7.      選取樣本的第二三維作為資料樣本:

8.      

9.      依然利用貝葉斯決策和訓練集設計分類器，寫出判別函式。

判別函式:

d(1)=27.59*x1 + 47.52*x2 - 0.5*x1*(8.048*x1 -0.161*x2) + 0.5*x2*(0.161*x1 - 32.51*x2) - 80.46

d(2)= 13.87*x1 + 16.25*x2 + 0.5*x2*(5.176*x1 - 7.096*x2) -0.5*x1*(13.03*x1 - 5.176*x2) - 52.48

d(3)= 15.63*x1 + 10.69*x2 + 0.5*x2*(2.429*x1 - 3.153*x2) -0.5*x1*(9.992*x1 - 2.429*x2) - 52.22

10.  將判別函式作用於測試集，在二維座標系下將測試集的資料標註出來，注意不同的類別用不同形狀的點加以區分，用深淺不同的灰度表示錯分與正確分類的樣本點，給出分類的準確率

測試樣本個數:  51

錯誤分類樣本個數:4

分類準確率: 92.16%

實驗結果分析

從實驗結果來看分類效果還是很理想的，對樣本資料進行K-L變換後，再進行貝葉斯分類可以發現樣本的分類效果比原來更加高效了。但是K-L最佳變換效能雖好，但實現起來不易。樣本不同，協方差矩陣就不同，因此對應的變換T也不一樣。為了得到最佳變換，每來一個樣本就要重複上述步驟，運算相當繁瑣。

五、      系統設計（12分）

1.      任務描述：設計一個模式識別系統，可以識別筆、直尺與橡皮。

通過採集筆、直尺與橡皮長寬高等資料先利用訓練集正態分佈貝葉斯決策設計判別函式，然後用設計好的分類器對測試樣本進行識別分類。

2.      基本要求：按照統計類模式識別的基本步驟，從資料採集開始，到特徵提取、分類器設計，直至測試結果，給出每一步的具體實現過程。

(1).數據採集：通過對筆、直尺與橡皮的圖片提取圖片中物體的長、寬、高或者介面半徑的畫素個數。

(2).特徵提取: 這裡提取樣本的長寬比和最小截面面積來作為樣本特徵資料。其中樣本長寬比等於樣本長度的畫素個數除以樣本寬度的畫素個數。對於圓形的樣本用用截面半徑取代樣本高度。對最後的樣本資料進行歸一化。

二維(x1,x2)T，x1表示樣本長寬比，x2表示樣本最小橫截面積

筆	直尺	橡皮
(10.0 , 2.25)	(7.0 , 0.3)	(1.8 ,6.6)
(12.2 , 2.0)	(9.4 , 0.25)	(2.3 , 8.0)
(9.8 , 1.8)	(6.0 , 0.4)	(1.5 , 4.8)
(19.0 , 2.0)	(10.5 , 0.32)	(2.0 , 5.8)
(11.7 , 1.44)	(8.5 , 0.34)	(3.4 , 4.2)
(14.5 , 2.5)	(6.2 , 0.5)	(2.8 , 3.8)
(10.8 , 1.9)	(5.8 , 0.28)	(1.2 , 5.0)

(3).分類器設計:通過正態分佈貝葉斯決策對訓練集樣本進行訓練，樣本資料為二維(x1,x2)T，x1表示樣本長寬比，x2表示樣本最小橫截面積。計算各個類別的類內散佈矩陣和類間散佈矩陣。然後獲得三個類別的判別函式D(1),D(2),D(3)分別為筆、直尺和橡皮的判別函式。如果D(1)>D(2)&D(1)>D(3)，則該樣本屬於筆；如果D(1)>D(2)&D(1)>D(3)，則該樣本屬於直尺；如果D(1)>D(2)&D(1)>D(3)，則該樣本屬於橡皮。

              (4).測試結果：

                     筆: 100%

                     直尺: 100%

                     橡皮: 100%

六、      心得體會（8分）

通過這一個多月的學習，模式識別這門課給我帶來的收穫還是很多的。特別是這個在完成這個大作業的過程中自己學習到了很多東西。第一個是書本上的理論內容，我梳理了下自己總結的知識點。模式識別過程實際上是計算機通過對採集到的資料進行特徵提取後，用模式分類演算法構造判別函式（也就是分類器）.大作業中有對樣本進行分類的k均值演算法按最小距離原則依靠聚類中心分類。然後是感知器演算法，通過迭代不斷修改權向量，最後得到最優的判別函式。Fisher演算法對多維資料進行降維處理，找出最優的投影面。貝葉斯決策利用概率論中的知識，通過計算類間離散度和類內離散度，最後得到權向量來判斷樣本類別。

其次是為了實現演算法中複雜的矩陣計算，除了第一題和第二題的第一小問我是用C語言實現的，後面的程式都是用matlab實現的。當然自己也是一邊學習一邊使用。程式設計中遇到了很多困難，都是通過百度查詢來解決的。經過幾天的摸索，雖然不能說寫了這幾個演算法自己就多精通matlab了，但至少自己已經入門，知道怎麼用這個強大的科學計算工具來處理複雜的數學問題。

當然學習一門新的語言對學習計算機專業的學生來說並不是什麼很難的事情，其實最重要的是弄清書中的演算法實現。所以在編寫程式之前我總要花幾個小時來用心瞭解書中的理論知識。所以學習這門課也培養我鑽研知識的能力，以及培養耐心。因為很多程式總是伴隨著很多bug，所以你不得不花大量的時間來除錯程式。下面我把自己相關的程式附上作為附錄。

作業相關程式:

貝葉斯演算法:

clc;clear all;d=2;c=3;N=50;errnum = 0;N1=33;N2=40;N3=26;

p1 = N1/N;p2 =N2/N;p3=N3/N;D=load('Iris_data.txt');

data=zeros(150,d);G1=zeros(N1,d);G2=zeros(N2,d);G3=zeros(N3,d);

hG1=zeros(50-N1,d);hG2=zeros(50-N2,d);hG3=zeros(50-N3,d);

 data(:,1)=D(:,2); data(:,2)=D(:,4);

for i=1:1:N1 G1(i,:)=data(i,:);end

for i=1:1:N-N1 hG1(i,:)=data(i+N1,:);end

for i=1:1:N2  G2(i,:)=data(N+i,:);end

for i=1:1:N-N2  hG2(i,:)=data(N+N2+i,:);end

for i=1:1:N3  G3(i,:)=data(2*N+i,:);end

for i=1:1:N-N3  hG3(i,:)=data(2*N+N3+i,:);end

miu1=mean(G1,1)miu2=mean(G2,1)miu3=mean(G3,1)

sigma1=zeros(d,d);sigma2=zeros(d,d);sigma3=zeros(d,d);

for i=1:1:N1  sigma1=sigma1+(G1(i,:)-miu1)'*(G1(i,:)-miu1);  end

for i=1:1:N2   sigma2=sigma2+(G2(i,:)-miu2)'*(G2(i,:)-miu2);end

for i=1:1:N3   sigma3=sigma3+(G3(i,:)-miu3)'*(G3(i,:)-miu3);end

sigma1=sigma1/N1;sigma2=sigma2/N2;sigma3=sigma3/N3;

R=zeros(150,3);syms x1 x2 positive;

temp = [x1    x2];format short;

d1=log(p1)-1/2*log(det(sigma1))-1/2*(temp*inv(sigma1)*temp')+miu1*inv(sigma1)*temp'-1/2*miu1*inv(sigma1)*miu1';

d2=log(p2)-1/2*log(det(sigma2))-1/2*(temp*inv(sigma2)*temp')+miu2*inv(sigma2)*temp'-1/2*miu2*inv(sigma2)*miu2';

d3=log(p3)-1/2*log(det(sigma3))-1/2*(temp*inv(sigma3)*temp')+miu3*inv(sigma3)*temp'-1/2*miu3*inv(sigma3)*miu3';

fun1 [email protected](temp)(log(p1)-1/2*(temp-miu1)*inv(sigma1)*((temp-miu1)')-1/2*log(det(sigma1)));

fun2 [email protected](temp)(log(p2)-1/2*(temp-miu2)*inv(sigma2)*((temp-miu2)')-1/2*log(det(sigma2)));

fun3 [email protected](temp)(log(p3)-1/2*(temp-miu3)*inv(sigma3)*((temp-miu3)')-1/2*log(det(sigma3))');

sd1 = vpa(d1,4);sd2 =vpa(d2,4);sd3 = vpa(d3,4);

format shortaxis([0,13,1.8,14]);

for i=1:1:N-N1  hold on   x1= data(i+N1,1);  x2= data(i+N1,2);

   if fun1([x1,x2])>fun2([x1,x2])&&fun1([x1,x2])>fun3([x1,x2])

        scatter(x2,x1,40,'r','v','filled');else

scatter(x2,x1,40,[0.3,0,0],'v','filled'); errnum = errnum+1;end  end

for i=1:1:N-N2  hold on  x1= data(i+N2+N,1);  x2 = data(i+N2+N,2);

    if fun2([x1,x2])>fun1([x1,x2])&&fun2([x1,x2])>fun1([x1,x2])

        scatter(x2,x1,40,'b','^','filled'); else

       scatter(x2,x1,40,[0,0,0.3],'v','filled');

       errnum = errnum+1; end

end

for i=1:1:N-N3 hold on   x1 = data(i+N3+2*N,1);  x2 = data(i+N3+2*N,2);

    if fun3([x1,x2])>fun1([x1,x2])&&fun3([x1,x2])>fun2([x1,x2])

        scatter(x2,x1,40,'g','o','filled');else  scatter(x2,x1,40,[0,0.3,0],'o','filled');

       errnum = errnum+1;end

end errnum

 fisher演算法:

 num1= 26;num2 = 40; num3=33;sum=num1+num2+num3;

p1 =num1/sum;p2=num3/sum;p3=num3/sum; temp = [0,0];

for i=1:length(data1)  temp = temp+data1(i,:); end

meandata1 =temp/length(data1); meandata1'temp = [0,0];

 for i=1:length(data2)  temp = temp+data2(i,:);end

 meandata2 = temp/length(data2); meandata2'

 temp = [0,0];for i=1:length(data3)   temp = temp+data3(i,:);end

 meandata3 = temp/length(data3); meandata3'

 S1 = zeros(2);S2 = zeros(2);S3 = zeros(2);

for i=1:length(data1)

   S1 =S1+(meandata1-data1(i,:))'*(meandata1-data1(i,:));end

for i=1:length(data2) S2 =S2+(meandata2-data2(i,:))'*(meandata2-data2(i,:));end

for i=1:length(data3)  S3 =S3+(meandata3-data3(i,:))'*(meandata3-data3(i,:));end

Sw = S1+S2+S3; meandata0= meandata1*p1+meandata2*p2+meandata3*p3;

Sb =(meandata1-meandata0)'*(meandata1-meandata0)*p1+...(meandata2-meandata0)'*(meandata2-meandata0)*p2+(meandata3-meandata0)'*(meandata3-meandata0)*p3;W= Sw\Sb; ins1 = W\S1*W;ins2 =W\S2*W;ins3 = W\S3*W;

、