模式識別---貝葉斯分類

1.數據庫

• Dataset1.txt：328個同學的身高、體重、性別數據（78個女生250個男生）
• Dataset2.txt：124個同學的數據（40女、84男）
• Dataset3.txt：90個同學的數據（16女，74男）

數據集：提取碼：e8ph 

2.需要完成的工作

（1）以dataset1為訓練數據庫，假設身高與體重滿足高斯分布，進行高斯分布的參數估計，並進行基於最小錯誤率的貝葉斯分類，分別考慮男女的先驗概率，0.5-0.5；0.6-0.4；0.7-0.3，0.8-0.2，並以dataset2和dataset3為測試數據庫分析分類性能，並探討先驗概率對分類性能的影響。

（2）以dataset1為訓練數據庫，進行基於pazen窗方法的概率密度估計，並進行0.5-0.5先驗概率條件下的最小錯誤率並進行基於最小錯誤率的貝葉斯分類並以dataset2和dataset3為測試數據庫分析分類性能。

（3）基於線性分類器進行性別分類，同樣以dataset1作為訓練樣本，dataset2和dataset3作為測試樣本。

 3.方法與思路

　　由前兩問可以看出，其實判別的函數就是下圖，就是給定一個待測向量X，它是類別Wi的概率。

　　先驗概率p（wi）和類條件概率p（xlwi）都是未知的。根據僅有的樣本數據進行分類時，一種可行的辦法是我們需要先對先驗概率和類條件概率進行估計，然后再套用貝葉斯分類器。

先驗概率的估計較簡單：

• 每個樣本所屬的自然狀態都是已知的（有監督學習）；
• 依靠經驗；
• 用訓練樣本中各類出現的頻率估計。

類條件概率的估計（非常難），原因包括：

• 概率密度函數包含了一個隨機變量的全部信息；
• 概率密度函數可以是滿足p（x）≥0；fp（x）dx=1的任何函數；
• 在很多情況下，已有的訓練樣本數總是太少；
• 當用於表示特征的向量x的維數較大時，就會產生嚴重的計算復雜度問題（算法的執行時間，系統資源開銷等）。

　　這里討論的主要是對類條件概率的估計。估計的方法分為兩大類：
參數估計（parametric）：

　　參數估計法，樣本所屬的類別和類條件概率密度函數形式已知，而表征概率密度函數的某些參數是未知的。要求由已知類別的訓練數據樣本集，對概率密度的某些參數進行統計估計。如：最大似然估計，貝葉斯估計等。
　　思路：x->0->p（xlwi），將概率密度估計轉化為參數估計對於本次作業來說，我們假設類條件概率密度的函數形式是高斯型，第一想法的是下列式子：

　　然而，上面常見的高斯概率密度函數只是針對一維的參數X，對於大多數情況，輸入參數會是多維的，這里我們針對二元變量的概率密度函數求解為：

　　說明：上邊的參數p是多變量間的相關系數，設定值應小於1。

非參數估計（non-parametric）：

　　已知樣本所屬類別，但未知概率密度函數的形式，要求不用模型，而只利用訓練數據本身對概率密度做估計。如：parzen窗函數辦法，k近鄰估計等。

　　思路：x->p（xlwi），直接用樣本來估計概率密度函數，能處理任意的概率分布，而不必假設概率密度函數的形式已知。

parzen窗的原理如下：

第三問，參考這篇文章。

4.程序實現

第一問：兩種方法，基於sklearn與不基於sklearn，如下：

# -*- coding: utf-8 -*-

#Import Library of Gaussian Naive Bayes model
from sklearn.metrics import accuracy_score 
from sklearn.naive_bayes import GaussianNB


def getData(filename):
    fr1=open(filename)
    
    arrayOLines1=fr1.readlines() #讀取文件
    
    #特征矩陣
    X=[]
    y=[]
    #男同學
    for line in arrayOLines1:
        line=line.strip() #strip()去掉首尾空格
        listFromLine=line.split() #按空白字符分割成列表
       
        #將m，f變成1,0
        for i in range(len(listFromLine)):
            if listFromLine[2] == 'm':
               listFromLine[2] = 1
            elif listFromLine[2] == 'f':
                listFromLine[2] = 0
            else:
                listFromLine[2]=listFromLine[2]
        #string變float
        for i in range(len(listFromLine)):
            listFromLine[i]=float(listFromLine[i])

        #X是指特征向量；y是標簽
        X.append([listFromLine[0],listFromLine[1]])
        y.append(listFromLine[2])

    return X,y

#訓練數據
features_train,labels_train = getData('dataset1.txt')
#print(features_train,labels_train)
#創建一個高斯型分布的分類器
#model = GaussianNB(priors=[0.2,0.8])
model = GaussianNB()

# Train the model using the training sets 
model.fit(features_train, labels_train)

# 用訓練好的分類器去預測測試集的標簽值 
features_test,labels_test = getData('dataset2.txt')
pred =model.predict(features_test) 

# 計算並返回在測試集上的准確率 
y_pred =pred 
y_true =labels_test
accuracy_score(y_true, y_pred) 
  
print(accuracy_score(y_true, y_pred,normalize=False))

# -*- coding: utf-8 -*-

#-*-encoding:utf-8-*-
import numpy
import math
 
def importdata(filename = 'dataset1.txt') :
    '''
    導入訓練集
    '''
    f = open(filename,'r')
    dataset = []
    #arr = []
    for item in f :
        vars = item.split()
        dataset.append([float(vars[0]), float(vars[1]), vars[2].upper()])
    return dataset 
 
def getParameters(dataset) :
    '''
    從訓練集分別獲取不同類別下的期望、方差、標准差、類別的先驗概率以及變量間相關系數
    '''
    class1 = []
    class2 = []
    class_sum = []
    for item in dataset :
        class_sum.append([item[0],item[1]])
        if item[-1] == 'F' :
            class1.append([item[0],item[1]])
        if item[-1] == 'M' :
            class2.append([item[0],item[1]])
    class1 = numpy.array(class1)
    class2 = numpy.array(class2)
    class_total = numpy.array(class_sum)
    mean1 = numpy.mean(class1,axis=0)
    variance1 = numpy.var(class1,axis=0)
    stand_deviation1 = numpy.std(class1,axis=0)
    mean2 = numpy.mean(class2,axis=0)
    variance2 = numpy.var(class2,axis=0)
    stand_deviation2 = numpy.std(class2,axis=0)
    class_total = (len(class1) + len(class2)) * 1.0
     
    mean = numpy.mean(class_sum, axis=0)
    stand_deviation = numpy.std(class_sum, axis=0)
     
    new_arr = [ ((item[0] - mean[0]) * (item[1] - mean[1]) / stand_deviation[0] / stand_deviation[1])  for item in dataset]
    coefficient = numpy.mean(new_arr)
 
    return (mean1,mean2),(variance1,variance2),(stand_deviation1, stand_deviation2),(len(class1)/class_total,len(class2)/class_total),coefficient
      
 
def GaussianFunc(mean, variance, stand_deviation, coefficient) :
    '''
    根據指定參數（期望、方差、標准差、多元向量間的相關性）生成高斯函數
    多元變量的高斯函數
    '''
    def func(X) :
        X = [X[0] - mean[0], X[1] - mean[1]]
        B = [[variance[0], coefficient * stand_deviation[0] * stand_deviation[1]],[coefficient * stand_deviation[0] * stand_deviation[1], variance[1]]]
        inv_B = numpy.linalg.inv(B)
        A = inv_B
        B_val = (1.0 - coefficient**2) * variance[0] * variance[1]
        tmp1 = 2*math.pi * (B_val ** 0.5)
        X = numpy.array([X])
        tmp2 = (-0.5) * numpy.dot(numpy.dot(X, A), X.T)
        res = 1.0 / tmp1 * (math.e ** tmp2)
        return res
    return func
 
 
def f(X, funcs, class_ps, index) :
    '''
    貝葉斯概率計算函數
    '''
    tmp1 = funcs[index](X) * class_ps[index]
    tmp2 = funcs[0](X) * class_ps[0] + funcs[1](X) * class_ps[1]
    return tmp1 / tmp2
 
 
def classify(X,funcs,class_ps,labels) :
    '''
    基於最小錯誤率的貝葉斯判別分類。對於二類分類問題，簡化了。
    '''
    res1 = f(X,funcs,class_ps,0)
    res2 = f(X,funcs,class_ps,1) 
    if res1 > res2 :
        return labels[0]
    else :
        return labels[1]
 
def test(dataset, funcs,class_ps,labels) :
    '''
    測試
    '''
    positive0 = 0
    positive1 = 0
    F = [item for item in dataset if item[-1] == 'F']
    len_F = len(F)
    len_M = len(dataset) - len_F
    for item in dataset :
        res = classify([item[0],item[1]], funcs, class_ps,labels)
        if res == item[-1] and res == 'F' :
            positive0 += 1
        if res == item[-1] and res == 'M' :
            positive1 += 1
    print ('F', positive0 * 1.0 / len_F)
    print ('M', positive1 * 1.0 / len_M)
 
 
if __name__ == '__main__' :
    dataset = importdata()
    (mean1,mean2),(variance1,variance2),(stand_deviation1, stand_deviation2), (class1_p, class2_p), coefficient = getParameters(dataset)
    func1 = GaussianFunc(mean1, variance1, stand_deviation1,coefficient)
    func2 = GaussianFunc(mean2, variance2, stand_deviation2,coefficient)
    #print func1([160,45])
    #print func1([170,50])
    #print func1([175,50])
    #print func1([190,20])
    funcs = []
    funcs.append(func1)
    funcs.append(func2)
    class_ps = []
    class_ps.append(class1_p)
    class_ps.append(class2_p)
 
    classs = [class_ps]
    '''
    手工指定先驗概率
    '''
    classs.append([0.5,0.5])
    classs.append([0.4,0.6])
    classs.append([0.3,0.7])
    classs.append([0.2,0.8])
 
    labels = ['F', 'M']
    for class_ps in classs :
        print( '-' * 24)
        print (class_ps)
        
        print ('-'*10,'dataset1','-'*10)
        testset0 = importdata('dataset1.txt')
        test(testset0, funcs, class_ps, labels)
        
        print ('-'*10,'dataset2','-'*10)
        testset1 = importdata('dataset2.txt')
        test(testset1, funcs, class_ps, labels)
        
        print ('-'*10,'dataset3','-'*10)
        testset2 = importdata('dataset3.txt')
        test(testset2, funcs, class_ps, labels)

第二問：

# -*- coding: utf-8 -*-

import numpy as np
 
 
# 讀取數據
def getData(filename):
    fr1=open(filename)
    
    arrayOLines1=fr1.readlines() #讀取文件
    
    #特征矩陣
    X=[]
    #y=[]
    #男同學
    for line in arrayOLines1:
        line=line.strip() #strip()去掉首尾空格
        listFromLine=line.split() #按空白字符分割成列表
       
        #將m，f變成1,0
        for i in range(len(listFromLine)):
            if listFromLine[2] == 'm':
               listFromLine[2] = 1
            elif listFromLine[2] == 'f':
                listFromLine[2] = 0
            else:
                listFromLine[2]=listFromLine[2]
        #string變float
        for i in range(len(listFromLine)):
            listFromLine[i]=float(listFromLine[i])

        #X是指特征向量；y是標簽
        #X.append([listFromLine[0],listFromLine[1]])
        #y.append(listFromLine[2])
        X.append(listFromLine)
    return X
 
 
def get_phi(x, xi, h):
    x = np.mat(x)
    xi = np.mat(xi)
    phi = np.exp(-(x - xi) * (x - xi).T / (2 * h * h))
    return phi
 
 
def get_px(x, xi, h):
    phi = 0
    n = len(xi)
    for i in range(n):
        # print("xi[i]", xi[i])
        phi += get_phi(x, xi[i], h)
    px = phi * 3 / (4 * np.pi * n * np.power(h, 3))
    return px
 
 
def parzen(h, test, data):
    positive0 = 0
    positive1 = 0
    
    px = [0, 0]
    print("h =", h)
    
    for j in range(len(test)):
        #print("x =", test[j])
        for i in range(len(px)):
            xi = [x[:3] for x in filter(lambda x: x[2] == i , data)]
        
            # print("xi", xi)
            # print("len xi", len(xi))
            px[i] = get_px(test[j], xi, h)
           # print("px", i, "=", px[i])
            
        if px[0] > px[1] :
            pedict0 = 0
            if pedict0  == test[j][2]:
                positive0 += 1
            #print("belong to w1")
            
        if px[0] < px[1]:
            #print("belong to w2")
            pedict1 = 1
            if pedict1  == test[j][2]:
                positive1 += 1
                
    print ('F', positive0 * 1.0 / 40 )
    print ('M', positive1 * 1.0 / 84)
 
def main():
    data = getData('dataset1.txt')
    # print(np.mat(data))
    test = getData('dataset2.txt')
    h1 = 1
    h2 = 0.1
 
    parzen(h1, test, data)
    parzen(h2, test, data)
 
 
if __name__ == '__main__':
    main()

第三問：

# -*- coding: utf-8 -*-

import matplotlib.pyplot as plt
import numpy as np
 

def importdata(filename) :
    '''
    導入訓練集
    '''
    f = open(filename,'r')
    dataset = []
    #arr = []
    for item in f :
        vars = item.split()
        dataset.append([float(vars[0]), float(vars[1]), vars[2].upper()])
    return dataset 

def getparameters(dataset) :
    '''
    從訓練集分別獲取不同類別下的期望、方差、標准差、類別的先驗概率以及變量間相關系數
    '''
    class1 = []
    class2 = []
    class_sum = []
    for item in dataset :
        class_sum.append([item[0],item[1]])
        if item[-1] == 'F' :
            class1.append([item[0],item[1]])
        if item[-1] == 'M' :
            class2.append([item[0],item[1]])
    class1 = np.array(class1)
    class2 = np.array(class2)
    class_total = np.array(class_sum)
    return class1,class2,class_total

#fisher算法實現
#為下一個函數fisher中的一個語句做定義
def cal_cov_and_avg(samples):
    """
    給定一個類別的數據，計算協方差矩陣和平均向量
    :param samples: 
    :return: 
    """
    u1 = np.mean(samples, axis=0)
    cov_m = np.zeros((samples.shape[1], samples.shape[1]))
    for s in samples:
        t = s - u1
        cov_m += t * t.reshape(2, 1)
    return cov_m, u1

def fisher(c_1, c_2):
    """
    fisher算法實現(請參考上面推導出來的公式，那個才是精華部分)
    :param c_1: 
    :param c_2: 
    :return: 
    """
    cov_1, u1 = cal_cov_and_avg(c_1)#調用上一個函數
    cov_2, u2 = cal_cov_and_avg(c_2)#調用上一個函數
    s_w = cov_1 + cov_2
    u, s, v = np.linalg.svd(s_w)  # 奇異值分解
    s_w_inv = np.dot(np.dot(v.T, np.linalg.inv(np.diag(s))), u.T)
    return np.dot(s_w_inv, u1 - u2)

def judge(sample, w, c_1, c_2):
    """
    true 屬於1
    false 屬於2
    :param sample:
    :param w:
    :param center_1:
    :param center_2:
    :return:
    """
    u1 = np.mean(c_1, axis=0)
    u2 = np.mean(c_2, axis=0)
    center_1 = np.dot(w.T, u1)
    center_2 = np.dot(w.T, u2)
    pos = np.dot(w.T, sample)
    return abs(pos - center_1) < abs(pos - center_2)





if __name__ == '__main__' :
    dataset = importdata('dataset1.txt')
    c_1,c_2,c_total = getparameters(dataset)
    #print('c_1',c_1)
    #print('c_2',c_2)
    (h1,w) =np.shape(c_1)
    #print(h1)
    (h2,w) =np.shape(c_2)
    #print(h2)
    
    
    
    w = fisher(c_1, c_2)# 調用函數，得到參數w
     
    for i in range(0,h1):
        out1 = judge(c_1[i], w, c_1, c_2)   # 判斷所屬的類別
        #print(out1)
    for i in range(0,h2):
        out2 = judge(c_2[i], w, c_1, c_2)   # 判斷所屬的類別
        #print(out2)




    plt.title("scatter diagram")
    plt.scatter(c_1[:, 0], c_1[:, 1], color="red")
    plt.scatter(c_2[:, 0], c_2[:, 1], color="black")
    line_x = np.arange(min(np.min(c_1[:, 0]), np.min(c_2[:, 0])),
                   max(np.max(c_1[:, 0]), np.max(c_2[:, 0])),
                   step=1)

    line_y = -(w[1] * line_x) / w[0]
    plt.plot(line_x, line_y)
    plt.xlim(xmax=200,xmin=150)
    plt.ylim(ymax=200,ymin=-100)
    plt.xlabel("x")
    plt.ylabel("y")