|
纯手写,刚写完,还没来得及写注释,已通过Matlab的polyfit验证(阶数高或者数据量太大会有double数据溢出的危险,低阶的都吻合),时间有点紧,程序注释,数学推导等时间充裕一点再贴上来
// 最小二乘法拟合.cpp : Defines the entry point for the console application.
//
#include "stdafx.h"
#include "stdlib.h"
#include "math.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;
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]);
return 0;
}
测试test.txt里的内容
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
matlab拟合m代码:
load test.txt;
x = test(:, 1)';
y = test(:, 2)';
para = polyfit(x, y, 5); |
阿莫论坛20周年了!感谢大家的支持与爱护!!
你熬了10碗粥,别人一桶水倒进去,淘走90碗,剩下10碗给你,你看似没亏,其实你那10碗已经没有之前的裹腹了,人家的一桶水换90碗,继续卖。说白了,通货膨胀就是,你的钱是挣来的,他的钱是印来的,掺和在一起,你的钱就贬值了。
|