最小二乘曲线拟合——C语言算法实现一

最小二乘曲线拟合

 

给定一组数据,我们要对这组数据进行曲线拟合。

假定要拟合的曲线方程为 y=a0 + a1*x^1 + a2*x^2 + a3*x^3 + ...+ an*x^n

    x             y

0.995119      -7.620000
2.001185      -2.460000
2.999068     108.760000
4.001035     625.020000
4.999859     2170.500000
6.004461     5814.580000
6.999335     13191.840000
7.999433     26622.060000
9.002257      49230.220000
10.003888    85066.500000
11.004076    139226.280000
12.001602     217970.140000
13.003390     328843.860000
14.001623    480798.420000
15.003034    684310.000000
16.002561    951499.980000
17.003010   1296254.940000
18.003897   1734346.660000
19.002563   2283552.120000
20.003530    2963773.500000

代码如下:

#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include <time.h>
#define ParaBuffer(Buffer,Row,Col) (*(Buffer + (Row) * (SizeSrc + 1) + (Col)))
 /***********************************************************************************
 ***********************************************************************************/
 static int GetXY(const char* FileName, double* X, double* Y, int* Amount)
 {
     FILE* File = fopen(FileName, "r");
     if (!File)
		return -1;
     for (*Amount = 0; !feof(File); X++, Y++, (*Amount)++)
         if (2 != fscanf(File, (const char*)"%lf %lf", X, Y))
            break;
     fclose(File);
     return 0;
 }
 /***********************************************************************************
 ***********************************************************************************/
 static int PrintPara(double* Para, int SizeSrc)
 {
     int i, j;
     for (i = 0; i < SizeSrc; i++)
     {
         for (j = 0; j <= SizeSrc; j++)
         printf("%10.6lf ", ParaBuffer(Para, i, j));
         printf("\r\n");
     }
     printf("\r\n");
     return 0;
 }
 /***********************************************************************************
 ***********************************************************************************/
 static int ParalimitRow(double* Para, int SizeSrc, int Row)
 {
     int i;
     double Max, Min, Temp;
     for (Max = abs(ParaBuffer(Para, Row, 0)), Min = Max, i = SizeSrc; i; i--)
     {
         Temp = abs(ParaBuffer(Para, Row, i));
         if (Max < Temp)
             Max = Temp;
         if (Min > Temp)
             Min = Temp;
     }
     Max = (Max + Min) * 0.000005;
     for (i = SizeSrc; i >= 0; i--)
         ParaBuffer(Para, Row, i) /= Max;
     return 0;
 }
 /***********************************************************************************
 ***********************************************************************************/
 static int Paralimit(double* Para, int SizeSrc)
 {
     int i;
     for (i = 0; i < SizeSrc; i++)
          if (ParalimitRow(Para, SizeSrc, i))
             return -1;
     return 0;
 }
 /***********************************************************************************
 ***********************************************************************************/
 static int ParaPreDealA(double* Para, int SizeSrc, int Size)
 {
     int i, j;
     for (Size -= 1, i = 0; i < Size; i++)
     {
         for (j = 0; j < Size; j++)
         ParaBuffer(Para, i, j) = ParaBuffer(Para, i, j) * ParaBuffer(Para, Size, Size) - ParaBuffer(Para, Size, j) * ParaBuffer(Para, i, Size);
         ParaBuffer(Para, i, SizeSrc) = ParaBuffer(Para, i, SizeSrc) * ParaBuffer(Para, Size, Size) - ParaBuffer(Para, Size, SizeSrc) * ParaBuffer(Para, i, Size);
         ParaBuffer(Para, i, Size) = 0;
         ParalimitRow(Para, SizeSrc, i);
     }
     return 0;
 }
 /***********************************************************************************
 ***********************************************************************************/
 static int ParaDealA(double* Para, int SizeSrc)
 {
     int i;
     for (i = SizeSrc; i; i--)
         if (ParaPreDealA(Para, SizeSrc, i))
             return -1;
     return 0;
 }
 /***********************************************************************************
 ***********************************************************************************/
 static int ParaPreDealB(double* Para, int SizeSrc, int OffSet)
 {
     int i, j;
     for (i = OffSet + 1; i < SizeSrc; i++)
     {
         for (j = OffSet + 1; j <= i; j++)
         ParaBuffer(Para, i, j) *= ParaBuffer(Para, OffSet, OffSet);
         ParaBuffer(Para, i, SizeSrc) = ParaBuffer(Para, i, SizeSrc) * ParaBuffer(Para, OffSet, OffSet) - ParaBuffer(Para, i, OffSet) * ParaBuffer(Para, OffSet, SizeSrc);
         ParaBuffer(Para, i, OffSet) = 0;
         ParalimitRow(Para, SizeSrc, i);
     }
     return 0;
 }
 /***********************************************************************************
 ***********************************************************************************/
 static int ParaDealB(double* Para, int SizeSrc)
 {
     int i;
     for (i = 0; i < SizeSrc; i++)
             if (ParaPreDealB(Para, SizeSrc, i))
                     return -1;
     for (i = 0; i < SizeSrc; i++)
     {
         if (ParaBuffer(Para, i, i))
         {
             ParaBuffer(Para, i, SizeSrc) /= ParaBuffer(Para, i, i);
             ParaBuffer(Para, i, i) = 1.0;
         }
     }
     return 0;
 }
 /***********************************************************************************
 ***********************************************************************************/
 static int ParaDeal(double* Para, int SizeSrc)
 {
     PrintPara(Para, SizeSrc);
     Paralimit(Para, SizeSrc);
     PrintPara(Para, SizeSrc);
     if (ParaDealA(Para, SizeSrc))
             return -1;
     PrintPara(Para, SizeSrc);
     if (ParaDealB(Para, SizeSrc))
             return -1;
     return 0;
 }
 /***********************************************************************************
 ***********************************************************************************/
 static int GetParaBuffer(double* Para, const double* X, const double* Y, int Amount, int SizeSrc)
 {
     int i, j;
     for (i = 0; i < SizeSrc; i++)
             for (ParaBuffer(Para, 0, i) = 0, j = 0; j < Amount; j++)
                     ParaBuffer(Para, 0, i) += pow(*(X + j), 2 * (SizeSrc - 1) - i);
     for (i = 1; i < SizeSrc; i++)
             for (ParaBuffer(Para, i, SizeSrc - 1) = 0, j = 0; j < Amount; j++)
                     ParaBuffer(Para, i, SizeSrc - 1) += pow(*(X + j), SizeSrc - 1 - i);
     for (i = 0; i < SizeSrc; i++)
             for (ParaBuffer(Para, i, SizeSrc) = 0, j = 0; j < Amount; j++)
                     ParaBuffer(Para, i, SizeSrc) += (*(Y + j)) * pow(*(X + j), SizeSrc - 1 - i);
     for (i = 1; i < SizeSrc; i++)
             for (j = 0; j < SizeSrc - 1; j++)
                     ParaBuffer(Para, i, j) = ParaBuffer(Para, i - 1, j + 1);
     return 0;
 }
 //***********************************************************************************
 //***********************************************************************************
 int Cal(const double* BufferX, const double* BufferY, int Amount, int SizeSrc, double* ParaResK)
 {
     double* ParaK = (double*)malloc(SizeSrc * (SizeSrc + 1) * sizeof(double));
     GetParaBuffer(ParaK, BufferX, BufferY, Amount, SizeSrc);
     ParaDeal(ParaK, SizeSrc);
     for (Amount = 0; Amount < SizeSrc; Amount++, ParaResK++)
             *ParaResK = ParaBuffer(ParaK, Amount, SizeSrc);
     free(ParaK);
     return 0;
 }
 /***********************************************************************************
 ***********************************************************************************/
 int main(int argc, char* argv[])
 {
     int Amount;
	 clock_t StartTime, FinishTime;	//	声明两个时间变量
	 double DiffTime;
	 StartTime = clock();		//	开始计时

     double BufferX[1024], BufferY[1024], ParaK[6];	//	5次拟合, 一共6个系数(包含常数项)
	 //	读入要拟合的数据
     if (GetXY((const char*)"test.txt", (double*)BufferX, (double*)BufferY, &Amount))
          return 0;
	
	  printf("Amount: %d\n", Amount);
     Cal((const double*)BufferX, (const double*)BufferY, Amount, sizeof(ParaK) / sizeof(double), (double*)ParaK);
	 printf("拟合系数为:\n");
	 printf("按升序排列\n");
     for (Amount = 0; Amount < sizeof(ParaK) / sizeof(double); Amount++)
             printf("ParaK[%d] = %lf\r\n", Amount, ParaK[Amount]);

	 FinishTime = clock();	//	结束计时
	 DiffTime = FinishTime - StartTime;	//拟合时间
	 printf("拟合时间为: %lf\n", DiffTime);
     return 0;
 }
 

输出结果为:

Amount: 20

拟合系数为:
按升序排列
ParaK[0] = 0.999157
ParaK[1] = -1.450258
ParaK[2] = -0.529332
ParaK[3] = 0.236626
ParaK[4] = 6.725930
ParaK[5] = -18.544115
拟合时间为: 9.000000ms

Matlab代码:

其中test.txt文件为上面的x,y向量

load test.txt;
 x = test(:, 1)';
 y = test(:, 2)';
 para = polyfit(x, y, 5)
 
 figure,
 plot(x,y,'b.')
 hold on
 plot(x,y,'r-')

Matlab拟合系数,降序排列:

para =   0.9992   -1.4503   -0.5293    0.2366    6.7259  -18.5441

Matlab拟合曲线结果图:

 


 

  • 24
    点赞
  • 173
    收藏
    觉得还不错? 一键收藏
  • 14
    评论
GSL(GNU Scientific Library)提供了最小二乘曲线拟合的功能,通过该功能可以根据给定的数据点,找到最符合数据的曲线模型。 最小二乘曲线拟合是一种常用的数据处理方法,旨在通过多项式、指数函数或其他数学模型,找到一个最优的拟合曲线,以描述数据的趋势或规律。在实际应用中,最小二乘曲线拟合常用于数据分析、信号处理、图像处理等领域。 使用GSL进行最小二乘曲线拟合的步骤如下: 1. 导入GSL库并初始化拟合模型参数,例如选择多项式拟合模型的阶数。 2. 提供待拟合的数据点,包括横坐标和纵坐标。 3. 调用GSL提供的函数,传入数据点和拟合模型参数,进行曲线拟合计算。 4. 根据计算结果,得到最优的拟合曲线模型的参数,例如多项式的系数。 5. 根据得到的拟合曲线模型参数,可以进行预测或者进一步分析。 需要注意的是,使用最小二乘曲线拟合时,可能会遇到过拟合或欠拟合的问题。过拟合指的是拟合曲线过于复杂,过度拟合了数据的噪声;欠拟合则是拟合曲线过于简单,无法很好地描述数据的特征。为了避免这些问题,选择合适的拟合模型和拟合参数非常重要。 总之,GSL最小二乘曲线拟合功能提供了一种方便、快捷和可靠的方法,用于对给定数据点进行最优曲线拟合。这种拟合方法在科学研究和工程实践中具有广泛的应用前景。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 14
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值