bayes.c

良好的代码实现

#include "stdio.h"
#include "math.h"
#define N 2   //模式样本的维数
#define N1 4  //第一个模式类中的模式样本的个数
#define N2 4  //第二个模式类中的模式样本的个数
/*变量说明: 
a[][]:第一个模式类的样本输入数据
b[][]:第二个模式类的样本输入数据
omega1:第一个模式类的先验概率
omega2:第二个模式类的先验概率
m1[1][N]: 第一个模式类的均值向量
m2[1][N]:第二个模式类的均值向量
c1[]:第一个模式类的协方差矩阵
c2[]:第二个模式类的协方差矩阵
c[][]:用来存放均值向量与其转置之积的N*N矩阵 (暂时使用一下,过后可重新赋值,另做它用)  
delta1:第一个模式类协方差矩阵的行列式值
delta2:第二个模式类协方差矩阵的行列式值
const1:第一个判别函数的常数项
const2:第二个判别函数的常数项*/
void main()
{int i,j,l;           //循环变量
 float  const1,const2,omega1,omega2,delta1,delta2,result1[1][1]={0},result2[1][1]={0};
 float a[N1][N],b[N2][N],c1[N][N],c2[N][N],c[N][N],m1[1][N],m2[1][N],am[N],an[N],aw[N];
 float  cinver(float *a,int n); // 矩阵求逆及求行列式值的函数
 void   swap(float *,float *);
 void matrixmultiply(float *,float *,float *,int ,int ,int );  //矩阵相乘函数
 FILE *in,*out;               //定义输入、输出文件指针
 if((out=fopen("output.txt","w"))==NULL)
	 exit(0);
 if((in=fopen("input.txt","r"))==NULL)
	 exit(0);
 fscanf(in,"%f",&omega1);     //读入第一个模式类的先验概率
 fscanf(in,"%f",&omega2);     //读入第二个模式类的先验概率
 for(i=0;i<N1;i++)            //读入第一个模式类的样本数据
	 for(j=0;j<N;j++)
	 {fscanf(in,"%f,",&a[i][j]);}
 for(i=0;i<N2;i++)            //读入第二个模式类的样本数据
    for(j=0;j<N;j++)
	 {fscanf(in,"%f,",&b[i][j]);}
for(i=0;i<N;i++)
{    {m1[0][i]=0.0;  m2[0][i]=0.0;
      am[i]=0.0;  an[i]=0.0;  aw[i]=0.0;
     }
     for(j=0;j<N;j++)
	     {c1[i][j]=0.0;
         c2[i][j]=0.0;
         c[i][j]=0.0;
	     }
}   
 //求均值向量m1[],m2[]
 for(i=0;i<N;i++)
    for(j=0;j<N1;j++)
	{m1[0][i]+=a[j][i]/N1;}
 for(i=0;i<N;i++)
    for(j=0;j<N2;j++)
    {m2[0][i]+=b[j][i]/N2;}
  //求协方差矩阵c1[][]
   for(l=0;l<N1;l++)
      for(i=0;i<N;i++)
        for(j=0;j<N;j++)
		 {c[i][j]=0;
		  c1[i][j]+=(a[l][i]*a[l][j])/N1;
          c[i][j]=m1[0][i]*m1[0][j];
         }
 for(i=0;i<N;i++)
   for(j=0;j<N;j++)
   {c1[i][j]=c1[i][j]-c[i][j];}
 //求协方差矩阵c2[][]
for(l=0;l<N2;l++)
      for(i=0;i<N;i++)
	     for(j=0;j<N;j++)
		 {c[i][j]=0;
		  c2[i][j]+=(b[l][i]*b[l][j])/N2;
          c[i][j]=m2[0][i]*m2[0][j];
		 }
for(i=0;i<N;i++)
  for(j=0;j<N;j++)
  {c2[i][j]=c2[i][j]-c[i][j];}
  //求协方差矩阵c1[][]的逆矩阵c1[][]和行列式的值delta1,协方差矩阵c2[][]的逆矩阵c2[][]和行列式的值delta2
   delta1=cinver(c1,N);
   delta2=cinver(c2,N);
   if((delta1<=0)||(delta2<=0))
   exit(0);
  //求第一个模式类判别函数的常数项const1
  for(i=0;i<N;i++)
     for(j=i;j<N;j++)
		 {c[i][j]=0.0;}
   matrixmultiply(m1,c1,c,1,N,N);
   matrixmultiply(c,m1,result1,1,N,1);
   const1=log(omega1)-0.5*log(delta1)-0.5*result1[0][0];
  //求第一个模式类判别函数的一次项系数
  for(i=0;i<N;i++)
      {am[i]=0;aw[i]=0;
        {for(j=0;j<N;j++)
		 am[i]+=c1[i][j]*m1[0][j];
	     }
      }
   for(i=0;i<N;i++)
        {for(j=0;j<N;j++)
           {aw[i]+=c1[j][i]*m1[0][j];}
            aw[i]=(am[i]+aw[i])/2;
        }
  //求第一个模式类判别函数的二次项系数
    for(i=0;i<N;i++)
		 for(j=i;j<N;j++)
		   if(i!=j)
		   {c1[i][j]+=(c1[j][i]);
		    c1[i][j]=-0.5*c1[i][j];
		   }
 //------------------------------输出第一个模式类别的判别函数-----------------------------
		fprintf(out,"第一个模式类的判别函数为:");
		    //二次项系数
	   	   for(i=0;i<N;i++)
		       for(j=i;j<N;j++)
			   {fprintf(out,"(%3.1fx%d*x%d)+",c1[i][j],i+1,j+1);}
	    	 //一次项系数
             for(i=0;i<N;i++)
			    {fprintf(out,"(%3.1fx%d)+",aw[i],i+1);}
    		 //常数项系数
		   fprintf(out,"(%3.1f)",const1);
//求第二个模式类别的常数项const2
   for(i=0;i<N;i++)
		 for(j=i;j<N;j++)
		 {c[i][j]=0;}
  matrixmultiply(m2,c2,c,1,N,N);
  matrixmultiply(c,m2,result2,1,N,1);
  const2=log(omega2)-0.5*log(delta2)-0.5*result2[0][0];
//求第二个模式类判别函数的一次项系数
 	for(i=0;i<N;i++)
	{am[i]=0; an[i]=0;
	   for(j=0;j<N;j++)
	   {am[i]+=c2[i][j]*m2[0][j];}
    }
   for(i=0;i<N;i++)
	     for(j=0;j<N;j++)
		    {an[i]+=c2[j][i]*m2[0][j];}
   for(i=0;i<N;i++)
       {an[i]=(am[i]+an[i])/2;}
   //求第二个模式类别的二次项的系数
	 for(i=0;i<N;i++)
		 for(j=i;j<N;j++)
		   if(i!=j)
		   {c2[i][j]+=(c2[j][i]);
		    c2[i][j]=-0.5*c2[i][j];
		   }
/*------------------- 输出第二个模式类的判别函数----------------------------------*/
		    //二次项系数
	   	   fprintf(out,"\n\n\n\n第二个模式类的判别函数为:");
		   for(i=0;i<N;i++)
		       for(j=i;j<N;j++)
			   {fprintf(out,"(%3.1fx%d*x%d)+",c2[i][j],i+1,j+1);}
            //一次项系数
             for(i=0;i<N;i++)
			 {fprintf(out,"(%3.1fx%d)+",an[i],i+1);}
		    //常数项
		    fprintf(out,"(%3.1f)",const2);
//----------------------- 输出判别界面函数-----------------------------------------------
		  	//二次项系数
	   	    fprintf(out,"\n\n\n\n判别界面为:");
	    	for(i=0;i<N;i++)
		       for(j=i;j<N;j++)
			   {fprintf(out,"(%3.1fx%d*x%d)+",(c1[i][j]-c2[i][j])/2,i+1,j+1);}
            //一次项系数
             for(i=0;i<N;i++)
			 {fprintf(out,"(%3.1fx%d)+",(aw[i]-an[i]),i+1);}
		    //常数项
		    fprintf(out,"(%3.1f)=0",(const1-const2));
fclose(in);
fclose(out);
}
void swap(float *a,float *b)
{float temp;
temp=*a;*a=*b;*b=temp;
}  	  
//  函数cinver用来求矩阵的逆矩阵和行列式的值,其中,a是矩阵名,a最初是原矩阵,经过该函数处理后,
//  a用来存放处理的结果,即逆矩阵.n是矩阵的行(列)的值.函数inverse的返回值为矩阵的行列值.
float  cinver(float *a,int n)
{	float d,delta=1.0;
    int *is,*js,i,j,k,l,m,flag=1;       
	is=malloc(n*sizeof(int));
	js=malloc(n*sizeof(int));
 for(k=0;k<n;k++)
 {d=0.0;
    for(i=k;i<n;i++)
	   {for(j=k;j<n;j++)
			  if(fabs(*(a+i*n+j))>d)
			  {d=*(a+i*n+j);
			   is[k]=i;
			   js[k]=j;
			  }
	    }
    if(is[k]!=k)
	{l=is[k];
	 flag*=-1;
	for(i=0;i<n;i++)
	  {swap(a+l*n+i,a+k*n+i);}
	}
if(js[k]!=k)
	{m=js[k];
	 flag*=-1;
		for(i=0;i<n;i++)
		{swap(a+i*n+m,a+i*n+k);}
	}
	delta=delta*(*(a+k*n+k));
	*(a+k*n+k)=1/(*(a+k*n+k));
	for(j=0;j<n;j++)
	{if(j!=k)
	  *(a+k*n+j)=*(a+k*n+j)*(*(a+k*n+k));
	}
	for(i=0;i<n;i++)
	 {if(i!=k)
	    {for(j=0;j<n;j++)
		  if(j!=k)
		  {*(a+i*n+j)=*(a+i*n+j)-(*(a+i*n+k))*(*(a+k*n+j));}
	    }
	 }
	for(i=0;i<n;i++)
	 {if(i!=k)
		   *(a+i*n+k)=(-1)*(*(a+i*n+k))*(*(a+k*n+k));
	}
}
	for(k=n-1;k>=0;k--)
	   { for(j=0;j<n;j++)
		 {i=js[k];
		  swap(a+k*n+j,a+i*n+j);
		 }
		for(i=0;i<n;i++)
		{j=is[k];
		 swap(a+i*n+k,a+i*n+j);
        }
	}
   free(is);
   free(js);
   return(delta*flag);	 
}
//matrixmultiply函数用来求解两个矩阵的乘积.A是第一个矩阵名,B是第二个矩阵名,C是结果矩阵名.
//A是a*b维矩阵, B是b*c维矩阵,C是a*c维矩阵;
void matrixmultiply(float *A,float *B,float *C,int a,int b,int c)
{int m,n,j;
  for(m=0;m<a;m++)
     for(n=0;n<c;n++)
       for(j=0;j<b;j++)
	   {*(C+m*c+n)+=(*(A+m*b+j))*(*(B+j*c+n));}
}