/***********************************************************************************
打印系数矩阵,只用于调试,不具备运算功能
对于一个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;
}
二阶的系数矩阵方程组:
(目标函数y = A * x^2 + b * x^1 + C *x^0)
注意!
(X^N表示的是X1^N + X2^N + X3^N + X4^N + …… + Xn^N)
(Y * X^N表示的是Y1 * X1^N + Y1 * X2^N + Y3 * X3^N + Y4 * X4^N + …… + Yn * Xn^N)
X^4 * A + X^3 * B + X^2 * C = Y * X ^2
X^3 * A + X^2 * B + X^1 * C = Y * X ^1
X^2 * A + X^1 * B + X^0 * C = Y * X ^0
三阶的系数矩阵方程组:
注意!
(X^N表示的是X1^N + X2^N + X3^N + X4^N + …… + Xn^N)
(Y * X^N表示的是Y1 * X1^N + Y1 * X2^N + Y3 * X3^N + Y4 * X4^N + …… + Yn * Xn^N)
(目标函数y = A * x^3 + b * x^2 + C *x^1 + D * x^0)
X^6 * A + X^5 * B + X^4 * C + X^3 * D = Y * X ^3
X^5 * A + X^4 * B + X^3 * C + X^2 * D = Y * X ^2
X^4 * A + X^3 * B + X^2 * C + X^1 * D = Y * X ^1
X^3 * A + X^2 * B + X^1 * C + X^0 * D = Y * X ^0
更高阶的规律也跟上面一样,所以第一步就是获取系数矩阵,需要算的值有X^0, X^1, …… X^(2N), Y * X^0, Y * X^1,……Y* X^N,N表示阶数
在本程序中用如下代码实现:
/***********************************************************************************
最小二乘法的第一步就是从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;
//先把第0行的搞定
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);
//把第(SizeSrc - 1)列的搞定
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);
//把第(SizeSrc)列的搞定
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;
}
/***********************************************************************************
系数矩阵行列式变换,与ParaPreDealA配合
完成第一次变换,变换成三角矩阵
***********************************************************************************/
static int ParaDealA(double* Para, int SizeSrc)
{
int i;
for (i = SizeSrc; i; i--)
if (ParaPreDealA(Para, SizeSrc, i))
return -1;
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;
}