最小二乘法线性拟合函数C语言实现

本方法转自https://github.com/natedomin/polyfit

polyfit.c：

/*------------------------------------------------------------
METHOD:  polyfit
INPUTS:  dependentValues[0..(countOfElements-1)] 	//即xData
         independentValues[0...(countOfElements-1)] //即yData
         countOfElements							//即元素个数
         order - Order of the polynomial fitting	//要拟合的最高阶数
OUTPUTS: coefficients[0..order] - indexed by term	
               (the (coef*x^3) is coefficients[3])
         //拟合的各项系数的结果，x^3 的系数为coefficients[3]
-------------------------------------------------------------*/
int polyfit(const double* const dependentValues ,
            const double* const independentValues ,
            unsigned int        countOfElements ,
            unsigned int        order ,
            double* coefficients)
{
    // Declarations...
    // ----------------------------------
    enum { maxOrder = 5 };

    double B[maxOrder + 1] = {0.0f};
    double P[((maxOrder + 1) * 2) + 1] = {0.0f};
    double A[(maxOrder + 1) * 2 * (maxOrder + 1)] = {0.0f};

    double x , y , powx;

    unsigned int ii , jj , kk;

    // Verify initial conditions....
    // ----------------------------------

    // This method requires that the countOfElements > 
    // (order+1) 
    if (countOfElements <= order)
        return -1;

    // This method has imposed an arbitrary bound of
    // order <= maxOrder.  Increase maxOrder if necessary.
    if (order > maxOrder)
        return -1;

    // Begin Code...
    // ----------------------------------

    // Identify the column vector
    for (ii = 0; ii < countOfElements; ii++)
    {
        x = dependentValues[ii];
        y = independentValues[ii];
        powx = 1;

        for (jj = 0; jj < (order + 1); jj++)
        {
            B[jj] = B[jj] + (y * powx);
            powx = powx * x;
        }
    }

    // Initialize the PowX array
    P[0] = countOfElements;

    // Compute the sum of the Powers of X
    for (ii = 0; ii < countOfElements; ii++)
    {
        x = dependentValues[ii];
        powx = dependentValues[ii];

        for (jj = 1; jj < ((2 * (order + 1)) + 1); jj++)
        {
            P[jj] = P[jj] + powx;
            powx = powx * x;
        }
    }

    // Initialize the reduction matrix
    //
    for (ii = 0; ii < (order + 1); ii++)
    {
        for (jj = 0; jj < (order + 1); jj++)
        {
            A[(ii * (2 * (order + 1))) + jj] = P[ii + jj];
        }

        A[(ii * (2 * (order + 1))) + (ii + (order + 1))] = 1;
    }

    // Move the Identity matrix portion of the redux matrix
    // to the left side (find the inverse of the left side
    // of the redux matrix
    for (ii = 0; ii < (order + 1); ii++)
    {
        x = A[(ii * (2 * (order + 1))) + ii];
        if (x != 0)
        {
            for (kk = 0; kk < (2 * (order + 1)); kk++)
            {
                A[(ii * (2 * (order + 1))) + kk] =
                    A[(ii * (2 * (order + 1))) + kk] / x;
            }

            for (jj = 0; jj < (order + 1); jj++)
            {
                if ((jj - ii) != 0)
                {
                    y = A[(jj * (2 * (order + 1))) + ii];
                    for (kk = 0; kk < (2 * (order + 1)); kk++)
                    {
                        A[(jj * (2 * (order + 1))) + kk] =
                            A[(jj * (2 * (order + 1))) + kk] -
                            y * A[(ii * (2 * (order + 1))) + kk];
                    }
                }
            }
        }
        else
        {
            // Cannot work with singular matrices
            return -1;
        }
    }

    // Calculate and Identify the coefficients
    for (ii = 0; ii < (order + 1); ii++)
    {
        for (jj = 0; jj < (order + 1); jj++)
        {
            x = 0;
            for (kk = 0; kk < (order + 1); kk++)
            {
                x = x + (A[(ii * (2 * (order + 1))) + (kk + (order + 1))] *
                    B[kk]);
            }
            coefficients[ii] = x;
        }
    }

    return 0;
}

验证该函数的正确性，代码如下：
main.c

#include <stdbool.h>
#include <stdint.h>
#include <stdio.h>

#define ELE_COUNT 5
int main(void)
{
  const unsigned int order = 3;
  int result;
  // These inputs should result in the following approximate coefficients:
  //         0.5           2.5           1.0        3.0
  //    y = (0.5 * x^3) + (2.5 * x^2) + (1.0 * x) + 3.0
  double    xData[ELE_COUNT] = {12.0,    77.8,      44.1,     23.6,    108.2};
  double    yData[ELE_COUNT] = {1239.00, 250668.38, 47792.19, 7991.13, 662740.98};
  double coefficients[order + 1]; // resulting array of coefs

  // Perform the polyfit
  result = polyfit(xData ,
                   yData ,
                   ELE_COUNT ,
                   order ,
                   coefficients);
  printf("\n------------------------------\npolyfit result is:\n");
  printf("%lfx^3+%lfx^2+%lfx+%f" , coefficients[3],coefficients[2],coefficients[1],
  								   coefficients[0]);
  printf("\n------------------------------\n");
  return 0;
}

进行编译运行后，结果如下，可看出结果基本一致，也可自行拿其他函数进行验证。