fisher.cpp

使用C++编写的经典数值计算算法,包括NEWTON迭带法,SIMPSON,二分法,二分法求根,还附加了模式识别的FISHER算法

#include<iostream.h>
#include<math.h>
#include<ostream.h>
#define N 4
//输出矩阵
void PrintMat(double a[][N],int m,int n)
{
	for(int i=0;i<m;i++){
		cout<<endl;
		for(int j=0;j<n;j++)
			cout<<a[i][j]<<"  ";
	}
}
//求矩阵a的逆矩阵
//调用函数matinv1之前，a[][]保存原来矩阵A
//调用函数matinv1之后，a[][]保存矩阵A的逆矩阵
//n表示矩阵的维数
int matinv1(double a[][N],int n)
{
	int i,j,k,i0;
	double ep=1e-12,w,z,cc;
	int p[N];
	double av[N][N];	
	for(i=1;i<=n;i++)
		for(j=1;j<=n;j++)
			av[i-1][j-1]=0;
	for(i=1;i<=n;i++)
		av[i-1][i-1]=1;
	for(k=1;k<=n;k++){
		w=0;
		for(i=k;i<=n;i++)
			if(fabs(a[i-1][k-1]) > fabs(w)){
					w=a[i-1][k-1];
					i0=i;
					}
		if(fabs(w)<ep)
			return 0;
		if(i0!=k)
			for(j=1;j<=n;j++){
				z=a[i0-1][j-1];
				a[i0-1][j-1]=a[k-1][j-1];
				a[k-1][j-1]=z;

				z=av[i0-1][j-1];
				av[i0-1][j-1]=av[k-1][j-1];
				av[k-1][j-1]=z;
				}
		p[k-1]=i0;
		for(j=1;j<=n;j++){
			a[k-1][j-1]=a[k-1][j-1]/w;
			av[k-1][j-1]=av[k-1][j-1]/w;
		}
		for(i=1;i<=n;i++)
			if(i!=k){
				cc=a[i-1][k-1];
				for(j=1;j<=n;j++){
					a[i-1][j-1]=a[i-1][j-1] - cc*a[k-1][j-1];
					av[i-1][j-1]=av[i-1][j-1] - cc*av[k-1][j-1];
				}
			}
	}// for (k...)
	for(i=1;i<=n;i++)
		for(j=1;j<=n;j++)
			a[i-1][j-1]=av[i-1][j-1];
	return 1;
}
void main()
{
	double X[10][4]={0.276,0.18,0.446,0.683,
		0.378,0.2,0.764,0.673,
		0.325,0.2,0.8,0.633,
		0.138,0.21,0.75,0.728,
		0.29,0.241,0.87,0.649,
		0.27,0.19,1.73,0.613,
		0.45,0.23,2.66,0.544,
		0.076,0.26,0.85,0.733,
		0.346,0.27,1.32,0.621,
		0.186,0.3,0.56,0.796};
    double Y[9][4]={0.62,0.24,6.22,0.544,
		0.61,0.25,1.42,0.494,
		0.62,0.27,1.46,0.51,
		0.56,0.13,1.3,0.372,
		0.56,0.2,3.0,0.221,
		0.29,0.25,4.66,0.395,
		0.302,0.22,3.18,0.25,
		0.347,0.19,17.9,0.23,
		0.269,0.25,8.7,0.145};
    double m1[4]={0,0,0,0};
    double m2[4]={0,0,0,0};
	double S1[4][4]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
	double S2[4][4]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
	double sum=0;
	int i,j,k;
    double x[5][4]={0.302,0.23,1.78,0.59,
		            0.344,0.24,3.4,0.618,
					0.358,0.21,1.37,0.619,
					0.432,0.215,0.9,0.214,
					0.47,0.2,2.9,0.22};
	double SW[4][4];
	double R=0;
	double W[4]={0,0,0,0};	
	for(i=0;i<4;i++)
	{
		for(j=0;j<10;j++)
		{
			m1[i]+=X[j][i];
		}
		m1[i]=m1[i]/10;                      //第一类样本均值
	}
    for(i=0;i<4;i++)
	{
		for(j=0;j<9;j++)
		{
			m2[i]+=Y[j][i];
		}
		m2[i]=m2[i]/9;                       //第二类样本均值
	}
	
    for(j=0;j<10;j++)
	{
		for(i=0;i<4;i++)
		{
			for(k=0;k<4;k++)
				S1[i][k]+=(X[j][i]-m1[i])*(X[j][k]-m1[k]);//第一类离散度		
		}
		
	}
    for(j=0;j<9;j++)
	{
		for(i=0;i<4;i++)
		{
			for(k=0;k<4;k++)
				S2[i][k]+=(Y[j][i]-m2[i])*(Y[j][k]-m2[k]);//第二类离散度	
		}	
	}
	for(i=0;i<4;i++)
	{
		for(j=0;j<4;j++)
		{
			SW[i][j]=S1[i][j]+S2[i][j];            //类内离散度
		}
	}	
    matinv1(SW,4);
	cout<<"SW的逆";
	PrintMat(SW,4,4);
    cout<<'\n';
	for(i=0;i<4;i++)
	{
		m1[i]=m1[i]-m2[i];	
	}
    for(i=0;i<4;i++)
	{
		for(j=0;j<4;j++)
		{
			W[i]+=SW[i][j]*m1[j];
		}
	}
	cout<<"判别向量W："<<'\n';
    for(i=0;i<4;i++)
	{
		cout<<W[i]<<'\n';	
	}
	for(j=0;j<6;j++)
	{
		for(i=0;i<4;i++)
		{
			R+=W[i]*x[j][i];
		}
		if(R>0)cout<<"油层"<<'\n';
		else cout<<"水层"<<'\n';
	}
}