簡單思路如下：
1，采用目標函數(shù)對多項式系數(shù)求偏導(dǎo)，得到最優(yōu)值條件，組成一個方程組；
2，方程組的解法采用行列式變換（兩次變換：普通行列式——三角行列式——對角行列式——求解），行列式的求解算法上優(yōu)化過一次了，目前還沒有更好的思路再優(yōu)化運算方法，限幅和精度準備再修改修改
目前存在的問題：
1，代碼還是太粗糙
2，數(shù)學(xué)原理可行，但是計算機運算有幅度溢出和精度問題，這方面欠考慮，導(dǎo)致高階大數(shù)據(jù)可能擬合失敗
下面附簡單注釋版代碼：
#include "stdafx.h"
#include "stdlib.h"
#include "math.h"
//把二維的數(shù)組與一維數(shù)組的轉(zhuǎn)換，也可以直接用二維數(shù)組，只是我的習(xí)慣是不用二維數(shù)組
#define ParaBuffer(Buffer,Row,Col) (*(Buffer + (Row) * (SizeSrc + 1) + (Col)))
//
/***********************************************************************************
從txt文件里讀取double型的X，Y數(shù)據(jù)
txt文件里的存儲格式為
X1 Y1
X2 Y2
X3 Y3
X4 Y4
X5 Y5
X6 Y6
X7 Y7
X8 Y8
函數(shù)返回X，Y，以及數(shù)據(jù)的數(shù)目（以組為單位）
***********************************************************************************/
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;
}
/***********************************************************************************
打印系數(shù)矩陣，只用于調(diào)試，不具備運算功能
對于一個N階擬合，它的系數(shù)矩陣大小是（N + 1）行（N + 2）列
double* Para：系數(shù)矩陣存儲地址
int SizeSrc：系數(shù)矩陣大小（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;
}
/***********************************************************************************
系數(shù)矩陣的限幅處理，防止它溢出，目前這個函數(shù)很不完善，并不能很好地解決這個問題
原理：矩陣解行列式，同一行乘以一個系數(shù)，行列式的解不變
當然，相對溢出問題，還有一個精度問題，也是同樣的思路，現(xiàn)在對于這兩塊的處理很不完善，有待優(yōu)化
以行為單位處理
***********************************************************************************/
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;
}
/***********************************************************************************
系數(shù)矩陣行列式變換
***********************************************************************************/
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;
}
/***********************************************************************************
系數(shù)矩陣行列式變換，與ParaPreDealA配合
完成第一次變換，變換成三角矩陣
***********************************************************************************/
static int ParaDealA(double* Para, int SizeSrc)
{
int i;
for (i = SizeSrc; i; i--)
if (ParaPreDealA(Para, SizeSrc, i))
return -1;
return 0;
}
/***********************************************************************************
系數(shù)矩陣行列式變換
***********************************************************************************/
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;
}
/***********************************************************************************
系數(shù)矩陣行列式變換，與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;
}
/***********************************************************************************
系數(shù)矩陣變換
***********************************************************************************/
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數(shù)據(jù)里面獲取系數(shù)矩陣
double* Para：系數(shù)矩陣地址
const double* X：X數(shù)據(jù)地址
const double* Y：Y數(shù)據(jù)地址
int Amount：XY數(shù)據(jù)組數(shù)
int SizeSrc：系數(shù)矩陣大�。⊿izeSrc）行（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函數(shù)
數(shù)據(jù)組數(shù)：20
階數(shù)：5
***********************************************************************************/
int main(int argc, char* argv[])
{
//數(shù)據(jù)組數(shù)
int Amount;
//XY緩存，系數(shù)矩陣緩存
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);
//輸出系數(shù)
for (Amount = 0; Amount < sizeof(ParaK) / sizeof(double); Amount++)
printf("ParaK[%d] = %lf!\r\n", Amount, ParaK[Amount]);
//屏幕暫留
system("pause");
return 0;
}
|
