pls.h

最小二乘回归程序 并且举例测试 希望对大家有帮助

#include	<stdio.h>
#include	<math.h>
#include	<string.h>
#include	<stdlib.h>

/* Local Include Declaration */
//#include	"AnalysisParameterCommon.h"				/* common data header */

/******************************************************************
*
*  MODULE	  : lineq
*  ABSTRACT	  : Automatic Regression Calculation
*  FUNCTION	  : Regression calculation function
*  NOTE
*  RETURN	  :  0: Normal Return
*				-1: Abnormal Return
*
*  CREATE	  : 2003/09/08 H.Tanaka
*  UPDATE	  : 
*
******************************************************************/

/*  (C) 2002-2003 MITSUBISHI ELECTRIC CORPORATION */


  bool gaussj(float A[33][33], float B[33][33],int n)
{
 int i,j,k;
 float lMax,temp;

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

 //初始化B为单位阵
 for(i=0;i<n;i++)
 {
  for(j=0;j<n;j++)
  {
   if(i!=j)B[i][j] = 0;
   else B[i][j] = 1;
  }
 }
 for(i=0;i<n;i++)
 {
  //寻找主元
  lMax = tt[i][i];
  k = i;
  for(j=i+1;j<4;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<n;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<n;j++)
  {
   tt[i][j] = tt[i][j] / temp; //主对角线上元素变成1
   B[i][j] = B[i][j] / temp; //伴随计算
  }

  for(j=0;j<n;j++) //行 0 -> n
  {
   if(j!=i)  //不是第i行
   {
    temp = tt[j][i];
    for(k=0;k<4;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;
}





int	lineq(
			float	aa[32][32],				/* (I) X */
			float	*bb,					/* (I) Y  (O) A */	//  X * A = Y --> A = X(-1) * Y
			int		*ip,					/* (I) working buffer */
			int		n						/* (I) count of analysis parameter */
	)                                     
{
	int		i;
//	int		k;
	int		kp1;
	int		m;
	int		j;
//	int		nm1;
//	int		km;
//	int		p;
	 float rp[33][33];
//	float	s;
//	float	t;
//	float	s;
//	float	t;
	float	a[33][33];  //it's x'y
	float	b[32];
	float  aaa[32][33];//it's x
	float  bbb[32];    //it`s y
	float  ab[33][32]; //it`s x'
	float  abb[33][32];    //it's (x'x)^-1*x'
	//x左面需要加[1111……1]'
	memset( &rp, 0, sizeof(rp) );
	memset( &a, 0, sizeof(a) );
	memset( &b, 0, sizeof(b) );
	memset( &aaa, 0, sizeof(aaa) );
	memset( &bbb, 0, sizeof(bbb) );
	memset( &ab, 0, sizeof(ab) );
	memset( &abb, 0, sizeof(abb) );

	for(i=0;i<n+1;i++){
		aaa[i][0]=1;
	}


	for ( i=0; i<n; i++ ) {
		for ( j=0; j<n; j++ ) {
			aaa[i][j+1] = aa[i][j];
		}
		bbb[i] = bb[i];
	}


	for ( i=0; i<n; i++ )         //calc x'
		for ( j=0; j<n+1; j++ ) 
			ab[j][i] = aaa[i][j];


    for(i=0;i<n+1;i++)                //calc x'x
	{
		for(m=0;m<n+1;m++)
		{   
			for(kp1=0;kp1<n;kp1++)
			a[i][m]+=ab[i][kp1]*aaa[kp1][m];
			
		}
	}

///////////////////

	    for(i=0;i<4;i++)
        	for(j=0;j<4;j++)
		{

		  printf("%f ",a[i][j]);
		  if(j==3)
			  printf("\n");
		}
printf("\n\n");
//////////////////////

gaussj(a, rp,n+1);
	//Converse(aaaaaa, rp,n+1) ;//现在rp是 x'x的逆，之后求出rp*x'，再*y就ok了

	//////////////////////////////////
	    for(i=0;i<4;i++)
        	for(j=0;j<4;j++)
		{
			printf("%f   ",rp[i][j]);
	        if(j==3)
	     	printf("\n");
		}

	printf("\n\n");

///////////////////////////////////
memset(&abb,0,sizeof(abb));
    for(i=0;i< n+1;i++)                //calc (x'x)^-1*x'
	{
		for(j=0;j<n;j++)
		{   
			for(kp1=0;kp1<n+1;kp1++)
			abb[i][j]+=rp[i][kp1]*ab[kp1][j];
			
		}

	}

	
memset(&b,0,sizeof(b));

	//算好abb*bbb就完事了
	for(i=0;i<n+1;i++)
		for(j=0;j<n;j++)
		  b[i]+=abb[i][j]*bbb[j];




	for ( i=0; i<n+1; i++ ) 
	bb[i] = (float)b[i];
	

	return(0);
}