最小二乘法用于多项式的拟合及C语言程序实现

改写自：

https://blog.csdn.net/tutu1583/article/details/82111060

https://blog.csdn.net/piaoxuezhong/article/details/54973750

将C++代码改为C语言

 

一元线性方程可以看做是多元函数的特例，现在用矩阵形式表述多元函数情况下，最小二乘的一般形式：

 设拟合多项式为：

各店到这条曲线的距离之和，即偏差平方和如下：

对等式右边求ai的偏导数，得到：

                       ......

把这些等式表示成矩阵的形式，就可以得到下面的矩阵：

（3）

进行化简计算：

上面公式（3）可以写为：

 

 XY.txt内容：

0    1.0
0.25 1.28
0.5  1.65
0.75 2.12
1    2.72

 

 

 

注意要做防止除数为0的保护

#include <stdio.h>
#include "stdlib.h"
#include "math.h"
//#include <vector>
//using namespace std;

#define MAX_EXP 4 /* x最高次幂 */
#define MATRIX_DIM (MAX_EXP + 1)   /* 矩阵阶数 */
#define SMPNUM 5    /* 采样个数 */
 
struct point
{
    double x;
    double y;
};

/* 采样点想， y */
float sampleX[SMPNUM] = {0.0};
float sampleY[SMPNUM] = {0.0};

void getFile(char *File);  //获取文件数据
void getCoeff(float sampleX[SMPNUM], float sampleY[SMPNUM], float coeff[MATRIX_DIM]);   //获取系数
 
int main()
{
    int i;
    char *File = "XY.txt";
    //vector<point> sample;
    //doubleVector  coefficient;
    //sample = getFile(File);
    getFile(File);
    printf("拟合多项式阶数n=");
    //scanf("%d", &n);
    // coefficient = getCoeff(sample, n);
    //n = 3;
    
    float coeff[MATRIX_DIM] = {0};
    getCoeff(sampleX, sampleY, coeff);
    printf("\n拟合矩阵的系数为：\n");
    for (i = 0; i < MATRIX_DIM; i++)
    {
        printf("a%d = %lf\n", i, coeff[i]);
    }

    return 0;
}
//矩阵方程
void getCoeff(float sampleX[SMPNUM], float sampleY[SMPNUM], float coeff[MATRIX_DIM])
{
    double sum;
    int i, j, k;

    float matX[MATRIX_DIM][MATRIX_DIM];  //公式3左矩阵
    float matY[MATRIX_DIM][MATRIX_DIM];  //公式3右矩阵
    float temp2[MATRIX_DIM];
    //公式3左矩阵
    for (i = 0; i < MATRIX_DIM; i++)
    {
        for (j = 0; j < MATRIX_DIM; j++)
        {
            sum = 0.0;
            for (k = 0; k < SMPNUM; k++)
            {
                sum += pow(sampleX[k], j + i);
            }
            matX[i][j] = sum;
        }
    }
 
    //打印matFunX矩阵
    printf("matFunX矩阵：\n");
    for (i = 0; i < MATRIX_DIM; i++)
    {
        for (j = 0; j < MATRIX_DIM; j++)
        {
            printf("%f\t", matX[i][j]);
        }
        printf("\n");
    }

    //公式3右矩阵
    for (i = 0; i < MATRIX_DIM; i++)
    {
        //temp.clear();
        sum = 0;
        for (k = 0; k < SMPNUM; k++)
        {
            sum += sampleY[k] * pow(sampleX[k], i);
        }
        matY[i][0] = sum;
    }
    //printf("matFunY.size=%d\n", matFunY.size());
    //打印matFunY的矩阵
    printf("matFunY的矩阵：\n");
    for (i = 0; i < MATRIX_DIM; i++)
    {
        printf("%f\n", matY[i][0]); 
    }
    //矩阵行列式变换，将matFunX矩阵变为下三角矩阵,将matFunY矩阵做怎样的变换呢？
    //AX=Y在将X变换为下三角矩阵X'时，是左右两边同乘ratio
    double num1, num2, ratio;
    for (i = 0; i < MATRIX_DIM; i++)
    {
        num1 = matX[i][i];
        for (j = i + 1; j < MATRIX_DIM; j++)
        {
            num2 = matX[j][i];
            ratio = num2 / num1;
            for (k = 0; k < MATRIX_DIM; k++)
            {
                matX[j][k] = matX[j][k] - matX[i][k] * ratio;
            }
            matY[j][0] = matY[j][0] - matY[i][0] * ratio;
        }
    }
    //打印matFunX行列变化之后的矩阵
    printf("matFunX行列变换之后的矩阵：\n");
    for (i = 0; i < MATRIX_DIM; i++)
    {
        for (j = 0; j < MATRIX_DIM; j++)
        {
            printf("%f\t", matX[i][j]);
        }
        printf("\n");
    }
    //打印matFunY行列变换之后的矩阵
    printf("matFunY行列变换之后的矩阵：\n");
    for (i = 0; i < MATRIX_DIM; i++)
    {
        printf("%f\n", matY[i][0]);
    }
    //计算拟合曲线的系数
    //doubleVector coeff(n + 1, 0);
    for (i = MATRIX_DIM - 1; i >= 0; i--)
    {
        if (i == MATRIX_DIM - 1)
        {
            coeff[i] = matY[i][0] / matX[i][i];
        }
        else
        {
            for (j = i + 1; j < MATRIX_DIM; j++)
            {
                matY[i][0] = matY[i][0] - coeff[j] * matX[i][j];
            }
            coeff[i] = matY[i][0] / matX[i][i];
        }
    }
    return;
}
 
//获取文件数据
void getFile(char *File)
{
    int i = 0, j = 0, k = 0;
    //vector<point> dst;

    FILE *fp = fopen(File, "r");
 
    if (fp == NULL)
    {
        printf("Open file error!!!\n");
        exit(0);
    }
 
    point temp;
    double num;
 
    while (fscanf(fp, "%lf", &num) != EOF)
    {
        if (i % 2 == 0)
        {
            temp.x = num;
            sampleX[j++] = num;

            
        }
        else
        {
            temp.y = num;
            sampleY[k++] = num;
        }
        i++;
    }
    fclose(fp);
    return;
    //return dst;
}