C语言写的算法 最小二乘法

[复制链接]
194|22
 楼主 | 2017-12-27 15:37 | 显示全部楼层 |阅读模式
拟合直线(包含测试函数):
#include "stdafx.h"
#include "process.h"
 楼主 | 2017-12-27 15:37 | 显示全部楼层
/*******************************************************************************
阶乘函数宏定义
*******************************************************************************/
#define     my_pow2(x)      ((x) * (x))
 楼主 | 2017-12-27 15:38 | 显示全部楼层
/*******************************************************************************
定义一个2行3列的系数矩阵
*******************************************************************************/
static double Para[2][3] = {0.0};
 楼主 | 2017-12-27 15:38 | 显示全部楼层
/*******************************************************************************
系数矩阵的初始化,没有撒谎,真的只有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;
}
 楼主 | 2017-12-27 15:39 | 显示全部楼层
/*******************************************************************************
系数矩阵的运算,没有撒谎,真的只有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;
}
 楼主 | 2017-12-27 15:39 | 显示全部楼层
/***********************************************************************************
从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;
}
 楼主 | 2017-12-27 15:39 | 显示全部楼层
/*******************************************************************************
*******************************************************************************/
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;
}
 楼主 | 2017-12-27 15:39 | 显示全部楼层
int main(int argc, char* argv[])
{
        Test((const char*)"test.txt");
        return 0;
}
 楼主 | 2017-12-27 15:40 | 显示全部楼层
拟合抛物线(包含测试函数):
#include "stdafx.h"
#include "process.h"
 楼主 | 2017-12-27 15:40 | 显示全部楼层
/*******************************************************************************
阶乘函数宏定义
*******************************************************************************/
#define     my_pow2(x)      ((x) * (x))
#define     my_pow3(x)      ((x) * my_pow2(x))
#define     my_pow4(x)      (my_pow2(x) * my_pow2(x))
 楼主 | 2017-12-27 15:40 | 显示全部楼层
/*******************************************************************************
定义一个2行3列的系数矩阵
*******************************************************************************/
static double Para[3][4] = {0.0};
 楼主 | 2017-12-27 15:41 | 显示全部楼层
/*******************************************************************************
系数矩阵的初始化
*******************************************************************************/
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;
}
 楼主 | 2017-12-27 15:41 | 显示全部楼层
/*******************************************************************************
系数矩阵的运算
*******************************************************************************/
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;
}
 楼主 | 2017-12-27 15:41 | 显示全部楼层
/***********************************************************************************
从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;
}
 楼主 | 2017-12-27 15:43 | 显示全部楼层
/*******************************************************************************
*******************************************************************************/
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;
}
 楼主 | 2017-12-27 15:43 | 显示全部楼层
int main(int argc, char* argv[])
{
        Test((const char*)"test.txt");
        return 0;
}
 楼主 | 2017-12-27 15:43 | 显示全部楼层
低阶拟合算法的优化:
当阶数比较低时,如直线或者抛物线,解行列式可以优化,不需要做成递推模式,直接强行解就是,而且一般不会出现溢出问题。
还有pow(double x, double y)函数应该自己写替换,原因是在我们的应用中y肯定是整数,而这个函数按照double运算的,效率应该会差很多,所以来个阉割版的pow(double x, int y).在结束低的时候,直接写两个宏定义。
#define     my_pow0(x)      (1.0)
#define     my_pow1(x)      (x)
#define     my_pow2(x)      (my_pow1(x) * my_pow1(x))
#define     my_pow2(x)      (my_pow1(x) * my_pow2(x))
#define     my_pow4(x)      (my_pow2(x) * my_pow2(x))
#define     my_pow5(x)      (my_pow2(x) * my_pow3(x))
#define     my_pow6(x)      (my_pow3(x) * my_pow3(x))
直线拟合的简化:
假设最后拟合直线:
y = a * x + b
残差方程:
(a * x + b – y)^2 = x^2 * a^2 + b^2 + y^2 + 2 * (x * a * b – x * y * a – y * b)
残差方程对a偏导:
x^2 * a + x^1 * b – y * x^1
残差方程对a偏导:
x^1 * a + x^0 * b – y * x^0
当残差方程之和取极值时,偏导为零。
a * ∑x^2 + b * ∑x^1 = ∑(y * x^1)
a * ∑x^1 + b * ∑x^0 = ∑(y * x^0)
 楼主 | 2017-12-27 15:44 | 显示全部楼层
(1)首先是定义一个2行3列的系数矩阵,并初始化
double Para[2] [3] = {0.0};
int ParaInit(const double* X, const double* Y, int Amount)
{
        Para[1][1] = Amount;
        for ( ; Amount; Amount--, X++, Y++)
{
        Para[0][0]  +=  my_pow2(*X);
Para[0][1]  +=  my_pow1(*X);
//Para[1][1]  +=  my_pow0(*Y);
Para[0][2]  +=  my_pow1(*X) * my_pow1(*Y);
//Para[1][2]  +=  my_pow0(*X) * my_pow1(*Y);
Para[1][2]  +=  my_pow1(*Y);
}
Para[1][0] = Para[0][1];
return 0;
}
 楼主 | 2017-12-27 15:44 | 显示全部楼层
(2)解行列式:
由方程组
a * ∑x^2 + b * ∑x^1 = ∑(y * x^1)
a * ∑x^1 + b * ∑x^0 = ∑(y * x^0)
首先可以看到:
A,Para[0][0]和Para[1][1]肯定大于0
B,Para[0][1] = Para[1][0]
C,        可以证明[2][2]行列式的值
Para[0][0] * Para[1][1] - Para[0][1] * Para[1][0]
(∑x^2) * (∑x^0) - (∑x^1) * (∑x^1)> 0
也即是说,方阵正定,偏导为零的点就是最优值的点
可以看出,第1行的系数比第0行的系数要小很多,所以先采用第1行消第0行,再采用第0行消第1行。
int ParaDeal(void)
{
/*先把Para[0][1]消成0,由于有Para[1][1]肯定大于0,所以用乘法(除法),除法更好一些,数据溢出的问题几乎不存在,适应性会更好*/
Para[0][0] = Para[0][0] - Para[1][0] * (Para[0][1] / Para[1][1];
Para[0][2] = Para[0][2] - Para[1][2] * (Para[0][1] / Para[1][1];
Para[0][1] = 0;
//再把Para[1][0] 消成0
Para[1][2] = Para[1][2] - Para[0][2] * (Para[1][0] / Para[0][0]);
Para[1][0] = 0;
//把Para[0][0],Para[1][1] 消成1
Para[0][2] /= Para[0][0];//这个就是最后a的值
Para[0][0] = 1.0;
Para[1][2] /= Para[1][1]; //这个就是最后b的值
Para[1][1] = 1.0;
return 0;
}
至此,一阶直线拟合的推导过程和C代码全部完毕,可以看出运算量可以说很小,两个函数均只需调用一次即可得到结果,在源数据的个数和大小不是特别大的时候,完全可以移植到MCU中进行运算,实用价值比较高。
再来个二阶抛物线拟合的优化:
前期推导和一阶直线差不多,略去,系数矩阵方程组:
a * ∑x^4 + b * ∑x^3 + c * ∑x^2 = ∑(y * x^2)
a * ∑x^3 + b * ∑x^2 + c * ∑x^1 = ∑(y * x^1)
a * ∑x^2 + b * ∑x^1 + c * ∑x^0 = ∑(y * x^0)
 楼主 | 2017-12-27 15:44 | 显示全部楼层
(1)首先是定义一个3行4列的系数矩阵,并初始化
double Para[3] [4] = {0.0};
int ParaInit(const double* X, const double* Y, int Amount)
{
        Para[1][1] = 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]  += my_pow1(*X);
//Para[2][2] +=  my_pow0(*X);
Para[0][3]  +=  my_pow2(*X) * my_pow1(*Y);
Para[1][3]  +=  my_pow1(*X) * my_pow1(*Y);
Para[2][3]  +=  my_pow0(*X) * my_pow1(*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;
}
扫描二维码,随时随地手机跟帖
您需要登录后才可以回帖 登录 | 注册 手机登录

本版积分规则

快速回复

您需要登录后才可以回帖
登录 | 注册 手机登录
高级模式

论坛热帖

关闭

热门推荐上一条 /4 下一条

分享 快速回复 返回顶部 返回列表