fisher线性判别+KNN处理iris数据集

fisher线性判别+KNN处理iris数据集

---

问题概述

iris数据集3类，4维，150个数据
sonar数据2类，60维，208个样本
样本分类方法：留一法

思路概述

fisher概述

类均值
$${\mu }_i={\frac{1}{N}}_i\sum _{x\in {\omega }_i}x$$
总均值
$$\mu =\frac{1}{N}\sum _{x\in \omega i}{N}_i{\mu }_i$$
各类内离散度矩阵
$$S_i=\sum _{x_j\in {\Omega }_i}^n(x_j-{\mu }_i)(x_j-{\mu }_i{)}^T$$
总类内离散度矩阵
$$S_{\omega }=\sum _{i=1}^c{S}_{\omega i}$$
二分类时的样本类间离散度
$$S_b=({\mu }_1-{\mu }_2)({\mu }_1-{\mu }_2{)}^T$$
如果是多分类的话样本类间离散度矩阵应该是这样
$$S_b=\sum _{i=1}^c{N}_i({\mu }_i-\mu )({\mu }_i-\mu {)}^T$$
fisher准则函数(广义Rayleigh商)为
$$ma{x}_{\omega }{J}_F(\omega )=\frac{{\omega }^T{S}_b\omega }{{\omega }^T{S}_{\omega }\omega }$$
为求满足上式的${\omega }^{\ast }$，可以引入拉格朗日乘子求解，假设分母为非零常数，即令
$${\omega }^T{S}_{\omega }\omega =c\ne 0$$

则可定义拉格朗日函数为
$$L(\omega ,\lambda )={\omega }^T{S}_b\omega -\lambda ({\omega }^T{S}_{\omega }\omega -c)$$

上式对$\omega$求偏导数得

$$\frac{\partial L(\omega ,\lambda }{\partial \omega }=S_b\omega -\lambda S_{\omega }\omega$$

令上式等于0化简得

$$S_{\omega }^{-1}{S}_b{\omega }^{\ast }=\lambda {\omega }^{\ast }$$

对于iris数据集这种三分类问题来说,我们要求出最大的n个特征值$\lambda $,求出其对应的$\omega_n$让数据分别映射在这n个方向上从而达到降维的目的。
如果是像sonar这种二分类问题我们可以直接进一步化简为
$${\omega }^{\ast }=S_{\omega }^{-1}({\mu }_1-{\mu }_2)$$
从而得到我们想要的方向然后使其降到一维

iris数据集

利用fisher线性判别将数据降维
使用k近邻等算法分类验证
用留一法统计出识别率

sonar

由于数据维数较高，所以先进行数据归一化处理，然后进行fisher判别，并用留一法统计出识别率

iris 代码分析

降维环节

数据读取部分

利用fgetc函数逐个读取字符，iris数据集较简单。
每一行中第一个和第三个字符为第一维属性的个位和小数部分，
每一行中第五个和第七个字符为第一维属性的个位和小数部分，
每一行中第九个和第十一个字符为第一维属性的个位和小数部分，
每一行中第十三个和第十五个字符为第一维属性的个位和小数部分。
第一类Iris-setosa数据每一行数据的换行符在第二十七个符号位。
第二类Iris-versicolor数据每一行数据的换行符在第三十一个符号位。
第三类Iris-virginica数据每一行数据的换行符在第三十个符号位。
基于以上几点特征将数据集中所有信息读入整形数组data中，data第一维为数据序号，data第二维的五个数分别是数据集的四个属性和一个类别标识位(第一类为0，第二类为1，第三类为2)

void fld::load_pram()
{
    FILE *f;
    	if ((f = fopen("/home/jojo/jason_cpp/iris.data", "rb+")) == NULL)
	{
		perror("fail to read");
	}
     if(f != NULL)
    {
        ch = 'a';
        while(ch != EOF)
        {
            putchar(ch);
            ch = fgetc(f);
            if(ch!='\n')ch_1[tag_a++]=ch;
            else
            {
                data[tag_b++][0]=(ch_1[0]-'0')*10+(ch_1[2]-'0');
                data[tag_b-1][1]=(ch_1[4]-'0')*10+(ch_1[6]-'0');
                data[tag_b-1][2]=(ch_1[8]-'0')*10+(ch_1[10]-'0');
                data[tag_b-1][3]=(ch_1[12]-'0')*10+(ch_1[14]-'0');
                if(tag_a==27)data[tag_b-1][4]=0;
                else if(tag_a==31)data[tag_b-1][4]=1;
                else data[tag_b-1][4]=2;
                tag_a=0;
            }
        }
    }
    else
        printf("fail to open! \n");
    fclose(f);
}

将整形数据data导入double数组data_db中

void fld::deal_data()
{
    for(tag_b=0;tag_b<150;tag_b++)
    {
        for(tag_a=0;tag_a<4;tag_a++)
        {
            //写的有点罗嗦（读的时候是两位数的整形，这里转换为双精度浮点数）
            data_db[tag_b][tag_a]=(data[tag_b][tag_a]%10)*0.1+data[tag_b][tag_a]/10;
        }
        data_db[tag_b][4]=data[tag_b][4];
    }
}

留一法载入训练集和测试集

函数参数只有一个，即将这一标号的数据作为测试数据，其他所有样本都作为训练样本
将data_db的第m个样本传到data_test中其他放入data_train中

void fld::load_datat(int m)
{
    key = m;
    for(tag_b=0,tag_a=0;tag_b<150;tag_b++,tag_a++)
    {
        if(m==tag_b)
        {
            for(tag_c=0;tag_c<5;tag_c++)
                data_test[0][tag_c]=data_db[tag_b][tag_c];
            tag_b++;
        }
        for(tag_c=0;tag_c<5;tag_c++)
            data_train[tag_a][tag_c]=data_db[tag_b][tag_c];
    }
}

计算训练集均值

首先将各类的总数求出并存入数组num_type中，
然后将每一维，每一类属性求均值放入二维数组u中
定义 double u[3][4],u_all[4]={0,0,0,0};
u对应三维样本的四个元素
u_all是所有样本的四个均值

void fld::load_u_num()
{
    for(tag_b=0;tag_b<3;tag_b++)
    num_type[tag_b]=0;
    for(tag_a=0;tag_a<3;tag_a++)
    for(tag_b=0;tag_b<4;tag_b++)
    u[tag_a][tag_b]=0;

    for(tag_a=0;tag_a<num[0];tag_a++)
    num_type[int(data_train[tag_a][4])]++;
    for(tag_a=0;tag_a<num[0];tag_a++)
    {
        for(tag_b=0;tag_b<4;tag_b++)
        {
            //for accuracy that may be change
            u[int(data_train[tag_a][4])][tag_b]+=(data_train[tag_a][tag_b]/num_type[int(data_train[tag_a][4])]);
        }
    }
    for(tag_a=0;tag_a<num[0];tag_a++)
    for(tag_b=0;tag_b<4;tag_b++)
    u_all[tag_b]+=data_train[tag_a][tag_b];
    for(tag_b=0;tag_b<4;tag_b++)
    u_all[tag_b]/=num[0];

}

计算$S_{\omega }$和$S_b$矩阵

其中sw[3]存放的是sw[0],sw[1],sw[2]的和，即三类类内散度的和
num_type是每一类的数量

void fld::load_sw_sb()
{
    for(tag_a=0;tag_a<4;tag_a++)
    for(tag_c=0;tag_c<4;tag_c++)
    for(tag_b=0;tag_b<4;tag_b++)
        sw[tag_a][tag_b][tag_c]=0;
    for(tag_a=0;tag_a<num[0];tag_a++)
    for(tag_c=0;tag_c<4;tag_c++)
    for(tag_b=0;tag_b<4;tag_b++)
        sw[int(data_train[tag_a][4])][tag_b][tag_c]+=((data_train[tag_a][tag_b]-u[int(data_train[tag_a][4])][tag_b])*(data_train[tag_a][tag_c])-u[int(data_train[tag_a][4])][tag_c]);
    for(tag_a=0;tag_a<3;tag_a++)
    for(tag_c=0;tag_c<4;tag_c++)
    for(tag_b=0;tag_b<4;tag_b++)
        sw[3][tag_c][tag_b]+=sw[tag_a][tag_c][tag_b];

    for(tag_c=0;tag_c<4;tag_c++)
    for(tag_b=0;tag_b<4;tag_b++)
        sb[tag_b][tag_c]=0;
    for(tag_a=0;tag_a<3;tag_a++)
    for(tag_c=0;tag_c<4;tag_c++)
    for(tag_b=0;tag_b<4;tag_b++)
        sb[tag_b][tag_c]+=((u[tag_a][tag_b]-u_all[tag_b])*(u[tag_a][tag_c]-u_all[tag_c])*num_type[tag_a]);
}

计算矩阵$S_{\omega }^{-1}{S}_b$

结果放入la数组中，brinv函数用来计算矩阵的逆

void fld::value()
{
    for(tag_a=0;tag_a<4;tag_a++)
    for(tag_b=0;tag_b<4;tag_b++)
        sw_[tag_a][tag_b]=sw[3][tag_a][tag_b];
    double (*p)[4]=sw_;
    brinv(*p,4);
    for(tag_a=0;tag_a<4;tag_a++)
    for(tag_b=0;tag_b<4;tag_b++)
        la[tag_a][tag_b]=0;
    for(tag_a=0;tag_a<4;tag_a++)
    for(tag_b=0;tag_b<4;tag_b++)
    for(tag_c=0;tag_c<4;tag_c++)
        la[tag_a][tag_b]+=sw_[tag_a][tag_c]*sb[tag_c][tag_b];
        
}

计算$S_{\omega }^{-1}{S}_b$特征值和特征向量

调用了eigen库求特征值的函数es.pseudoEigenvalueMatrix()
和求特征向量的函数es.pseudoEigenvectors()
然后找出特征值最大的两个特征向量存入w数组

void fld::ev_deal()
{

    Matrix4d ma;
    for(tag_a=0;tag_a<4;tag_a++)
    for(tag_b=0;tag_b<4;tag_b++)
    ma(tag_a,tag_b)=la[tag_a][tag_b];
    EigenSolver<Matrix4d> es(ma);
	MatrixXd D = es.pseudoEigenvalueMatrix();
	MatrixXd V = es.pseudoEigenvectors();
    
    double tag[4]={-9990,-9990,-9990,-9990};
    //第一个是最大特征值，第二个是特征值对应序号，第三个是第二大特征值，第四个 是序号
    for(tag_a=0;tag_a<4;tag_a++)
    if(D(tag_a,tag_a)>tag[0])
    {
        tag[0]=D(tag_a,tag_a);
        tag[1]=tag_a;
    }
    D(tag[1],tag[1])=-9990;
    for(tag_a=0;tag_a<4;tag_a++)
        if(D(tag_a,tag_a)>tag[2])
    {
        tag[2]=D(tag_a,tag_a);
        tag[3]=tag_a;
    }
    for(tag_a=0;tag_a<4;tag_a++)
        w[0][tag_a]=V(tag[1],tag_a);
    for(tag_a=0;tag_a<4;tag_a++)
        w[1][tag_a]=V(tag[3],tag_a);
}

将数据进行投影降维

将data_db*w得到的结果放入data_new
tag存放的是test测试数据的号，在这里把降维后的测试数据放入data_newtest
降维后的测试数据放入data_newtrain

void fld::last_deal()
{
    for(tag_a=0;tag_a<150;tag_a++)
    for(tag_b=0;tag_b<2;tag_b++)
        data_new[tag_a][tag_b]=0;
    for(tag_a=0;tag_a<150;tag_a++)
        data_new[tag_a][2]=data[tag_a][4];
    for(tag_a=0;tag_a<150;tag_a++)
    for(tag_b=0;tag_b<4;tag_b++)
    {
        data_new[tag_a][0]+=data_db[tag_a][tag_b]*w[0][tag_b];
        data_new[tag_a][1]+=data_db[tag_a][tag_b]*w[1][tag_b];    
    }
    for(tag_b=0,tag_a=0;tag_b<150;tag_b++,tag_a++)
    {

        if(key==tag_b)
        {
            for(tag_c=0;tag_c<3;tag_c++)
                data_newtest[0][tag_c]=data_new[tag_b][tag_c];
            tag_a--;
            continue;
        }
        for(tag_c=0;tag_c<3;tag_c++)
        data_newtrain[tag_a][tag_c]=data_new[tag_a][tag_c];
    }

}

fld对象

class fld 
{
    private:
    //test在data中的序号
    int key ;
    //tag用来循环计数，data用于存放数据集中的数据（个位，十位）
    int data[150][5],tag_a = 0,tag_b=0,tag_c=0;
    double data_train[149][5];
    double data_test[1][5];
    double w[2][4];
    double data_db[150][5];
    double u[3][4],u_all[4]={0,0,0,0};
    //sw[3] is sum of sw[i]
    double sw[4][4][4],sb[4][4],sw_[4][4];
    double la[4][4];
    public:
    double data_newtrain[149][3];
    double data_newtest[1][3];
    int num_type[3];
    double data_new[150][3];
    //load data
    void load_pram();
    //data deal
    void deal_data();
    //show data
    void show_data();
    //load data_test and data_train
    //半成品
    void load_datat(int m);
    //calculate all u(perhaps be bug);
    void load_u_num();

    void load_sw_sb();

    void value();

    void brinv(double a[],int n);

    void ev_deal();

    void last_deal();
};

分类环节

向fld_2d载入降维后的均值（用于归一化处理）

void fld_2d::load_u(fld a)
{

    for(tag_c=0;tag_c<3;tag_c++)
    for(tag_a=0;tag_a<2;tag_a++)
        u[tag_c][tag_a]=0;
    for(tag_a=0;tag_a<149;tag_a++)
    for(tag_b=0;tag_b<2;tag_b++)
        u[int(a.data_newtrain[tag_a][2])][tag_b]+=a.data_newtrain[tag_a][tag_b];
    for(tag_c=0;tag_c<3;tag_c++)
    for(tag_a=0;tag_a<2;tag_a++)
        u[tag_c][tag_a]/=a.num_type[tag_c];

}

向fld_2d载入数据并进行归一化处理使降维后的数据尽量能够在800*800的图像中显示出其分布

void fld_2d::load_data( fld a)
{
    for(tag_a=0;tag_a<149;tag_a++)
    for(tag_b=0;tag_b<2;tag_b++)
    data_train[tag_a][tag_b] = 3*(a.data_newtrain[tag_a][tag_b]-((u[1][tag_b]+u[2][tag_b]+u[0][tag_b])/3))+400;
    data_test[0]=3*(a.data_newtest[0][0]-((u[1][0]+u[2][0]+u[0][0])/3))+400;
    data_test[1]=3*(a.data_newtest[0][1]-((u[1][1]+u[2][1]+u[0][1])/3))+400;
    data_test[2]=a.data_newtest[0][2];
      for(tag_a=0;tag_a<149;tag_a++)
    data_train[tag_a][2]=a.data_newtrain[tag_a][2];
}

KNN分类

找出降维后距离test最近的k个样本，统计k个样本中哪种最多，将test也划分为这一类
然后判断knn和test原本分类是否相符合以此来验证fisher的效果

bool fld_2d::knn(int k)
{
    double data_knn[149][4];
    for(tag_a=0;tag_a<149;tag_a++)
    for(tag_b=0;tag_b<3;tag_b++)
        data_knn[tag_a][tag_b]=data_train[tag_a][tag_b];
    for(tag_a=0;tag_a<149;tag_a++)
        data_knn[tag_a][3]=pow((data_knn[tag_a][0]-data_test[0]),2)+pow((data_knn[tag_a][1]-data_test[1]),2);

    double tag[2]={99999,0},point[50][3];
    for(tag_b=0;tag_b<k;tag_b++)
    {           
        tag[0]=99999;
        for(tag_a=0;tag_a<149;tag_a++)
        {
            if(data_knn[tag_a][3]<tag[0])
            {
                tag[0]=data_knn[tag_a][3];
                tag[1]=tag_a;
            }
        }
        data_knn[int (tag[1])][3]=99999;
        point[tag_b][0]=data_knn[int (tag[1])][0];
        point[tag_b][1]=data_knn[int (tag[1])][1];
        point[tag_b][2]=data_knn[int (tag[1])][2];
    }
    int list[3]={0,0,0};
    for(tag_b=0;tag_b<k;tag_b++)
        list[int(point[tag_b][2])]++;

    Mat a=Mat(Size(800,800), CV_8UC3);
    for(tag_a=0;tag_a<800;tag_a++)
    for(tag_b=0;tag_b<800;tag_b++)
    for(tag_c=0;tag_c<3;tag_c++)
        a.at<Vec3b>(tag_a,tag_b)[tag_c]=0;

    for(tag_a=0;tag_a<700;tag_a++)
    {
        a.at<Vec3b>(tag_a,700)[0]=200;
        a.at<Vec3b>(tag_a,700)[1]=200;
        a.at<Vec3b>(tag_a,700)[2]=200;
        a.at<Vec3b>(700,tag_a)[0]=200;
        a.at<Vec3b>(700,tag_a)[1]=200;
        a.at<Vec3b>(700,tag_a)[2]=200;
    }
    for(tag_a=0;tag_a<149;tag_a++)
    {
        if(data_train[tag_a][2]==0)
        {
            for(tag_b=-1;tag_b<1;tag_b++)   
            for(tag_c=-1;tag_c<1;tag_c++)   
                a.at<Vec3b>(data_train[tag_a][0]+tag_b,data_train[tag_a][1]+tag_c)[0]=200;
                //printf("%lf %lf %lf \n",data_new[tag_a][0],data_new[tag_a][1],data_new[tag_a][2]);

        }
        else if(data_train[tag_a][2]==1)
        {
            for(tag_b=-1;tag_b<1;tag_b++)       
            for(tag_c=-1;tag_c<1;tag_c++)   
                a.at<Vec3b>(data_train[tag_a][0]+tag_b,data_train[tag_a][1]+tag_c)[1]=200;
                    //printf("%lf %lf %lf \n",data_new[tag_a][0],data_new[tag_a][1],data_new[tag_a][2]);

        }
        else if(data_train[tag_a][2]==2)
        {
            for(tag_b=-1;tag_b<1;tag_b++)   
            for(tag_c=-1;tag_c<1;tag_c++)   
                a.at<Vec3b>(data_train[tag_a][0]+tag_b,data_train[tag_a][1]+tag_c)[2]=200;
                    //printf("%lf %lf %lf \n",data_new[tag_a][0],data_new[tag_a][1],data_new[tag_a][2]);

        }

    }
            
    for(tag_b=-3;tag_b<4;tag_b++)
    for(tag_a=-3;tag_a<4;tag_a++)
    {

        a.at<Vec3b>(data_test[0]+tag_b,data_test[1]+tag_a)[0]=200;
        a.at<Vec3b>(data_test[0]+tag_b,data_test[1]+tag_a)[1]=200;
    }

    for(tag_a=0;tag_a<k;tag_a++)  
    for(tag_b=-2;tag_b<3;tag_b++)   
    for(tag_c=-2;tag_c<3;tag_c++)  
        a.at<Vec3b>(point[tag_a][0]+tag_b,point[tag_a][1]+tag_c)[int(point[tag_a][2])]=255;
    imshow("a",a);
    waitKey(0);
    if(list[int (data_test[2])]>=list[0]&&list[int (data_test[2])]>=list[1]&&list[int (data_test[2])]>=list[2])
    return true;
    else return false;
}

分类环节fld_2d类

class fld_2d
{
    private:
    double data_train[150][3],data_test[3];
    int tag_a,tag_b,tag_c,sum[3];
    double data_db[150][3];
    double sw_[2][2],la[4][4],u[3][2],s[149][4],sw[2][2]={0,0,0,0};
    public:
    void load_data( fld a);
    void load_u(fld a);
    bool test();
    bool knn(int k);
    int test_type();
};

主函数

int main()
{

    int m=150,tag;
    fld_2d b;
    fld a;
    a.load_pram();
    a.deal_data();//两位数*0.1转换浮点数
    float all=0,r=0,AA[3][2]={0,0,0,0,0,0};
    double A[3];
    
    //留一法（将序号为m的数据作为测试集）
    while(m--)
    {
        //将测试训练集分开
        a.load_datat(m);
        //输出载入数据
        //a.show_data();
        //计算均值

        a.load_u_num();
        //计算sw，sb
        a.load_sw_sb();
        //计算sw的逆sw_，然后将sb*sw_矩阵赋值给la
        a.value();
        //求la矩阵的特征值
        a.ev_deal();
        //测试，训练样本分别投影
        a.last_deal();
        //画降维后的分布图
        //计算均值
        b.load_u(a);
        //载入数据
        b.load_data(a);
        if(b.knn(19))
        {
            AA[b.test_type()][0]+=1.0;
            r++;
        }
        AA[b.test_type()][1 ]+=1.0;
        all++;

    }
    for(tag=0;tag<3;tag++)
        A[tag]=AA[tag][0]/AA[tag][1];
    printf("\nOA=%.3f\n",r/all);
    printf("AA=%.3f\n",(A[0]+A[1]+A[2])/3);
    return 0;
}

部分数据更新流程图：

部分数据降维后19NN分布情况

三种颜色（红绿蓝）分别对应三类，白色是测试数据，放大点为knn样本数据

sonar 代码分析

主函数

int main ()
{
    float all=0,r=0,AA[2][2]={0,0,0,0};
    double A[2];
    int m=208,tag;
    fld a;
    //读取数据
    a.load_pram();
    while(m--)
    {
    //数据集分类
    a.load_datat(m);
    //计算均值
    a.load_u_num();
    //计算sw和sb
    a.load_sw_sb();
    //降维，分类
    if(a.last_deal())
    {
        r++;
    AA[a.type()][0]++;
    }
    all++;
    AA[a.type()][1]++;
    }

    for(tag=0;tag<2;tag++)
    A[tag]=AA[tag][0]/AA[tag][1];
    
    printf("\nOA=%.3f\n",r/all);
    printf("AA=%.3f\n",(A[0]+A[1])/2);
    return 0;
}

数据读取，分测试训练集，载入均值的环节和iris一样，故不再赘述

求sw和sb矩阵

二分类时类间散度公式
$$S_b=({\mu }_1-{\mu }_2)({\mu }_1-{\mu }_2{)}^T$$

void fld::load_sw_sb()
{
    for(tag_a=0;tag_a<60;tag_a++)
    for(tag_b=0;tag_b<60;tag_b++)
    sw[tag_a][tag_b]=0;
    for(tag_a=0;tag_a<60;tag_a++)
    for(tag_b=0;tag_b<60;tag_b++)
    for(tag_c=0;tag_c<207;tag_c++)
    sw[tag_a][tag_b]+=((data_train[tag_c][tag_a]-u[int(data_train[tag_c][60])][tag_a])*(data_train[tag_c][tag_b]-u[int(data_train[tag_c][60])][tag_b]));
    for(tag_a=0;tag_a<60;tag_a++)
    for(tag_b=0;tag_b<60;tag_b++)
    sw_[tag_a][tag_b]=sw[tag_a][tag_b];
    double (*p)[60]=sw_;
    brinv(*p,60);
       for(tag_a=0;tag_a<60;tag_a++)
       w[tag_a]=0;
    for(tag_a=0;tag_a<60;tag_a++)
    for(tag_b=0;tag_b<60;tag_b++)
    w[tag_a]+=sw_[tag_a][tag_b]*(u[0][tag_b]-u[1][tag_b]);
}

将测试样本归类并返回是否正确

分类阈值选各维度两类中心的平均值

bool fld::last_deal()
{
    //均值
    double fin;
    //降维后的数据
    double finnn[207];
    for(tag_a=0;tag_a<207;tag_a++)
    {
        fin=0;
        for(tag_b=0;tag_b<60;tag_b++)
        {
            fin+=w[tag_b]*data_train[tag_a][tag_b];
        }
        finnn[tag_a]=fin;
    }
    fin=0;
    //两类在某一维度的均值
    int fin1[2]={0,0};
            for(tag_a=0;tag_a<207;tag_a++)
            fin1[int(data_train[tag_a][60])]+=finnn[tag_a]/207;
        fin=(fin1[0]+fin1[1])*0.5;
    double finn[2]={0,1};
            for(tag_b=0;tag_b<60;tag_b++)
        {
            finn[0]+=w[tag_b]*data_test[0][tag_b];
        }
        finn[1]=data_test[0][60];

        if((finn[0]<fin&&finn[1]==1)||(finn[0]>fin&&finn[1]==0))
        return true;
        else return false;
}