bayes.cpp

贝叶斯算法实现

#include<iostream.h>
#include<math.h>
#include"bayes.h"
#include"matrix.h"
#include<stdlib.h>


//计算均值向量U
Matrix getU(const Matrix &X)
{
	int d = X.getCol();
	int n = X.getRow();
	Matrix U(d, 1);        //U=(u1,u2,...,ud)^T，
	//用 d x 1 的矩阵表示
	for(int k=0; k<n; k++)    //k 代表第 k 个样本
		for(int j=0; j<d; j++)   //j 代表样本的第 j 个特征值
			U.set(j, 0, U.get(j,0)+X.get(k,j));
		//这里要注意对矩阵元素的访问下标是从 0 开始的
		return 1.0/n * U; //U 是所有样本的均值，故加和后除以样本总数
}

//计算协方差矩阵E
Matrix getE(const Matrix &X, const Matrix &U)
{
	int d = X.getCol(); //X 的列数代表特征空间维数
	int n = X.getRow(); //X 的每一行代表一个样本
	Matrix E(d, d);         //E 是 d x d 维矩阵，d 为每个样本特征空间数
	//故此处可以用 X 列数来初始化它
	double *ar = new double[d];   //此数组用于从矩阵中提取行，创建新矩阵
	//从而可以进行下面的矩阵运算
	for(int k=0; k<n; k++){     //k 代表第 k 个样本
		for(int j=0; j<d; j++)
			ar[j] = X.get(k, j);    //将 X 的一行元素值放入数组 ar 中
		
		Matrix Xk(d, 1, ar);     //利用 ar 创建矩阵
		E = E + (Xk - U) * trv(Xk - U); //(X-U) * (X-U)^T
	}
	delete[] ar; //释放动态分配的数组空间
	//E 是 n 个矩阵( (Xk - U)与其转置的乘积 )的算术平均
	return 1.0/n * E;
}

//贝叶斯判别函数
int bayesFun(const Matrix U[], const Matrix E[],
			 const double Pw[],const Matrix &X, int c)
{
	double testPw = 0;
	for(int w=0; w<c; w++){
		if(Pw[w] <= 0){
			cout << "Wrong Pw[i]\n"; exit(0); }
		testPw += Pw[w];
		if(det(E[w]) == 0){
			cout << "Check input, E[" << w << "]=0!\n";
			exit(0); }
	}
	if(testPw != 1){
		cout << "Wrong Pw[i]\n"; exit(0); }
	
	//对于贝叶斯判别函数中的 -d/2 * ln(2*PI)项因与 i 无关故可去掉
	//另外(X-U)(E^-1)(X-U)^T 的计算结果虽然是一个数，但在这里
	//其表示为一个一行一列的矩阵，因此需要对其求行列式值后才能
	//赋给左边的 double maxg
	double maxg = -1.0/2 * det( trv(X-U[0]) * inv(E[0]) * (X-U[0]) )
		-1.0/2 * log(det(E[0])) + log(Pw[0]);
	int code = 0; //用来返回类别代号
	
	//求使得贝叶斯判别函数取得最大值类
	for(int i=1; i<c; i++){
		double g =   -1.0/2 * det( trv(X-U[i]) * inv(E[i]) * (X-U[i]) )
			-1.0/2 * log(det(E[i])) + log(Pw[i]);
		if(g > maxg){
			maxg = g;
			code = i;
		}
	}//for i
	
	return code;
}//bayesFun()

//贝叶斯分类器
void bayesDepart(const Matrix X[], const double Pw[],
				 const Matrix &x, int c, char* name[])
{
	Matrix* U = new Matrix[c];
	Matrix* E = new Matrix[c];
	
	//先求出每个类的均值向量和协方差矩阵
	for(int i=0; i<c; i++){
		U[i] = getU(X[i]);
		E[i] = getE(X[i], U[i]);
	}
	//调用贝叶斯判别函数，判断X所属类别，输出类别名称
	cout << name[bayesFun(U,E,Pw,x,c)];
}