最小二乘曲线拟合——C语言算法实现一
来源:互联网 发布:淘宝摄影师 编辑:程序博客网 时间:2024/06/10 19:53
最小二乘曲线拟合
给定一组数据,我们要对这组数据进行曲线拟合。
假定要拟合的曲线方程为 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.000000msMatlab代码:
其中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拟合曲线结果图:
- 最小二乘曲线拟合——C语言算法实现一
- 最小二乘曲线拟合——C语言算法实现二
- 最小二乘算法 C 语言实现
- 最小二乘曲线拟合matlab实现
- 最小二乘曲线拟合matlab实现
- 非线性曲线拟合函数 lsqcurvefit 最小二乘
- 最小二乘曲线拟合的MATLAB仿真
- 最小二乘C
- 最小二乘法曲线拟合 C语言实现
- CART—最小二乘回归树Python实现
- C语言实现矩阵连乘算法
- PLS 偏最小二乘算法 demo <一>
- 最小二乘算法之我见(一)
- 加权最小二乘滤波算法原理及实现
- C 实现最小二乘,步骤以及代码
- C语言——Prim算法实现最小生成树
- 算法导论—矩阵链乘C/C++实现
- 最小二乘回归树生成算法
- Java Swing界面编程(17)---单行文本输入组件:JTextField
- poj 2049
- hibernate回顾之缓存机制-一级缓存、二级缓存、查询缓存
- GNU风格ARM汇编编程实战之一 <C与汇编混合编程>
- input标签只能输入正整数数字
- 最小二乘曲线拟合——C语言算法实现一
- STM32学习之路-感觉自己走到了一个天大的坑里了!
- document.body.scrollTop用法
- POJ 3373 Changing Digits(记录路径的dp)
- ACM中java的使用(速成)
- 自己的音乐播放器
- 行逻辑联接的顺序表实现矩阵转置及乘法
- OCP 1Z0 051 143
- 产品做出来了,我们该怎么办?