最小二乘法C算法终极整理版本

天道出勤

      前两天发了一个关于最小二乘法C语言算法的帖子（http://www.amobbs.com/thread-5535084-1-1.html），这次重新发一个，但是绝对不是属于重复发帖，如果你对算法感兴趣，这个绝对适合你的口味，这次综合了最小二乘法高阶函数拟合（对之前存在的溢出问题有了比较好的解决办法，没有彻底解决），直线拟合（经过高度优化，初始化代码6行，运算代码7行！），抛物线拟合（经过高度优化，初始化代码12行，运算代码20行！）先上代码，工程在附件里面，后续会有进一步的分析（目前还有一个最小二乘法的算法正在研究，思路跟这个版本相差很大，尚不知道是否能实现）



最

小二

乘法

拟合直线（包含测试函数）：
#include "stdafx.h"
#include "process.h"

/*******************************************************************************
阶乘函数宏定义
*******************************************************************************/
#define     my_pow2(x)      ((x) * (x))

/*******************************************************************************
定义一个2行3列的系数矩阵
*******************************************************************************/
static double Para[2][3] = {0.0};

/*******************************************************************************
系数矩阵的初始化，没有撒谎，真的只有6行 + 一个循环
*******************************************************************************/
static int ParaInit(double* X, double* Y, int Amount)
{
        Para[1][1] = Amount;
        for ( ; Amount; Amount--, X++, Y++)
        {
                Para[0][0] += my_pow2(*X);
                Para[0][1] += (*X);
                Para[0][2] += (*X) * (*Y);
                Para[1][2] += (*Y);
        }
        Para[1][0] = Para[0][1];
        return 0;
}

/*******************************************************************************
系数矩阵的运算，没有撒谎，真的只有7行
*******************************************************************************/
static int ParaDeal(void)
{
        Para[0][0] -= Para[1][0] * (Para[0][1] / Para[1][1]);
        Para[0][2] -= Para[1][2] * (Para[0][1] / Para[1][1]);
        Para[0][1] = 0;
        Para[1][2] -= Para[0][2] * (Para[1][0] / Para[0][0]);
        Para[1][0] = 0;
        Para[0][2] /= Para[0][0];
        //Para[0][0] = 1.0;
        Para[1][2] /= Para[1][1];
        //Para[1][1] = 1.0;
        return 0;
}

/***********************************************************************************
从txt文件里读取double型的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;
}

/*******************************************************************************
*******************************************************************************/
int Test(const char* FileName)
{
        int Amount;
        double BufferX[1024], BufferY[1024];
        if (GetXY(FileName, (double*)BufferX, (double*)BufferY, &Amount))
        {
                printf("读取数据文件失败！\r\n");
                return -1;
        }
        printf("读取数据文件成功，共有%d个点！\r\n", Amount);
        ParaInit((double*)BufferX, (double*)BufferY, Amount);
        ParaDeal();
        printf("拟合数据成功，拟合直线为：\r\ny = (%lf) * x + (%lf);\r\n", Para[0][2], Para[1][2]);
        system("pause");
        return 0;
}

int main(int argc, char* argv[])
{
        Test((const char*)"test.txt");
        return 0;
}

拟合抛物线（包含测试函数）：
#include "stdafx.h"
#include "process.h"

/*******************************************************************************
阶乘函数宏定义
*******************************************************************************/
#define     my_pow2(x)      ((x) * (x))
#define     my_pow3(x)      ((x) * my_pow2(x))
#define     my_pow4(x)      (my_pow2(x) * my_pow2(x))

/*******************************************************************************
定义一个2行3列的系数矩阵
*******************************************************************************/
static double Para[3][4] = {0.0};

/*******************************************************************************
系数矩阵的初始化
*******************************************************************************/
static int ParaInit(const double* X, const double* Y, int Amount)
{
        Para[2][2] = Amount;
        for ( ; Amount; Amount--, X++, Y++)
        {
                Para[0][0]  +=  my_pow4(*X);
                Para[0][1]  += my_pow3(*X);
                Para[0][2]  += my_pow2(*X);
                Para[1][2]  += (*X);
                Para[0][3]  +=  my_pow2(*X) * (*Y);
                Para[1][3]  +=  (*X) * (*Y);
                Para[2][3]  +=  (*Y);
        }
        Para[1][0] = Para[0][1];
        Para[1][1] = Para[0][2];
        Para[2][0] = Para[1][1];
        Para[2][1] = Para[1][2];
        return 0;
}

/*******************************************************************************
系数矩阵的运算
*******************************************************************************/
static int ParaDeal(void)
{
        //用Para[2][2]把Para[0][2]消成0
        Para[0][0] -= Para[2][0] * (Para[0][2] / Para[2][2]);
        Para[0][1] -= Para[2][1] * (Para[0][2] / Para[2][2]);
        Para[0][3] -= Para[2][3] * (Para[0][2] / Para[2][2]);
        Para[0][2] = 0;
        //用Para[2][2]把Para[1][2] 消成0
        Para[1][0] -= Para[2][0] * (Para[1][2] / Para[2][2]);
        Para[1][1] -= Para[2][1] * (Para[1][2] / Para[2][2]);
        Para[1][3] -= Para[2][3] * (Para[1][2] / Para[2][2]);
        Para[1][2] = 0;
        //用Para[1][1]把Para[0][1]消成0
        Para[0][0] -= Para[1][0] * (Para[0][1] / Para[1][1]);
        Para[0][3] -= Para[1][3] * (Para[0][1] / Para[1][1] );
        Para[0][1] = 0;
        //至此，已经成成为三角矩阵
        //用Para[0][0]把Para[1][0]消成0
        Para[1][3] -= Para[0][3] * (Para[1][0] / Para[0][0]);
        Para[1][0] = 0;
        //用Para[0][0]把Para[2][0]消成0
        Para[2][3] -= Para[0][3] * (Para[2][0] / Para[0][0]);
        Para[2][0] = 0;
        //用Para[1][1]把Para[2][1]消成0
        Para[2][3] -= Para[1][3] * (Para[2][1] / Para[1][1]);
        Para[2][1] = 0;
        //至此，已经成成为对角矩阵
        //把Para[0][0]，Para[1][1]，Para[2][2]消成1
        Para[0][3] /= Para[0][0];//这个就是最后a的值
        //Para[0][0] = 1.0;
        Para[1][3] /= Para[1][1]; //这个就是最后b的值
        //Para[1][1] = 1.0;
        Para[2][3] /= Para[2][2]; //这个就是最后c的值
        //Para[2][2] = 1.0;
        return 0;
}

/***********************************************************************************
从txt文件里读取double型的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;
}

/*******************************************************************************
*******************************************************************************/
int Test(const char* FileName)
{
        int Amount;
        double BufferX[1024], BufferY[1024];
        if (GetXY(FileName, (double*)BufferX, (double*)BufferY, &Amount))
        {
                printf("读取数据文件失败！\r\n");
                return -1;
        }
        printf("读取数据文件成功，共有%d个点！\r\n", Amount);
        ParaInit((double*)BufferX, (double*)BufferY, Amount);
        ParaDeal();
        printf("拟合数据成功，拟合直线为：\r\ny = (%lf) * x^2 + (%lf) * x + (%lf);\r\n", Para[0][3], Para[1][3], Para[2][3]);
        system("pause");
        return 0;
}

int main(int argc, char* argv[])
{
        Test((const char*)"test.txt");
        return 0;
}












本帖转载于他站


仅供同学学习用