Parameters are set indicating the number of variables NVAR, and the number of data points NPTS to use in the curve-fit. The equation structure is defined in the aU array for each of the variables and points. The data to curve-fit is put in the dependent vector aV, generated by known coefficients aC0 which should be recovered by the analysis in array aC. The least-squares function is called, and then the data is printed out.
00001 00002 #include "Cpx.h" // complex variable definitions 00003 #include "Mtx.h" // matrix definitions 00004 00005 using namespace std; // standard C++ library namespace 00006 using namespace calc; // CalcLib namespace 00007 00008 typedef Cpx<float> CPXf; // alias for float type complex numbers 00009 typedef Mtx<CPXf > MTXc; // alias for complex matrices 00010 00011 // results print function 00012 void results(ostream& ostr, const char *cTle, const CPXf *a, int n) 00013 {ostr<<" "<<cTle<<" "; 00014 for(int i=0; i<n; i++) ostr<<" "<<a[i]; ostr<<endl;} 00015 00016 // least-squares curve-fit function 00017 void leastSquares(CPXf *c, const CPXf *u, const CPXf *v, int npts, int nvar) 00018 { 00019 // convert data arrays into matrix objects 00020 MTXc mU(u,npts,nvar); 00021 MTXc mV(v,npts,1); 00022 00023 // set up augmented matrix for Ax=B 00024 MTXc mA=(~mU)*mU; // build square matrix: transpose(U)*U = A 00025 MTXc mB=(~mU)*mV; // build coefficient vector: transpose(U)*V = B 00026 00027 // solve matrix with inverse of matrix A times B 00028 MTXc mC=(!mA)*mB; 00029 mC.copyArray(c); // best-fit coefficients in C as solution 00030 } 00031 00032 int main(void) 00033 { 00034 // define matrix size 00035 const int NPTS=40; // number of data points to curve-fit 00036 const int NVAR=4; // number of variables to use in curve-fit 00037 00038 // set coefficients for check 00039 const CPXf aC0[NVAR]= // least-squares method should recover these values 00040 {CPXf(-9.3f, 9.3f),CPXf( 1.1f,-1.1f), 00041 CPXf(-5.1f, 5.1f),CPXf( 0.0f, 0.0f)}; 00042 00043 // declare solution variables 00044 CPXf aU[NPTS][NVAR]; // independent variable data array 00045 CPXf aV[NPTS]; // dependent variable data vector 00046 CPXf aC[NVAR]; // least-squares coefficient array 00047 00048 // build data arrays 00049 for(int i=0; i<NPTS; i++) { 00050 const CPXf x((float)(i),(float)(-i)); 00051 aU[i][0]=1.0f; // matrix defines best-fit equation structure 00052 aU[i][1]=x; // equation is 3rd degree polynomial 00053 aU[i][2]=x*x; 00054 aU[i][3]=x*x*x; 00055 aV[i]=aC0[0]+aC0[1]*x+aC0[2]*x*x+aC0[3]*x*x*x; // dependent data vector 00056 } // usually generated from outside source 00057 00058 // curve-fit aV data to 3rd degree polynomial in aU 00059 leastSquares(aC,&aU[0][0],aV,NPTS,NVAR); 00060 00061 // print out coefficients 00062 cout.precision(5); 00063 cout.setf(ios::showpoint,ios::showpoint); 00064 cout<<endl<<" calc::Mtx Class Example Application"<<endl<<endl; 00065 results(cout,"Initial Coefficients: ",aC0,NVAR); 00066 results(cout,"Least-Squares Coefficients: ",aC ,NVAR); 00067 cout<<endl<<flush; 00068 00069 return 0; 00070 } 00071
The following output file shows the approximate agreement betweem the original coefficients aC0 and the least-squares produced coefficients aC. Of course, this is an example; finding the complete inverse of a square matrix is not the most efficient way to solve a least-squares problem.
00001 00002 calc::Mtx Class Example Application 00003 00004 Initial Coefficients: -9.3000,9.3000 1.1000,-1.1000 -5.1000,5.1000 0.0000,0.0000 00005 Least-Squares Coefficients: -9.1875,9.4062 1.1030,-1.1406 -5.1006,5.1008 7.6294e-06,1.6093e-06 00006
1.4.7