cpolyfit.cpp

多项式曲线拟合 任意介数 Purpose - Least-squares curve fit of arbitrary order working in C++ Builder 2007 as

//---------------------------------------------------------------------------------

template <class FloatType>
bool __fastcall CPolyFit<FloatType>::GaussJordan (Matrix &b,            
                                                  const vector<FloatType> &y,
                                                  vector<FloatType> &coef)
                                                  //b square matrix of coefficients
                                                  //y constant vector
                                                  //coef solution vector
                                                  //ncol order of matrix got from b.size()


{
/*
 { Gauss Jordan matrix inversion and solution }

 { B (n, n) coefficient matrix becomes inverse }
 { Y (n) original constant vector }
 { W (n, m) constant vector(s) become solution vector }
 { DETERM is the determinant }
 { ERROR = 1 if singular }
 { INDEX (n, 3) }
 { NV is number of constant vectors }
*/

 try
 {
  size_t ncol(b.size());
  size_t irow, icol;
  vector<vector<int> >index;
  Matrix w;

  Zeroise(w, ncol, ncol);
  Zeroise(index, ncol, 3);

  if(!GaussJordan2(b, y, w, index))
   return false;

  // Interchange columns
  size_t m;
  for (size_t i = 0; i <  ncol; ++i)
  {
    m = ncol - i - 1;
    if(index [m][0] != index [m][1])
    {
     irow = index [m][0];
     icol = index [m][1];
     for(size_t k = 0; k < ncol; ++k)
      Swap (b[k][irow], b[k][icol]);
    }
  }

  for(size_t k = 0; k < ncol; ++k)
  {
   if(index [k][2] != 0)
   {
    throw "Error in GaussJordan: matrix is singular";
   }
  }

  for(size_t i = 0; i < ncol; ++i)
   coef[i] = w[i][0];
 }
 catch(char *b)
 {
  Application->MessageBox(b, "PolyFit Exception",
   MB_OK | MB_ICONEXCLAMATION | MB_DEFBUTTON1 | MB_APPLMODAL | MB_SETFOREGROUND | MB_TOPMOST);
  return false;
 }
 return true;
}   // end;	{ procedure GaussJordan }
//----------------------------------------------------------------------------------------------

template <class FloatType>
bool __fastcall CPolyFit<FloatType>::GaussJordan2(Matrix &b,
                                                  const vector<FloatType> &y,
                                                  Matrix &w,
                                                  vector<vector<int> > &index)
{
 //GaussJordan2;         // first half of GaussJordan
  // actual start of gaussj
 try
 {

  FloatType big, t;
  FloatType pivot;
  FloatType determ;
  size_t irow, icol;
  size_t ncol(b.size());
  size_t nv = 1;                  // single constant vector
  for(size_t i = 0; i < ncol; ++i)
  {
   w[i][0] = y[i];      // copy constant vector
   index[i][2] = -1;
  }
  determ = 1.0;
  for(size_t i = 0; i < ncol; ++i)
  {
   // Search for largest element
   big = 0.0;
   for(size_t j = 0; j < ncol; ++j)
   {
    if(index[j][2] != 0)
    {
     for(size_t k = 0; k < ncol; ++k)
     {
      if(index[k][2] > 0)
       throw "Error in GaussJordan2: matrix is singular";

      if(index[k][2] < 0 && fabs(b[j][k]) > big)
      {
       irow = j;
       icol = k;
       big = fabs(b[j][k]);
      }
     } //   { k-loop }
    }
   }  // { j-loop }
   index [icol][2] = index [icol][2] + 1;
   index [i][0] = irow;
   index [i][1] = icol;

   // Interchange rows to put pivot on diagonal
   // GJ3
   if(irow != icol)
   {
    determ = -determ;
    for(size_t m = 0; m < ncol; ++m)
     Swap (b [irow][m], b[icol][m]);
    if (nv > 0)
     for (size_t m = 0; m < nv; ++m)
      Swap (w[irow][m], w[icol][m]);
   } // end GJ3

    // divide pivot row by pivot column
   pivot = b[icol][icol];
   determ *= pivot;
   b[icol][icol] = 1.0;

   for(size_t m = 0; m < ncol; ++m)
    b[icol][m] /= pivot;
   if(nv > 0)
    for(size_t m = 0; m < nv; ++m)
     w[icol][m] /= pivot;

   // Reduce nonpivot rows
   for(size_t n = 0; n < ncol; ++n)
   {
    if(n != icol)
    {
     t = b[n][icol];
     b[n][icol] = 0.0;
     for(size_t m = 0; m < ncol; ++m)
      b[n][m] -= b[icol][m] * t;
     if(nv > 0)
      for(size_t m = 0; m < nv; ++m)
       w[n][m] -= w[icol][m] * t;
    }
   }
  } // { i-loop }
 }
 catch(char *b)
 {
  Application->MessageBox(b, "PolyFit Exception",
   MB_OK | MB_ICONEXCLAMATION | MB_DEFBUTTON1 | MB_APPLMODAL | MB_SETFOREGROUND | MB_TOPMOST);
  return false;
 }
 return true;
 //end;		{ GaussJordan2 }
}
//----------------------------------------------------------------------------------------------

 // Least-squares fit to nrow sets of x and y pairs of points

template <class FloatType>
void __fastcall CPolyFit<FloatType>::LinFit
                                 (const vector<FloatType> &x,       // independent variable
                                  vector<FloatType> &y,             // dependent variable
                                  vector<FloatType> &y_calc,	       // calculated dep. variable
                                  vector<FloatType> &resid,         // array of residuals
                                  vector<FloatType> &coef,	         // coefficients
                                  vector<FloatType> &sig,           // error on coefficients
                                  const size_t nrow,             // length of variable array
                                  size_t ncol)                   // number of terms
{
// bool error;
 size_t i, j, nm;
 FloatType xi, yi, yc, srs, see, sum_y, sum_y2;
 Matrix xmatr;        // Data matrix
 Matrix a;
 vector<FloatType> g;      // Constant vector

 Zeroise(g, ncol);
 Zeroise(a, ncol, ncol);
 Zeroise(xmatr, nrow, ncol);
 for(i = 0; i < nrow; ++i)
 {
  // { setup x matrix }
  xi = x[i];
  xmatr[i][0] = 1.0;	      //  { first column }
  for(j = 1; j < ncol; ++j)
   xmatr [i][j] = xmatr [i][j - 1] * xi;
 }
 Square (xmatr, y, a, g, nrow, ncol);
 GaussJordan (a, g, coef);
 sum_y = 0.0;
 sum_y2 = 0.0;
 srs = 0.0;
 for(i = 0; i < nrow; ++i)
 {
  yi = y[i];
  yc = 0.0;
  for(j = 0; j < ncol; ++j)
   yc +=  coef [j] * xmatr [i][j];
  y_calc[i] = yc;
  resid[i] = yc - yi;
  srs += Sqr (resid [i]);
  sum_y +=  yi;
  sum_y2 += yi * yi;
 }

 if(nrow == ncol)
  nm = 1;
 else nm = nrow - ncol;

 see = sqrt (srs / nm);
 for(i = 0; i < ncol; ++i) //		{ errors on solution }
  sig[i] = see * sqrt (a [i][i]);
}
//------------------------------------------------------------------------------------

// Utility functions

// fills a vector with zeros.
template <typename X>
void NSUtility::Zeroise(vector<X> &array, size_t n)
{
 array.clear();
 for(size_t j = 0; j < n; ++j)
  array.push_back(0);
}
//--------------------------------------------------------------------------

// fills a (m by n) matrix with zeros.
template <typename X>
void NSUtility::Zeroise(vector<vector<X> > &matrix, size_t m, size_t n)
{
 vector<X> zero;
 Zeroise(zero, n);
 matrix.clear();
 for(size_t j = 0; j < m; ++j)
  matrix.push_back(zero);
}
//--------------------------------------------------------------------------

// calcs x^i
template <typename X>
X NSUtility::Pow(const X &x, int i)
{
 X ret = 1.0;
 X t = i < 0 ? 1.0 / x : x;
 int n = abs(i);

 for(i = 0; i < n; ++i)
   ret *= t;
 return ret;
}
//--------------------------------------------------------------------------