fisher.cpp

一个简单的fisher分类器代码,并含有实验用的两个数据集合.

#include "stdio.h"
#include "math.h"
#include "stdlib.h"

#define trainingnum  100
#define testingnum 100
#define feadim 2
#define N1 50
#define N2 50

double s1[feadim][feadim],s2[feadim][feadim];//两类样本类内离散度
double m1[feadim][1],m2[feadim][1];   //两类样本的均值
double insw[feadim][feadim];//存放sw的逆矩阵
double w[feadim][1];//存放解结果
double _m1,_m2,yo;  //存放~m1,~m2
double y[testingnum];//一维空间样本

int right,wrong;//存放识别的结果

struct sample
{
	char serial[4];  //序号
	float feature[feadim][1];//存储身高与体重
	int trueclass;   //真实类别
	int classifiedclass; //被识别的分类
};

//打印矩阵
void PrintOut(double M[feadim][feadim])
{
 int i,j;
 for(i = 0; i<feadim;  i++) //行
 {
  for(j=0; j<feadim; j++)
   printf("%8.2f",M[i][j]);
   printf("\n");
 }
}
//采用部分主元法的高斯消去法求方阵A的逆矩阵B
//
//A 方阵  (IN)
//n 方阵的阶数 (IN)
//B 方阵  (OUT)
//返回true 说明正确计算出逆矩阵
//返回false说明矩阵A不存在逆矩阵

bool gaussj(double A[feadim][feadim], double B[feadim][feadim])
{
 int i,j,k;
 double lMax,temp;

 //临时矩阵存放A
 double tt[feadim][feadim];
 for(i=0;i<feadim;i++)
 {
  for(j=0;j<feadim;j++)
  {
   tt[i][j] = A[i][j];
  }
 }

 //初始化B为单位阵
 for(i=0;i<feadim;i++)
 {
  for(j=0;j<feadim;j++)
  {
   if(i!=j)B[i][j] = 0;
   else B[i][j] = 1;
  }
 }
 for(i=0;i<feadim;i++)
 {
  //寻找主元
  lMax = tt[i][i];
  k = i;
  for(j=i+1;j<feadim;j++) //扫描从i+1到n的各行
  {
   if( fabs(tt[j][i]) > fabs(lMax))
   {
    lMax = tt[j][i];
    k = j;
   }
  }
  //如果主元所在行不是第i行，进行行交换
  if(k!=i)
  {
   for(j=0;j<feadim;j++)
   {
    temp = tt[i][j] ;
    tt[i][j] = tt[k][j];
    tt[k][j] = temp;
    //B伴随计算
    temp = B[i][j];
    B[i][j] = B[k][j];
    B[k][j] = temp;
   }
  }
  //判断主元是否是0，如果是，则矩阵A不是满秩矩阵，不存在逆矩阵
  if(tt[i][i] == 0) return false;
  //消去A的第i列除去i行以外的各行元素
  temp = tt[i][i];
  for(j=0;j<feadim;j++)
  {
   tt[i][j] = tt[i][j] / temp; //主对角线上元素变成1
   B[i][j] = B[i][j] / temp; //伴随计算
  }

  for(j=0;j<feadim;j++) //行 0 -> n
  {
   if(j!=i)  //不是第i行
   {
    temp = tt[j][i];
    for(k=0;k<feadim;k++) // j行元素 - i行元素* j行i列元素
    {
     tt[j][k] = tt[j][k] - tt[i][k] * temp;
     B[j][k] = B[j][k] - B[i][k] * temp;
    }
   }
  }
 }
 return true;
}
//计算两个矩阵的乘积
// AB -> C
void Multi(double A[feadim][feadim], double B[feadim][1], double C[feadim][1])
{
 int i,j,k;
 for(i=0;i<feadim;i++)
 {
  for(j=0;j<1;j++)
  {
   C[i][j] = 0;
   for(k=0;k<feadim;k++)
   {
    C[i][j] += A[i][k] * B[k][j];
   }
  }
 }
}
//主程序
void main()
{   sample x[trainingnum];
    sample X[testingnum];
    int i,j,k;
	int sflag;
    double insw[feadim][feadim],sw[feadim][feadim];//存放sw的逆矩阵及sw
	double pm[feadim][1];//存放m1-m2
	double inw[1][feadim];//存放w的转置
	char snum[4];
	float sheight,sweight;
	double temp1[feadim][1],temp2[feadim][1],temp11[1][feadim],temp22[1][feadim];//在求s1,s2时存放临时数据
	//读入训练样本
	FILE *fp;
	if((fp=fopen("StudentsData1.txt","r"))==NULL)
   {
     printf("Can't open file");
     exit(0);
    }
  for(j=0; j<trainingnum; j++)   
   {
	 fscanf(fp,"%s%d%f%f",&snum,&sflag,&sheight,&sweight);
		  x[j].feature[0][0]=sheight;
		  x[j].feature[1][0]=sweight;
		  x[j].trueclass=sflag;
		  for(k=0;k<4;k++)
		  x[j].serial[k]=snum[k];
		  
    }
  
  fclose(fp); 
 /*for(i=0;i<10;i++)  //测试输出
  {
	  printf("%s%3d%8.2f%8.2f\n",x[i].serial,x[i].trueclass,x[i].feature[0][0],x[i].feature[1][0]);
  }
  */
  for(i=0;i<N1;i++)     //计算m1,m2
  {  for(j=0;j<feadim;j++)
	  m1[j][0]+=x[i].feature[j][0]/N1;
  }
  for(i=N1;i<2*N1;i++)
  {
	 for(j=0;j<feadim;j++)
	  m2[j][0]+=x[i].feature[j][0]/N2;
  }
  printf("m1,m2的结果如下：\n");
  printf("----------------- \n");
   for(j=0;j<feadim;j++)
	   printf("%5.2f%8.2f\n",m1[j][0],m2[j][0]); //输出m1,m2计算结果

 for(i=0;i<N1;i++)     //计算s1,s2
 {  for(j=0;j<feadim;j++)
	  temp1[j][0]=x[i].feature[j][0]-m1[j][0];
	  for(j=0;j<feadim;j++)  //转置temp1
	     temp11[0][j]=temp1[j][0];
	  for(k=0;k<feadim;k++)//求s1
		  {
			  for(j=0;j<feadim;j++)
                s1[k][j]+=temp1[k][0]*temp11[0][j];
		  }
	  
  }
  for(i=N1;i<2*N1;i++)
  {
	  for(j=0;j<feadim;j++)
	  temp2[j][0]=x[i].feature[j][0]-m2[j][0];
	  for(j=0;j<feadim;j++)  //转置temp2
	     temp22[0][j]=temp2[j][0];
	    for(k=0;k<feadim;k++)//求s2
		  {
			  for(j=0;j<feadim;j++)
                s2[k][j]+=temp2[k][0]*temp22[0][j];
		  }
	  
  }
  /* for(j=0;j<2;j++)
	  {  
		  temp11[0][j]=temp1[j][0];
		  temp22[0][j]=temp2[j][0];
	  }*/
  printf("\n");
  printf("s1的结果如下：\n");
  printf("----------------- \n");
  PrintOut(s1);
  
  printf("\n");
  printf("s2的结果如下：\n");
  printf("----------------- \n");
  PrintOut(s2);

  //求sw
  for(i=0;i<feadim;i++)
	  for(j=0;j<feadim;j++)
		  sw[i][j]=s1[i][j]+s2[i][j];
  printf("\n");
  printf("sw：\n");
  printf("----------------- \n");
  PrintOut(sw);
 if(!gaussj(sw,insw))//求sw的逆矩阵
 {
  printf("没有逆矩阵！\n");
 }
 else
 {
	 printf("\n的逆矩阵是:\n");
     for(i=0;i<feadim;i++)
	 {
      for(j=0;j<feadim;j++)
	  printf("%11.7f",insw[i][j]);
      printf("\n");
	 }
 }
 for(i=0;i<feadim;i++)//计算m1-m2
    pm[i][0]=m1[i][0]-m2[i][0];
 Multi(insw,pm,w);//计算w=逆sw*(m1-m2)
 printf("\n解w的结果是:\n");
 for(i=0;i<feadim;i++)
	 printf("%11.7f\n",w[i][0]);
 for(i=0;i<feadim;i++)  //转置w
	     inw[0][i]=w[i][0];
 for(i=0;i<feadim;i++) //计算一维空间的样本均值~m1,~m2
 {
	 _m1+=inw[0][i]*m1[i][0];
	 _m2+=inw[0][i]*m2[i][0];
 }
 yo=(_m1+_m2)/2;//计算阈值yo

 //读入测试样本
   FILE *fr;
	if((fr=fopen("StudentsData2.txt","r"))==NULL)
   {
     printf("Can't open file");
     exit(0);
    }
  for(j=0; j<testingnum; j++)   //读入数据
   {
	 fscanf(fr,"%s%d%f%f",&snum,&sflag,&sheight,&sweight);
		  X[j].feature[0][0]=sheight;
		  X[j].feature[1][0]=sweight;
		  X[j].trueclass=sflag;
		  for(k=0;k<4;k++)
		  X[j].serial[k]=snum[k];
		  X[j].classifiedclass=-1;//初始化被识别的类别
    }
  
  fclose(fr); 
  
 //进行判别
  for(i=0;i<testingnum;i++)
	  for(j=0;j<feadim;j++)
	  {
		  y[i]+=inw[0][j]*X[i].feature[j][0];
          if(y[i]<yo)
			  X[i].classifiedclass=-1;
		  else
			  X[i].classifiedclass=1;
	  }
  for(i=0;i<testingnum;i++)
	  if(X[i].trueclass==X[i].classifiedclass)
		  right++;
	  else
		  wrong++;
 printf("\n测试结果如下：\n");
 printf("----------------- \n");
 printf("正确分类的个数为:%3d\n",right);
 printf("错误分类的个数为:%3d\n",wrong);
}