贝叶斯决策理论对CIFAR-10数据图像分类

1：数据集介绍

CIFAR-10 是一个用于普世物体识别的数据集，分为airplane、automobile、bird、cat、deer、
dog、frog、horse、ship、truck共10 类。
共60000张32*32大小的彩色RGB图像，分为10类，50000张用于训练，10000张用于测试。
训练集分为5个训练batches，测试集分为1个测试batch。每个batch有10类，每类大概有1000张图片。
本次实验中选用大约10000张图片作为训练集，这些图片类别为1或2或3
选用3000张图片作为测试集

2：算法说明

实验主要分为两步
第一步对数据集使用PCA进行降维，每张图片由1024x3=3072个像素点构成，降维以后每张图片保留20个特征；
第二步基于贝叶斯决策理论进行判决，基于本实验样本分布，选用高斯分布的朴素贝叶斯分类算法。

2.1 PCA（主成分分析）

PCA(Principal Component Analysis)，即主成分分析方法，是一种使用最广泛的数据降维算法。PCA的主要思想是将n维特征映射到k维上，这k维是全新的正交特征也被称为主成分，是在原有n维特征的基础上重新构造出来的k维特征。PCA的工作就是从原始的空间中顺序地找一组相互正交的坐标轴，新的坐标轴的选择与数据本身是密切相关的。其中，第一个新坐标轴选择是原始数据中方差最大的方向，第二个新坐标轴选取是与第一个坐标轴正交的平面中使得方差最大的，第三个轴是与第1,2个轴正交的平面中方差最大的。依次类推，可以得到n个这样的坐标轴。通过这种方式获得的新的坐标轴，我们发现，大部分方差都包含在前面k个坐标轴中，后面的坐标轴所含的方差几乎为0。于是，我们可以忽略余下的坐标轴，只保留前面k个含有绝大部分方差的坐标轴。事实上，这相当于只保留包含绝大部分方差的维度特征，而忽略包含方差几乎为0的特征维度，实现对数据特征的降维处理。
问题在于如何得到这些包含最大差异性的主成分方向
事实上，通过计算数据矩阵的协方差矩阵，然后得到协方差矩阵的特征值特征向量，选择特征值最大(即方差最大)的k个特征所对应的特征向量组成的矩阵。这样就可以将数据矩阵转换到新的空间当中，实现数据特征的降维。
由于得到协方差矩阵的特征值特征向量有两种方法：特征值分解协方差矩阵、奇异值分解协方差矩阵，所以PCA算法有两种实现方法：基于特征值分解协方差矩阵实现PCA算法、基于SVD分解协方差矩阵实现PCA算法

这就是PCA的优化目标
PCA算法描述如下：
输入：样本集D，低维空间维数d^l
过程：
1：对所有样本进行中心化；
2：计算样本的协方差矩阵；
3：对协方差矩阵做特征值分解；
4：取最大的d^l个特征值对应的特征向量 $w_1,w_2,……w_(d^l)$w1​,w2​,……w(​dl)
输出：投影矩阵 $W^*=(w_1,w_2,……w_(d^l ))$W∗=(w1​,w2​,……w(​dl))

2.2 贝叶斯决策

贝叶斯决策理论是主观贝叶斯派归纳理论的重要组成部分。 贝叶斯决策就是在不完全情报下，对部分未知的状态用主观概率估计，然后用贝叶斯公式对发生概率进行修正，最后再利用期望值和修正概率做出最优决策
基本思想：
已知类条件概率密度参数表达式和先验概率
利用贝叶斯公式转换成后验概率
根据后验概率大小进行决策分类
理论分析：
（1）如果我们已知被分类类别概率分布的形式和已经标记类别的训练样本集合，那我们就需要从训练样本集合中来估计概率分布的参数。在现实世界中有时会出现这种情况。（如已知为正态分布了，根据标记好类别的样本来估计参数，常见的是极大似然率和贝叶斯参数估计方法）
（2）如果我们不知道任何有关被分类类别概率分布的知识，已知已经标记类别的训练样本集合和判别式函数的形式，那我们就需要从训练样本集合中来估计判别式函数的参数。在现实世界中有时会出现这种情况。（如已知判别式函数为线性或二次的，那么就要根据训练样本来估计判别式的参数，常见的是线性判别式和神经网络）
（3）如果我们既不知道任何有关被分类类别概率分布的知识，也不知道判别式函数的形式，只有已经标记类别的训练样本集合。那我们就需要从训练样本集合中来估计概率分布函数的参数。在现实世界中经常出现这种情况。（如首先要估计是什么分布，再估计参数。常见的是非参数估计）
（4）只有没有标记类别的训练样本集合。这是经常发生的情形。我们需要对训练样本集合进行聚类，从而估计它们概率分布的参数。（这是无监督的学习）
（5）如果我们已知被分类类别的概率分布，那么，我们不需要训练样本集合，利用贝叶斯决策理论就可以设计最优分类器。但是，在现实世界中从没有出现过这种情况。这里是贝叶斯决策理论常用的地方。
问题：假设我们将根据特征矢量x 提供的证据来分类某个物体，那么我们进行分类的标准是什么？decide wj， if（p(wj|x)>p(wi|x)）(i不等于j)应用贝叶斯展开后可以得到p(x|wj)p(wj)>p(x|wi)p(wi)即或然率p(x|wj)/p(x|wi)>p(wi)/p(wj)，决策规则就是似然率测试规则。
结论：
对于任何给定问题，可以通过似然率测试决策规则得到最小的错误概率。这个错误概率称为贝叶斯错误率，且是所有分类器中可以得到的最好结果。最小化错误概率的决策规则就是最大化后验概率判据。
决策判据：
贝叶斯决策理论方法是统计模式识别中的一个基本方法。贝叶斯决策判据既考虑了各类参考总体出现的概率大小，又考虑了因误判造成的损失大小，判别能力强。贝叶斯方法更适用于下列场合：
(1) 样本(子样)的数量(容量)不充分大，因而大子样统计理论不适宜的场合。
(2) 试验具有继承性，反映在统计学上就是要具有在试验之前已有先验信息的场合。用这种方法进行分类时要求两点：
第一，要决策分类的参考总体的类别数是一定的。例如两类参考总体(正常状态Dl和异常状态D2)，或L类参考总体D1，D2，…，DL(如良好、满意、可以、不满意、不允许、……)。
第二，各类参考总体的概率分布是已知的，即每一类参考总体出现的先验概率P(Di)以及各类概率密度函数P(x/Di)是已知的。显然，0≤P(Di)≤1，(i=l，2，…，L)，∑P(Di)=1。
贝叶斯公式告诉我们：
P(A|B)=P(AB)/P(B)=P(B|A)*P(A)/P(B)
从公式中可知，如果要计算B条件下A发生的概率，只需要计算出后面等式的三个部分，B事件的概率（P(B)），是B的先验概率、A属于某类的概率（P(A)），是A的先验概率、以及已知A的某个分类下，事件B的概率（P(B|A)），是后验概率。
如果要确定某个样本归属于哪一类，则需要计算出归属不同类的概率，再从中挑选出最大的概率
我们把上面的贝叶斯公式写出这样，也许你能更好的理解：
MAX(P(Ai|B))=MAX(P(B|Ai)*P(Ai)/P(B))
而这个公式告诉我们，需要计算最大的后验概率，只需要计算出分子的最大值即可，而不同水平的概率P©非常容易获得，故难点就在于P(X|C)的概率计算。而问题的解决，正是聪明之处，即贝叶斯假设变量X间是条件独立的，故而P(X|C)的概率就可以计算为：
P(B|Ai) =P(B1/Ai)*P(B2/Ai)P(B3/Ai)…*P(Bn/Ai)
具体到本实验中
$p(c│x)=(P(c)P(x|c))/(P(x))=(P(c))/(P(x)) ∏_(i=1)^d▒〖P(x_i |c)〗$p(c│x)=(P(c)P(x∣c))/(P(x))=(P(c))/(P(x))∏(​i=1)d▒〖P(xi​∣c)〗
其中，d为属性数目， $x_i$xi​为x在第i个属性上的取值，c为类别
对连续属性可以考虑概率密度函数，本实验假定数据分布为高斯， $p(x_i│c)~N(u_(c,i),〖σ_(c,i)〗^2)$p(xi​│c) N(u(​c,i),〖σ(​c,i)〗2),其中 $u_(c,i)$u(​c,i) 和 $〖σ_(c,i)〗^2$〖σ(​c,i)〗2分别是第c类样本在第i个属性上取值的均值和方差。

3：实验分析

考虑到python语言的简单易用及在机器学习领域包含的众多的库函数可供调用，本实验使用python进行编程，本实验除使用sklearn库内置的PCA函数来对数据集降维外，其余处理判决等都自己编程实现。
用PCA保留20维情况下，分类结果如图：
可以看出，此时分类准确率为0.642
当然，准确率随PCA保留维度数量的不同而有些差异，维度数量由程序中的featureCount参数决定，欲实验不同维度下分类效果的不同，只要改变此参数即可

代码：

# -*- coding: UTF-8 -*-
import numpy as np
import math
import scipy.io as sio
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
'''
该数据集共有60000张彩色图像，这些图像是32*32，分为10个类，
每类6000张图。这里面有50000张用于训练，构成了5个训练批，
每一批10000张图；另外10000用于测试，单独构成一批。测试批的数据里，
取自10类中的每一类，每一类随机取1000张。抽剩下的就随机排列组成了训练批。
注意一个训练批中的各类图像并不一定数量相同，总的来看训练批，每一类都有5000张图
使用loadmat读取训练集，返回一个字典，有5个键值对，其中Data是9968*3072的图片数据numpy.ndarray格式
每一行是一张图片32*32*3=3072，一共9968张
'''
def loadData(file):
    #file='G:/lecture of grade one/pattern recognition/data/train_data.mat'
    trainImg=sio.loadmat(file)
    print trainImg["Data"][1,:]
    return trainImg

def dimenReduc(data,featureCount):
    pca = PCA(n_components=featureCount)
    newData = pca.fit_transform(data)
    #print(pca.explained_variance_ratio_)
    return newData

def extractClass(label,data,c):
    #提取出第c类图片放到listC
    listC=[]
    for i in range(len(label)):
        if label[i]==c:
            listC.append(list(data[i,:]))
    #pc=len(listC)/len(label)
    return listC

def calcEC(trainImg,featureCount):
    matrixE=np.empty((5,featureCount))#matrixE[i][j]表示第i类样本在第j个属性上的均值
    matrixC=np.empty((5,featureCount))
    pc=[]
    #newtrainData=dimenReduc(trainImg['Data'],featureCount)
    
    pca = PCA(n_components=featureCount)
    newtrainData=pca.fit_transform(trainImg['Data'])
    
    for i in range(5):
        listC=extractClass(trainImg['Label'],newtrainData,i+1)
        pc.append(float(len(listC))/float(len(trainImg['Label'])))
        for j in range(featureCount):
            list1=[]
            for k in range(len(listC)):
                list1.append(listC[k][j])
                matrixE[i,j]=np.mean(list1)
                matrixC[i,j]=np.cov(list1)
            list1=[]
    return matrixE,matrixC,pc,newtrainData,pca         
        
def calcPostrior(matrixE,matrixC,pc,testData,featureCount):
    #给定一张图片，返回p(c|x)的最大值的索引，即后验概率p(c|x)最大值对应的类
    post=[]
    for i in range(3): 
        p=1
        for j in range(featureCount):
            p=p*math.exp(-(testData[j]-matrixE[i][j])**2/(2*matrixC[i][j]))/math.sqrt(2*math.pi*matrixC[i][j])
        post.append(pc[i]*p)
    
    return post.index(max(post))+1
        

def extractTestImg(testImg):
    #从大的数据集中抽取测试集，就是挑出类别为1，2,3的图片
    newTestImg={'Label':[],'Data':[]}
    Label=[]

    index=[]
    for i in range(len(testImg["Label"])):
        if testImg["Label"][i]==1 or testImg["Label"][i]==2 or testImg["Label"][i]==3:
            index.append(i)
            Label.append(testImg["Label"][i])
    Data=np.empty((len(index),3072)) 
    Label=np.empty((len(index),1))       
    for i in range(len(index)):
        Data[i,:]=testImg["Data"][index[i],:]
        Label[i]=testImg["Label"][index[i]]
    newTestImg['Label']=Label
    newTestImg['Data']=Data
    return newTestImg

if __name__ == '__main__':
   featureCount=20
   file='G:/lecture of grade one/pattern recognition/data/train_data.mat'
   trainImg = loadData(file)
   matrixE,matrixC,pc,newTrainData,pca=calcEC(trainImg,featureCount)
   file='G:/lecture of grade one/pattern recognition/data/test_data.mat'
   testImg = extractTestImg(loadData(file))
   newTestData=pca.transform(testImg['Data'])
   #newTestData=dimenReduc(testImg['Data'],featureCount)
   preLabel=[]
   for i in range(len(testImg['Label'])):
       preLabel.append(calcPostrior(matrixE,matrixC,pc,newTestData[i,:],featureCount))
   errorCount=0.0
   for i in range(len(preLabel)):
       if preLabel[i] != testImg['Label'][i]:
           errorCount+=1
   print('Errorcount are %f'%errorCount)
   print('Total test count are %f'%len(testImg['Label']))
   print('Test accuracy is %f'%(1.0-errorCount/float(len(testImg['Label']))))