最小二乘法一阶线性拟合二阶曲线拟合的C语言程序实现_百

合集下载

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

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

最小二乘法曲线拟合C语言实现简单思路如下:1,采用目标函数对多项式系数求偏导,得到最优值条件,组成一个方程组;2,方程组的解法采用行列式变换(两次变换:普通行列式——三角行列式——对角行列式——求解),行列式的求解算法上优化过一次了,目前还没有更好的思路再优化运算方法,限幅和精度准备再修改修改目前存在的问题:1,代码还是太粗糙2,数学原理可行,但是计算机运算有幅度溢出和精度问题,这方面欠考虑,导致高阶大数据可能拟合失败下面附简单注释版代码:#include "stdafx.h"#include "stdlib.h"#include "math.h"//把二维的数组与一维数组的转换,也可以直接用二维数组,只是我的习惯是不用二维数组#define ParaBuffer(Buffer,Row,Col) (*(Buffer + (Row) * (SizeSrc + 1) + (Col)))///************************************************************** *********************从txt文件里读取double型的X,Y数据txt文件里的存储格式为X1 Y1X2 Y2X3 Y3X4 Y4X5 Y5X6 Y6X7 Y7X8 Y8函数返回X,Y,以及数据的数目(以组为单位)*************************************************************** ********************/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;}/************************************************************** *********************打印系数矩阵,只用于调试,不具备运算功能对于一个N阶拟合,它的系数矩阵大小是(N + 1)行(N + 2)列double* Para:系数矩阵存储地址int SizeSrc:系数矩阵大小(SizeSrc)行(SizeSrc + 1)列***********************************************************************************/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;}/************************************************************** *********************系数矩阵行列式变换,与ParaPreDealA配合完成第一次变换,变换成三角矩阵*************************************************************** ********************/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;}/************************************************************** *********************系数矩阵行列式变换,与ParaPreDealB配合完成第一次变换,变换成对角矩阵,变换完毕***********************************************************************************/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;}/************************************************************** *********************最小二乘法的第一步就是从XY数据里面获取系数矩阵double* Para:系数矩阵地址const double* X:X数据地址const double* Y:Y数据地址int Amount:XY数据组数int SizeSrc:系数矩阵大小(SizeSrc)行(SizeSrc + 1)列*************************************************************** ********************/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;}/************************************************************** *********************测试main函数数据组数:20阶数:5***********************************************************************************/int main(int argc, char* argv[]){//数据组数int Amount;//XY缓存,系数矩阵缓存double BufferX[1024], BufferY[1024], ParaK[6];if (GetXY((const char*)"test.txt", (double*)BufferX, (double*)BufferY, &Amount))return 0;//运算Cal((const double*)BufferX, (const double*)BufferY, Amount, sizeof(ParaK) / sizeof(double), (double*)ParaK);//输出系数for (Amount = 0; Amount < sizeof(ParaK) / sizeof(double); Amount++)printf("ParaK[%d] = %lf!\r\n", Amount, ParaK[Amount]);//屏幕暂留system("pause");return 0;}。

c语言 曲线拟合

c语言 曲线拟合

c语言曲线拟合曲线拟合(Curve Fitting)是数据处理的常用方法之一,其基本思想是通过已知的一组数据点,找到一条曲线,使得这条曲线尽可能地接近这些数据点。

在C语言中,可以使用最小二乘法进行曲线拟合。

以下是一个简单的C语言代码示例,用于实现二次多项式拟合:```c#include<stdio.h>#include<math.h>#define N5//数据点个数int main(){double x[N]={1,2,3,4,5};//自变量数据点double y[N]={2.2,2.8,3.6,4.5,5.1};//因变量数据点double a[3]={0,0,0};//二次多项式系数,初始化为0double sum=0,sumx=0,sumx2=0,sumxy= 0;int i;for(i=0;i<N;i++){sum+=y[i];sumx+=x[i];sumx2+=x[i]*x[i];sumxy+=x[i]*y[i];}double mean_y=sum/N;//计算y的平均值double mean_x=sumx/N;//计算x的平均值//计算二次多项式系数a[0]=(N*sumxy-sumx*sumy)/(N*sumx2 -sumx*sumx);a[1]=(mean_y-a[0]*mean_x)/N;a[2]=mean_y-a[0]*mean_x-a[1];printf("拟合曲线为:y=%.2fx^2+%.2fx+%.2f\n", a[0],a[1],a[2]);return0;}```在这个示例中,我们首先定义了5个数据点,然后使用最小二乘法计算了二次多项式的系数。

最后,我们输出了拟合曲线的公式。

用C语言实现的曲线拟合的最小二乘法

用C语言实现的曲线拟合的最小二乘法

实验名称:曲线拟合的最小二乘法 实验目的了解曲线拟合的最小二乘法实验类型设计型实验环境Windows XP TC实验内容相关知识:已知C [a,b ]中函数f (x )的一组实验数据(x i ,y i )(i=0,1,…,m),其中y i =f (x i )。

设);,,1,0)((m n n j x j <= ϕ是C [a ,b]上线性无关函数族。

在)}(,),(),({10x x x span n ϕϕϕφ =中找函数f(x) 曲线拟合的最小二乘解∑==nj j j x a x S 0*)()(ϕ,其法方程(组)为:),,1,0(),(0n k d a nj k j j k==∑=ϕϕ其中,∑==mi i k i j i k j x x x 0)()()(),(ϕϕωϕϕ∑=≡=mi k i k i i k d x x f x f 0)()()(),(ϕωϕ k=0,1,…,n特别是,求函数f(x ) 曲线拟合的线性最小二乘解b ax x S +=)(*的计算公式为:∑∑∑∑∑∑======-+-=mi i m i i mi i i mi i mi i mi i x x m y x x y x b 02202)()1())(())((∑∑∑∑∑=====-+-+=mi i m i i mi i m i i m i i i x x m y x y x m a 0220)()1())(()1(数据结构:两个一维数组或一个二维数组 算法设计:(略)实验用例: 已知函数y=f (x)的一张表:处的近似值,并画出实验数据和直线。

编写代码:#include〈stdio.h>#include<stdlib.h>#include<graphics.h〉double qiuhe1(double a[10][2],int p){int i;double y;y=0;for(i=0;i〈10;i++)y=y+a[i][p];return y;}double qiuhe2(double a[10][2],int p){int i;double y=0;for(i=0;i<10;i++)y=y+a[i][0]*a[i][p];return y;}double nihe(double a[10][2],double x){double a1,b,y;a1=(10*qiuhe2(a,1)-qiuhe1(a,0)*qiuhe1(a,1))/(10*qiuhe2(a,0)-qiuhe1(a,0)*qiuhe1(a,0));b=(qiuhe2(a,0)*qiuhe1(a,1)—qiuhe1(a,0)*qiuhe2(a,1))/(10*qiuhe2(a,0)—qiuhe1(a,0)*qiuhe1(a,0));y=a1*x+b;return y;}int main(){double a[10][2]={0,68,10,67.1,20,66.4,30,65.6,40,64.6,50,61.8, 60,61.0,70,60.8,80,60.4,90,60};double x,x1,q=1;c har c[12];i nt i;long n;int arw[6]={515,235,520,240,515,245};int arw1[6]={315,45,320,40,325,45};int gdriver=IBM8514;int gmode=IBM8514HI;initgraph(&gdriver, &gmode, ”c:\\TC20\\BGI");cleardevice();printf(”input x:\n");scanf("%lf",&x);printf("%f\n",nihe(a,x));n=nihe(a,x)*1000000+1;c[0]='y’;c[1]=’=';c[4]='。

C语言最小二乘法

C语言最小二乘法

函数逼近与曲线拟合,用最小二乘法进行曲线拟合的C或C++编写的完整程序!已知x 0 5 10 15 20 25 30 35 40 45 50 55y 0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64近似解析表达式为y=at+bt^2+ct^3求a,b,c曲线拟合:#include <stdio.h>#include <stdlib.h>#include <malloc.h>#include <math.h>Smooth(double *x,double *y,double *a,int n,int m,double *dt1,double *dt2,double*dt3);void main(){int i ,n ,m ;double *x,*y,*a,dt1,dt2,dt3,b;n = 12;// 12个样点m = 4; //3次多项式拟合b = 0; //x的初值为0/*分别为x,y,a分配存贮空间*/x = (double *)calloc(n,sizeof(double));if(x == NULL){printf("内存分配失败\n");exit (0);}y = (double *)calloc(n,sizeof(double));if(y == NULL){printf("内存分配失败\n");exit (0);}a = (double *)calloc(n,sizeof(double));if(a == NULL){printf("内存分配失败\n");exit (0);}for(i=1;i<=n;i++){x[i-1]=b+(i-1)*5;/*每隔5取一个点,这样连续取12个点*/}y[0]=0;y[1]=1.27;y[2]=2.16;y[3]=2.86;y[4]=3.44;y[5]=3.87;y[6]=4.15;y[7]=4.37;y[8]=4.51;y[9]=4.58;y[10]=4.02;y[11]=4.64;/*x[i-1]点对应的y值是拟合已知值*/Smooth(x,y,a,n,m,&dt1,&dt2,&dt3); /*调用拟合函数*/for(i=1;i<=m;i++)printf("a[%d] = %.10f\n",(i-1),a[i-1]);printf("拟合多项式与数据点偏差的平方和为:\n");printf("%.10e\n",dt1);printf("拟合多项式与数据点偏差的绝对值之和为:\n");printf("%.10e\n",dt2);printf("拟合多项式与数据点偏差的绝对值最大值为:\n");printf("%.10e\n",dt3);free(x); /*释放存储空间*/free(y); /*释放存储空间*/free(a); /*释放存储空间*/}Smooth(double *x,double *y,double *a,int n,int m,double *dt1,double *dt2,double*dt3)//(x,y,a,n,m,dt1,dt2,dt3 )//double *x; /*实型一维数组,输入参数,存放节点的xi值*///double *y; /*实型一维数组,输入参数,存放节点的yi值*///double *a; /*双精度实型一维数组,长度为m。

最小二乘法 c语言

最小二乘法 c语言

最小二乘法 c语言最小二乘法是一种常用的数学方法,用于通过已知数据点拟合出一条最佳拟合曲线。

在本文中,我们将讨论如何使用C语言实现最小二乘法。

我们需要明确最小二乘法的基本原理。

最小二乘法的目标是找到一条曲线,使得该曲线上的点到已知数据点的距离之和最小。

具体地,我们假设已知数据点的集合为{(x1, y1), (x2, y2), ..., (xn, yn)},我们需要找到一条曲线y = f(x),使得f(xi)与yi的差的平方和最小。

那么,如何在C语言中实现最小二乘法呢?首先,我们需要定义一个函数来计算拟合曲线上的点f(xi)。

在这个函数中,我们可以使用多项式的形式来表示拟合曲线。

例如,我们可以选择使用一次多项式y = ax + b来拟合数据。

然后,我们可以使用最小二乘法的公式来计算出最优的a和b的值。

接下来,我们需要编写一个函数来计算拟合曲线上每个点f(xi)与已知数据点yi的差的平方和。

通过遍历已知数据点的集合,并计算每个点的差的平方,最后将所有差的平方求和,即可得到拟合曲线的误差。

然后,我们可以使用梯度下降法来最小化误差函数。

梯度下降法是一种优化算法,通过不断迭代来找到误差函数的最小值。

在每次迭代中,我们通过计算误差函数对参数a和b的偏导数,来更新a和b的值。

通过多次迭代,最终可以找到最优的a和b的值,从而得到最佳拟合曲线。

我们可以编写一个主函数来调用以上的函数,并将最终的拟合曲线绘制出来。

在这个主函数中,我们可以读取已知数据点的集合,并调用最小二乘法函数来计算拟合曲线的参数。

然后,我们可以使用绘图库来绘制已知数据点和拟合曲线,并将结果输出到屏幕上。

通过以上的步骤,我们就可以使用C语言实现最小二乘法了。

当然,在实际应用中,我们可能会遇到更复杂的数据和更高阶的多项式拟合。

但是基本的原理和方法是相同的,只是需要做一些适当的调整。

总结一下,最小二乘法是一种常用的数学方法,用于通过已知数据点拟合出一条最佳拟合曲线。

最小二乘法C语言编程

最小二乘法C语言编程
最小二乘法c语言编程最小二乘法c语言程序c语言最小二乘法最小二乘法c语言代码c语言编程九九乘法表最小二乘法最小二乘法公式最小二乘法曲线拟合最小二乘法原理偏最小二乘法
#include <hidef.h> /* common defines and macros */
#include "derivative.h" /* derivative-specific definitions */
x2=x*x*x*x;
x3+=x2;
}
*(pp++)=x3;
}
}
pp=*(b+1);
p=a9;
ppp=a8;
x3=0;
for(i=0;i<10;i++) {
x=*(p++);
x2=*(ppp++);
x3+=x*x2;
}
*pp=x3;
pp=*(b+2);
p=a9;
ppp=a8;
x3=0;
for(i=0;i<10;i++) {
*(p++)=(x1*x5-x2*x4);
p=*(b1+0);
pp=*(b4+0);
for (i=0;i<3;i++) {
for(j=0;j<3;j++){
x1=*(p++);
x=x1*c;
*(pp++)=x;
}
}
//////////////////////////求两矩阵的商//////////////////////////////////////

最小二乘法一阶线性拟合二阶曲线拟合的C语言程序实现

最小二乘法一阶线性拟合二阶曲线拟合的C语言程序实现
m n i
当 n=1 时,为 1 阶拟合,又称直线拟合,即系数矩阵是一个 2*2 的矩阵,通过线性方程的求解运 算,求得线性回归方程的系数表达式为:
当 n=2 时,为 2 阶曲线拟合,所得到的系数矩阵是一个 3*3 的矩阵【用 aij(i,j=1,2,……)的 形式表达】 ,通过线性方程的求解运算,求得线性回归方程的系数表达式为:
二、1 阶 2 阶拟合功能子函数和计算表达式
通过分析以上系数计算式中各项计算式,写出全部需要用到的子函数:
通过对照系数表达式里各个项的计算表达,写入主函数进行拟合计算。设定输入的数据格式为(x[ i ],y[ i ]) ,用 户输入数据的个数为 c,计算表达式程序代码如下: 1 阶直线拟合:
2 阶曲线拟合:
一、最小二乘法原理与计算方法
对于 m 组测量数据,选取
( x) a0 a1 x a2 x 2
an x n 进行 n 阶拟合,按
照残差平方和最小原则,对各个待定系数求偏导数,使之都等于 0,通过数学运算可得到各个系数的 最小二乘的法方程为 运算式如下,求解这个线性方程组就可以得出各个系数的值。
三、主函数代码
四、用 MATLAB 验证程序的运行结果
第一组:选择 y=x+1 进行线性拟合检验,可见 2 阶拟合对于线性关系,二次项系数为 0
第二组:选择 y=x^2+1 进行 2 阶曲线拟合检验
第三组:进行常规数据组检验
数据输入部分的设计参考了[物理实验计算器.
Zhouzb .
zhouzb889@]的部分代码,在此表示感谢。
智能仪器设计作业——最小二乘法——高世浩 1223150078
m 1 i 1 m xi i 1 m x n i i 1

c语言 最小二乘法

c语言 最小二乘法

c语言最小二乘法最小二乘法是一种数据拟合方法,可以用来找到最优解的参数。

在 C 语言中,可以使用矩阵运算和线性代数的方法来实现最小二乘法。

首先,需要准备好数据集。

假设有一组数据集 (x1, y1), (x2, y2), ..., (xn, yn),我们要拟合的模型是 y = a*x + b。

这个模型可以写成矩阵形式为 Y = X*P,其中 Y 是一个 n*1 的列矩阵,X 是一个 n*2 的矩阵,P 是一个 2*1 的列矩阵,表示模型的参数 a 和 b。

接下来,可以使用矩阵运算来求解 P。

具体地,可以通过求解 X^T * X 的逆矩阵,再乘以 X^T 和 Y,得到 P = (X^T * X)^-1 * X^T * Y。

实现代码如下:```#include <stdio.h>#include <stdlib.h>#include <math.h>#define N 10 // 数据集中的数据个数#define M 2 // 模型中的参数个数// 数据集double x[N] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};double y[N] = {2.1, 4.5, 7.4, 9.5, 12.1, 14.5, 17.3, 19.5, 22.2, 24.5};int main(){// 构造矩阵 X 和 Ydouble X[N][M] = {0};double Y[N][1] = {0};for (int i = 0; i < N; i++){X[i][0] = x[i];X[i][1] = 1;Y[i][0] = y[i];}// 求解 Pdouble XtX[M][M] = {0};double XtY[M][1] = {0};for (int i = 0; i < N; i++){for (int j = 0; j < M; j++){for (int k = 0; k < N; k++){XtX[j][k] += X[k][j] * X[k][i]; }XtY[j][0] += X[i][j] * Y[i][0];}}// 求解 XtX 的逆矩阵double det = XtX[0][0] * XtX[1][1] - XtX[0][1] * XtX[1][0]; double invXtX[M][M] = {{XtX[1][1] / det, -XtX[0][1] / det},{-XtX[1][0] / det, XtX[0][0] / det}};// 计算 Pdouble P[M][1] = {0};for (int i = 0; i < M; i++){for (int j = 0; j < M; j++){P[i][0] += invXtX[i][j] * XtY[j][0];}}// 输出结果printf('a = %lf, b = %lf', P[0][0], P[1][0]);return 0;}```输出结果为:```a = 2.340000,b = -0.300000```表示拟合得到的模型为 y = 2.34*x - 0.3。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
相关文档
最新文档