基于最小二乘法的曲线拟合的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  Y1
X2  Y2
X3  Y3
X4  Y4
X5  Y5
X6  Y6
X7  Y7
X8  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;
}

展开阅读全文

没有更多推荐了,返回首页