fisher线性判别
来源:互联网 发布:日本原单淘宝店 编辑:程序博客网 时间:2024/05/22 03:05
Fisher线性判别(Fisher Linear Discrimination,FLD),也称线性判别式分析(Linear Discriminant Analysis, LDA)。FLD是基于样本类别进行整体特征提取的有效方法。它在使用PCA方法进行降维的基础上考虑到训练样本的类间信息。FLD的基本原理就是找到一个最合适的投影轴,使各类样本在该轴上投影之间的距离尽可能远,而每一类内的样本的投影尽可能紧凑,从而使分类效果达到最佳,即在最大化类间距离的同时最小化类内距离。FLD方法在进行图像整体特征提取方面有着广泛的应用。
在应用统计方法解决模式识别问题时,经常会遇到所谓的“维数灾难”的问题,在低维空间里适用的方法在高维空间里可能完全不适用。因此压缩特征空间的维数有时是很重要的。Fisher方法实际上涉及维数压缩的问题。如果把多维特征空间的点投影到一条直线上,就能把特征空间压缩成一维,这个在数学上是很容易办到的。但是,在高维空间里很容易分开的样品,把它们投影到任意一根直线上,有可能不同类别的样品就混在一起,无法区分,如图1(a)所示投影到xl或x2轴无法区分。若把直线绕原点转动一下,就有可能找到一个方向,样品投影到这个方向的直线上,各类样品就能很好地分开,如图1(b)所示。因此直线方向的选择很重要。一般地,总能够找到一个最好的方向,使样品投影到这个方向的直线上很容易分开。如何找到这个最好的直线方向以及如何实现向最好方向投影的变换,这正是Fisher算法要解决的基本问题,这个投影变换恰是我们所寻求的解向量w*。
图1 Fisher线性判别示意图
样品训练集以及待测样品的特征总数目为n。为了找到最佳投影方向,需要计算出各类样品均值,样品类内离散度矩阵Si和总类间离散度矩阵Sw,样品类间离散度矩阵Sb,根据Fisher准则,找到最佳投影向量,将训练集内所有样品进行投影,投影到一维Y空间,由于Y空间是一维的,则需要求出Y空间的划分边界点,找到边界点后,就可以对待测样品进行一维Y空间的投影,判断它的投影点与分界点的关系,将其归类。具体方法如下。
-
-
-
-
-
-
- int Classification::Fisher_2Classes(int Class0, int Class1)
- {
- double Xmeans[2][25];
- double S[2][25][25];
- double Sw[25][25];
- double Sw_[25][25];
- double W[25];
- double difXmeans[25];
- double X[25];
- double m0,m1;
- double y0;
- int i,j,k;
-
- for(i=0;i<2;i++)
- for(j=0;j<25;j++)
- Xmeans[i][j]=0;
- int num0,num1;
-
- double mode0[200][25],mode1[200][25];
-
- num0=40;
- num1=40;
- for(i=0;i<num0;i++)
- {
- for(j=0;j<25;j++)
- {
- Xmeans[0][j]+=pattern[Class0].feature[i][j];
- mode0[i][j]=pattern[Class0].feature[i][j];
- }
- }
-
- for(i=0;i<num1;i++)
- {
- for(j=0;j<25;j++)
- {
- Xmeans[1][j]+=pattern[Class1].feature[i][j];
- mode1[i][j]=pattern[Class1].feature[i][j];
- }
- }
-
- for(i=0;i<25;i++)
- {
- Xmeans[0][i]/=(double)num0;
- Xmeans[1][i]/=(double)num1;
- }
-
- for(i=0;i<25;i++)
- for(j=0;j<25;j++)
- {
- double s0=0.0,s1=0.0;
- for(k=0;k<num0;k++)
- s0=s0+(mode0[k][i]-Xmeans[0][i])*(mode0[k][j]-Xmeans[0][j]);
- s0=s0/(double)(num0-1);
- S[0][i][j]=s0;
- for(k=0;k<num1;k++)
- s1=s1+(mode1[k][i]-Xmeans[1][i])*(mode1[k][j]-Xmeans[1][j]);
- s1=s1/(double)(num1-1);
- S[1][i][j]=s1;
- }
-
- for(i=0;i<25;i++)
- for(j=0;j<25;j++)
- {
- Sw[i][j]=S[0][i][j]+S[1][i][j];
- }
-
- for(i=0;i<25;i++)
- for(j=0;j<25;j++)
- Sw_[i][j]=Sw[i][j];
- double(*p)[25]=Sw_;
- brinv(*p,25);
-
- for(i=0;i<25;i++)
- difXmeans[i]=Xmeans[0][i]-Xmeans[1][i];
- for(i=0;i<25;i++)
- W[i]=0.0;
- brmul(Sw_,difXmeans,25,W);
-
-
- m0=0.0;
- m1=0.0;
- for(i=0;i<num0;i++)
- {
- m0+=brmul(W,mode0[i],25);
- }
- for(i=0;i<num1;i++)
- {
- m1+=brmul(W,mode1[i],25);
- }
- m0/=(double)num0;
- m1/=(double)num1;
- y0=(num0*m0+num1*m1)/(num0+num1);
-
-
- for(i=0;i<25;i++)
- X[i]=testsample[i];
- double y;
- y=brmul(W,X,25);
- if (y>=y0)
- return Class0;
- else
- return Class1;
- }
-
-
-
-
-
-
- int Classification::Fisher()
- {
- int i,j,number,maxval,num[10];
- for(i=0;i<10;i++)
- num[i]=0;
- for(i=0;i<10;i++)
- for(j=0;j<i;j++)
- num[Fisher_2Classes(i,j)]++;
- maxval=num[0];
- number=0;
- for(i=1;i<10;i++)
- {
- if(num[i]>maxval)
- {
- maxval=num[i];
- number=i;
- }
- }
- return number;
- }
-
-
-
-
-
-
-
-
-
- void brmul(double a[],double b[][25],int n,double c[])
- {
- for(int i=0;i<n;i++)
- {
- for(int j=0;j<n;j++)
- c[i]+=a[j]*b[j][i];
- }
- return;
- }
-
-
-
-
-
-
- void brinv(double a[],int n)
- {
- int *is,*js,i,j,k,l,u,v;
- double d,p;
- is=new int[n];
- js=new int[n];
- for (k=0; k<=n-1; k++)
- {
- d=0.0;
- for (i=k; i<=n-1; i++)
- for (j=k; j<=n-1; j++)
- {
- l=i*n+j; p=fabs(a[l]);
- if (p>d)
- {
- d=p; is[k]=i; js[k]=j;
- }
- }
- if (d+1.0==1.0)
- {
- free(is); free(js); printf("err**not inv\n");
- return;
- }
- if (is[k]!=k)
- for (j=0; j<=n-1; j++)
- {
- u=k*n+j; v=is[k]*n+j;
- p=a[u]; a[u]=a[v]; a[v]=p;
- }
- if (js[k]!=k)
- for (i=0; i<=n-1; i++)
- {
- u=i*n+k; v=i*n+js[k];
- p=a[u]; a[u]=a[v]; a[v]=p;
- }
- l=k*n+k;
- a[l]=1.0/a[l];
- for (j=0; j<=n-1; j++)
- if (j!=k)
- {
- u=k*n+j; a[u]=a[u]*a[l];
- }
- for (i=0; i<=n-1; i++)
- if (i!=k)
- for (j=0; j<=n-1; j++)
- if (j!=k)
- {
- u=i*n+j;
- a[u]=a[u]-a[i*n+k]*a[k*n+j];
- }
- for (i=0; i<=n-1; i++)
- if (i!=k)
- {
- u=i*n+k; a[u]=-a[u]*a[l];
- }
- }
- for (k=n-1; k>=0; k--)
- {
- if (js[k]!=k)
- for (j=0; j<=n-1; j++)
- {
- u=k*n+j; v=js[k]*n+j;
- p=a[u]; a[u]=a[v]; a[v]=p;
- }
- if (is[k]!=k)
- for (i=0; i<=n-1; i++)
- {
- u=i*n+k; v=i*n+is[k];
- p=a[u]; a[u]=a[v]; a[v]=p;
- }
- }
- delete is;
- delete js;
- }