pls.cpp

偏最小二乘算法代码及相关说明,在机器学习,实时数值仿真中用得较多

#include <iostream>		
#include "Comm.h"		
#include "Matrix.h"		
#include <stdlib.h>
#include <time.h>
#include <vector>
#include <string.h>
#include <math.h>
using namespace std;	
const int N=21;
const int M=5;
const int L=6;
const int K=20;
const int tr=15;
struct arr_type
{
	int i; //row
	int j; //col
	double value;
	double hvalue()
	{ return value;}
}h[N][L];


void main()   
{ 
	clock_t start, finish; 
    double duration; 
	start = clock();//开始时间

	FILE *in=fopen("E:\\wenyufeng\\pls_regression\\result\\orignaldata.txt","r");
    FILE *out=fopen("E:\\wenyufeng\\pls_regression\\result\\pls_out.txt","w");
   int i,j,k,g,ch,ch1,ch2;
	double X[N][M];
    double Y[N];
    for(i=0;i<N;i++)
	  for(j=0;j<(M+1);j++)
		{
		  fscanf(in,"%lf",&(h[i][j].value));

          if(j==M)
		  {
			  Y[i]=h[i][j].hvalue();
			 
		  }
		  else 
		  {X[i][j]=h[i][j].hvalue();}

         }

	matrix<double> mX(&X[0][0],N,M);
    matrix<double> mY(&Y[0],N,1);
	matrix<double> E(N,M);
	matrix<double> F(N,1);
	matrix<double> trE1(tr,M);
	matrix<double> trF1(tr,1);
	matrix<double> teE1(N-tr,M);
	matrix<double> teF1(N-tr,1);



    cout.setf(ios::fixed);	//输出数据为定点法
	cout.precision(6);//精度6位

	double eps=DOUBLEERROR;
	double xmean[M],xstddev[M];
	double ymean,ystddev;
	for(j=0;j<M;j++)
	{
		xmean[j]=0.0;
        for(i=0;i<N;i++)
		{xmean[j] += mX(i,j);}
		xmean[j] /= (double)N;
	}
    for(j=0;j<M;j++)
	{
		xstddev[j]=0.0;
        for(i=0;i<N;i++)
		{xstddev[j] += (mX(i,j)-xmean[j])*(mX(i,j)-xmean[j]);}
		xstddev[j] /=(double)N;
        xstddev[j]=sqrt(xstddev[j]);
       if(xstddev[j]<=eps) xstddev[j]=1.0;
	}
   for(i=0;i<N;i++)
      for(j=0;j<M;j++)
	  {E(i,j)=(mX(i,j)-xmean[j])/xstddev[j];}

	  matrix<double> trE(tr,M);
	  matrix<double> teE(N-tr,M);
	  for(i=0;i<tr;i++)
	  {
		  for(j=0;j<M;j++)
		  {trE(i,j)=E(i,j);}
	  }
	  for(i=0;i<N-tr;i++)
	  {
		  for(j=0;j<M;j++)
		  {teE(i,j)=E((i+tr),j);}
	  }
	  
    ymean=0.0;
	for(i=0;i<N;i++)
	{ymean += mY(i,0);}
    ymean /=(double)N;

	ystddev=0.0;
	for(i=0;i<N;i++)
	{ystddev += (mY(i,0)-ymean)*(mY(i,0)-ymean);}
	ystddev /= (double)N;
	ystddev=sqrt(ystddev);
	for(i=0;i<N;i++)
	{F(i,0)=(mY(i,0)-ymean)/ystddev;}

	matrix<double> trF(tr,1);
	matrix<double> teF(N-tr,1);
	for(i=0;i<tr;i++)
	{trF(i,0)=F(i,0);}
	for(i=0;i<N-tr;i++)
	{teF(i,0)=F((i+tr),0);}
/*
	double oldsse,sse;
	oldsse=0.0;
	for(i=0;i<tr;i++)
	{
		oldsse += trF(i,0)*trF(i,0);
	}
*/	
	matrix<double>teY(N-tr,1);
	for(i=0;i<N-tr;i++)
	{teY(i,0)=0;}
	

loop3:
	matrix<double> u0(tr,1);
	for(i=0;i<tr;i++)
	{u0(i,0)=trF(i,0);}//u0=F

    matrix<double> w0t(1,tr);
	double sumu0=0.0;
	for(i=0;i<tr;i++)
	{sumu0 += u0(i,0)*u0(i,0);}

    matrix<double> u0t(1,tr);
    MatrixTranspose(u0,u0t);

    MatrixMultiply(w0t,u0t,trE);
	w0t /= sumu0;

    matrix<double> w0(1,tr);
    MatrixTranspose(w0t,w0);
    
    matrix<double> sumw0(0,0);
    MatrixMultiply(sumw0,w0t,w0);

	w0t /= sqrt(sumw0(0,0));

    MatrixTranspose(w0t,w0);

	matrix<double> sumww0(0,0);
   MatrixMultiply(sumww0,w0t,w0);

    matrix<double> t0(tr,1);
    MatrixMultiply(t0,trE,w0);
	t0 /= sumww0(0,0);

    double q=1.0;

    matrix<double> p0t(1,tr);
    matrix<double> t0t(1,tr);
    MatrixTranspose(t0,t0t);

    matrix<double> sumt0(0,0);
	MatrixMultiply(sumt0,t0t,t0);

	MatrixMultiply(p0t,t0t,trE);
	p0t /= sumt0(0,0);

    matrix<double> p0(tr,1);
    MatrixTranspose(p0t,p0);

    matrix<double> sump0(0,0);
	MatrixMultiply(sump0,p0t,p0);

	p0t /= sqrt(sump0(0,0));

	t0 *=sqrt(sump0(0,0));

	w0 *= sqrt(sump0(0,0));
   
	matrix<double> tt(0,0);
	MatrixTranspose(t0,t0t);
	MatrixMultiply(tt,t0t,t0);

	matrix<double> b0(0,0);
	MatrixMultiply(b0,u0t,t0);
	b0 /= tt(0,0);

	matrix<double> _Y(N-tr,1);
	MatrixMultiply(_Y,teE,w0);
	_Y *= q;
	_Y *=b0;
	teY += _Y;
 
	MatrixMultiply(trE1,t0,p0t);
	trE -= trE1;
	trF1 = t0*q*b0;
	trF -= trF1;
	MatrixMultiply(teE1,teE,w0);
	MatrixMultiply(teE1,teE1,p0t);
	teE -= teE1;
/*
	sse=0.0;
	for(i=0;i<tr;i++)
	{
		sse += trF1(i,0)*trF1(i,0);
	}

	double r=0.0;
	r=(oldsse-sse)/oldsse;
	if(r<0.99) goto loop3;
*/		

    double y[N-tr];
    for(i=0;i<N-tr;i++)
    {y[i]=teY(i,0)*ystddev+ymean;}
	double re=0.0;
	for(i=0;i<N-tr;i++)
	{
		double rer=abs(y[i]-Y[i+tr]);
		double rerr=rer/abs(Y[N-tr]);
		re += rerr;
	}
	re /= (double)(N-tr);
	re *= 100.0;
	if(re>500.0) goto loop3;
	int correct = 0;
	int total = 0;
	double error = 0,a_error=0,r_error=0;
	double sumv = 0, sumy = 0, sumvv = 0, sumyy = 0, sumvy = 0;
	double Z,hh,T,z0=0;


	for( i=0;i<N-tr;i++)
	{
		if((y[i]-Y[i+tr])==0)
	
             ++correct;

		   error += (y[i]-Y[i+tr])*(y[i]-Y[i+tr]);

	        if((y[i]-Y[i+tr])>0) Z=(y[i]-Y[i+tr]);
			else Z=-(y[i]-Y[i+tr]);  
            a_error += Z;
			if(Y[i+tr]>0){ hh= Z/Y[i+tr];}
			else {hh= Z/(-Y[i+tr]);}
			fprintf(out,"relative error=%g, y_predict(%d)=%g, y_target(%d)=%g\n",hh*100.0,i,y[i],i,Y[i+tr]);
			r_error += hh;

		sumv += y[i];
		sumy += Y[i+tr];
		sumvv += y[i]*y[i];
		sumyy += Y[i+tr]*Y[i+tr];
		sumvy += y[i]*Y[i+tr];
		++total;
	}
       fprintf(out,"\n");
		if(total==N-tr)
		{
	     fprintf(out,"Mean squared error = %g (regression)\n",error/total);
		 fprintf(out,"Squared correlation coefficient = %g (regression)\n",
		       ((total*sumvy-sumv*sumy)*(total*sumvy-sumv*sumy))/
		       ((total*sumvv-sumv*sumv)*(total*sumyy-sumy*sumy))
		       );
         fprintf(out,"Mean absoute error= %g (regression)\n",a_error/total);
         fprintf(out,"Mean relative error= %g (regression)\n",100.0*r_error/total);
		}
		fclose(in);
		fclose(out);

	finish = clock(); //结束时间
    duration = (double)(finish - start) / CLOCKS_PER_SEC; //运行时间
    printf( "\n%f seconds\n", duration );
	FILE *tm=fopen("E:\\wenyufeng\\pls_regression\\result\\time.txt","w");
	fprintf(tm,"the run-duration of program :\n");
    fprintf( tm,"%f seconds\n", duration );

}