zouzhichao 发表于 2013-5-22 11:40:35

最小二乘法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 = {0.0};

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

/*******************************************************************************
系数矩阵的运算,没有撒谎,真的只有7行
*******************************************************************************/
static int ParaDeal(void)
{
        Para -= Para * (Para / Para);
        Para -= Para * (Para / Para);
        Para = 0;
        Para -= Para * (Para / Para);
        Para = 0;
        Para /= Para;
        //Para = 1.0;
        Para /= Para;
        //Para = 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, BufferY;
        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, Para);
        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 = {0.0};

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

/*******************************************************************************
系数矩阵的运算
*******************************************************************************/
static int ParaDeal(void)
{
        //用Para把Para消成0
        Para -= Para * (Para / Para);
        Para -= Para * (Para / Para);
        Para -= Para * (Para / Para);
        Para = 0;
        //用Para把Para 消成0
        Para -= Para * (Para / Para);
        Para -= Para * (Para / Para);
        Para -= Para * (Para / Para);
        Para = 0;
        //用Para把Para消成0
        Para -= Para * (Para / Para);
        Para -= Para * (Para / Para );
        Para = 0;
        //至此,已经成成为三角矩阵
        //用Para把Para消成0
        Para -= Para * (Para / Para);
        Para = 0;
        //用Para把Para消成0
        Para -= Para * (Para / Para);
        Para = 0;
        //用Para把Para消成0
        Para -= Para * (Para / Para);
        Para = 0;
        //至此,已经成成为对角矩阵
        //把Para,Para,Para消成1
        Para /= Para;//这个就是最后a的值
        //Para = 1.0;
        Para /= Para; //这个就是最后b的值
        //Para = 1.0;
        Para /= Para; //这个就是最后c的值
        //Para = 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, BufferY;
        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, Para, Para);
        system("pause");
        return 0;
}

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

zouzhichao 发表于 2013-5-22 11:41:10

自己顶一下

zouzhichao 发表于 2013-5-22 11:43:12

上传比较完整的分析文档
低阶拟合算法的优化:
当阶数比较低时,如直线或者抛物线,解行列式可以优化,不需要做成递推模式,直接强行解就是,而且一般不会出现溢出问题。
还有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)
(1)首先是定义一个2行3列的系数矩阵,并初始化
double Para = {0.0};
int ParaInit(const double* X, const double* Y, int Amount)
{
        Para = Amount;
        for ( ; Amount; Amount--, X++, Y++)
{
        Para+=my_pow2(*X);
Para+=my_pow1(*X);
//Para+=my_pow0(*Y);
Para+=my_pow1(*X) * my_pow1(*Y);
//Para+=my_pow0(*X) * my_pow1(*Y);
Para+=my_pow1(*Y);
}
Para = Para;
return 0;
}

(2)解行列式:
由方程组
a * ∑x^2 + b * ∑x^1 = ∑(y * x^1)
a * ∑x^1 + b * ∑x^0 = ∑(y * x^0)
首先可以看到:
A,Para和Para肯定大于0
B,Para = Para
C,        可以证明行列式的值
Para * Para - Para * Para
(∑x^2) * (∑x^0) - (∑x^1) * (∑x^1)> 0
也即是说,方阵正定,偏导为零的点就是最优值的点
可以看出,第1行的系数比第0行的系数要小很多,所以先采用第1行消第0行,再采用第0行消第1行。
int ParaDeal(void)
{
/*先把Para消成0,由于有Para肯定大于0,所以用乘法(除法),除法更好一些,数据溢出的问题几乎不存在,适应性会更好*/
Para = Para - Para * (Para / Para;
Para = Para - Para * (Para / Para;
Para = 0;
//再把Para 消成0
Para = Para - Para * (Para / Para);
Para = 0;
//把Para,Para 消成1
Para /= Para;//这个就是最后a的值
Para = 1.0;
Para /= Para; //这个就是最后b的值
Para = 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)
(1)首先是定义一个3行4列的系数矩阵,并初始化
double Para = {0.0};
int ParaInit(const double* X, const double* Y, int Amount)
{
        Para = Amount;
        for ( ; Amount; Amount--, X++, Y++)
{
Para+=my_pow4(*X);
Para+= my_pow3(*X);
Para+= my_pow2(*X);
Para+= my_pow1(*X);
//Para +=my_pow0(*X);
Para+=my_pow2(*X) * my_pow1(*Y);
Para+=my_pow1(*X) * my_pow1(*Y);
Para+=my_pow0(*X) * my_pow1(*Y);
}
Para = Para;
Para = Para;
Para = Para;
Para = Para;
return 0;
}

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)
首先可以看到:
A,Para,Para, Para肯定大于0
B,左边3*3的方阵是对称矩阵
D,        可以证明行列式的值 > 0,也即是说,方阵正定,偏导为零的点就是最优值的点
可以看出,第k + 1行的系数比第k行的系数要小很多,所以先采用第k +1行消第k行,再采用第k行消第k + 1行。

int ParaDeal(void)
{
//用Para把Para消成0
Para = Para - Para * (Para / Para);
Para = Para - Para * (Para / Para);
Para = Para - Para * (Para / Para);
Para = 0;
//用Para把Para 消成0
Para = Para - Para * (Para / Para);
Para = Para - Para * (Para / Para);
Para = Para - Para * (Para / Para);
Para = 0;
//用Para把Para消成0
Para = Para - Para * (Para / Para);
Para = Para - Para * (Para / Para );
Para = 0;
//至此,已经成成为三角矩阵
//用Para把Para消成0
Para = Para - Para * (Para / Para);
Para = 0;
//用Para把Para消成0
Para = Para - Para * (Para / Para);
Para = 0;
//用Para把Para消成0
Para = Para - Para * (Para / Para);
Para = 0;
//至此,已经成成为对角矩阵
//把Para,Para,Para消成1
Para /= Para;//这个就是最后a的值
Para = 1.0;
Para /= Para; //这个就是最后b的值
Para = 1.0;
Para /= Para; //这个就是最后c的值
Para = 1.0;
return 0;
}

zouzhichao 发表于 2013-5-22 11:54:03


上传拟合高阶曲线的VC6.0工程,经测试在20个数据,6阶的拟合中与matlab拟合结果基本吻合,没有出现数据溢出,由于取消了限幅函数,数据量大(线性增加),阶数高(觅增长),数据大(指数增长)的时候估计还是会溢出,欢迎测试

无级电工 发表于 2013-5-22 11:57:04

LZ有没有指数曲线拟合的C程序?

zouzhichao 发表于 2013-5-22 11:59:41

无级电工 发表于 2013-5-22 11:57 static/image/common/back.gif
LZ有没有指数曲线拟合的C程序?

目前没有,你要的话可以试着帮你推一个,这个多项式拟合也就前天才开始推导,但指数比这个难一点,需要比较严谨的数学方法证明才行

zouzhichao 发表于 2013-5-22 12:09:49

无级电工 发表于 2013-5-22 11:57 static/image/common/back.gif
LZ有没有指数曲线拟合的C程序?

如果你仔细看那个一阶或二阶的推导过程,应该你也能把指数的推出来,只不过把多项式换成指数函数,(最好经过严格的数学证明,或者有参考工具可以校验,比如matlab)

zouzhichao 发表于 2013-5-22 12:13:32

zouzhichao 发表于 2013-5-22 12:09 static/image/common/back.gif
如果你仔细看那个一阶或二阶的推导过程,应该你也能把指数的推出来,只不过把多项式换成指数函数,(最好 ...

说错了,指数跟这个还是很不相同的,第一步无法得到一个确定的系数矩阵,所以得换思路才行

kmani 发表于 2013-5-22 12:20:19

太感谢楼主了,雪中送炭。
不过楼主你能否弄个 递推抛物拟合。就是不断的有新数据进入,老数据推出;或者对老数据遗忘处理?

zouzhichao 发表于 2013-5-22 12:54:18

kmani 发表于 2013-5-22 12:20 static/image/common/back.gif
太感谢楼主了,雪中送炭。
不过楼主你能否弄个 递推抛物拟合。就是不断的有新数据进入,老数据推出;或者对 ...

不太明白你的意思,是动态的拟合吗?每次的XY数据有类似FIFO的更新?

zouzhichao 发表于 2013-5-22 12:55:45

kmani 发表于 2013-5-22 12:20 static/image/common/back.gif
太感谢楼主了,雪中送炭。
不过楼主你能否弄个 递推抛物拟合。就是不断的有新数据进入,老数据推出;或者对 ...

抛物线拟合的代码和数学推导都比较简单,你试着改一改吧

kmani 发表于 2013-5-22 12:56:31

zouzhichao 发表于 2013-5-22 12:54 static/image/common/back.gif
不太明白你的意思,是动态的拟合吗?每次的XY数据有类似FIFO的更新?

嗯,是的。
因为我想递推拟合的话,计算量可能会小一些。
直接拟合每次都重新计算了。

kmani 发表于 2013-5-22 12:58:06

zouzhichao 发表于 2013-5-22 12:55 static/image/common/back.gif
抛物线拟合的代码和数学推导都比较简单,你试着改一改吧

直接拟合的我会,就是递推的不行,不太会推公式...

pengshipower 发表于 2013-5-22 13:01:16

顶一个先

zouzhichao 发表于 2013-5-22 13:02:35

kmani 发表于 2013-5-22 12:56 static/image/common/back.gif
嗯,是的。
因为我想递推拟合的话,计算量可能会小一些。
直接拟合每次都重新计算了。 ...

你的意思是,假设这里有1000个数据,但是这1000个不是立马就有的,而是一边运行,一边产生的?
拟合1——100个数
拟合2——101个数
拟合3——102个数
拟合4——103个数

是这样吗?那么当前的拟合结果和上一次拟合结果需要关联吗?

zouzhichao 发表于 2013-5-22 13:05:36

kmani 发表于 2013-5-22 12:58 static/image/common/back.gif
直接拟合的我会,就是递推的不行,不太会推公式...

明白了,你要拟合1000个数据,但一次只拟合100个,但是之前的拟合结果会影响这次的拟合结果,最后的效果是可以做到实时拟合,并且不需要计算1000个数的拟合,只需要计算100个数的拟合,是这样吧?

zouzhichao 发表于 2013-5-22 13:06:40

那这个的逻辑关系就复杂了,暂时没有好的思路

zouzhichao 发表于 2013-5-22 13:07:12

至少先数学推导才敢下笔

kmani 发表于 2013-5-22 17:45:08

本帖最后由 kmani 于 2013-5-22 18:02 编辑

zouzhichao 发表于 2013-5-22 13:05 static/image/common/back.gif
明白了,你要拟合1000个数据,但一次只拟合100个,但是之前的拟合结果会影响这次的拟合结果,最后的效果 ...

嗯,就是这样的,每次我只用100个数据拟合。下一次计算时进来一个新的数据,同时去掉一个最旧的数据,再次拟合。
最终的目的就是实时拟合,新数据不断地进来,旧数据不断地出去,保证每次拟合都是100个数据。
我为啥想用递推算法呢,因为假如使用楼主的算法,每次都计算了100个点的拟合,下一次拟合还是计算100个点,就没有使用上一次计算的部分结果,这样增大了计算量。
就跟拉格朗日插值一样,只要新增加点之后,所有的项都要重新计算。但是牛顿插值就只用计算新的项就可以了,不用全部从头算。但是我的想法是有了新的数据后,还有自动去掉旧的数据。
假如用递推方式计算,我估计只能采用加遗忘因子的方法了,不能准确的使用100个数据。

jacobson 发表于 2013-5-22 18:21:29

加精吧。很好{:lol:}

zouzhichao 发表于 2013-5-22 18:25:15

kmani 发表于 2013-5-22 17:45 static/image/common/back.gif
嗯,就是这样的,每次我只用100个数据拟合。下一次计算时进来一个新的数据,同时去掉一个最旧的数据,再 ...

要实现精确的数学模型来实现有点难度,可能模糊地实现(不保证数学上的最优点,但是更接近最优点)可能好实现一些

zouzhichao 发表于 2013-5-22 18:25:38

jacobson 发表于 2013-5-22 18:21 static/image/common/back.gif
加精吧。很好

多谢夸奖

zouzhichao 发表于 2013-5-22 18:34:15

kmani 发表于 2013-5-22 17:45 static/image/common/back.gif
嗯,就是这样的,每次我只用100个数据拟合。下一次计算时进来一个新的数据,同时去掉一个最旧的数据,再 ...

多交流哈,递推的拟合你要是有什么好的思路教教我,突然对这个很有兴趣

zouzhichao 发表于 2013-5-22 18:38:30

目前临近毕业了,杂七杂八的事情多了,沉不下心看书,等忙完了再思考思考上面提到的指数曲线拟合和递推拟合

xiaowenshao 发表于 2013-5-22 23:22:26

厉害,向楼主学习

zouzhichao 发表于 2013-5-22 23:26:48

xiaowenshao 发表于 2013-5-22 23:22 static/image/common/back.gif
厉害,向楼主学习

{:smile:}

fjourdev 发表于 2013-5-23 07:41:35

厉害,向楼主学习{:lol:}最佩服有知识的人了

Eric2013 发表于 2013-5-23 08:29:35

收下      

wugang_1213 发表于 2013-5-23 08:31:57

收藏 谢谢。

jsszdfdn 发表于 2013-5-23 09:40:16

楼主很佩服你啊!学问做的这么好!还这么谦虚!

zouzhichao 发表于 2013-5-23 09:46:03

jsszdfdn 发表于 2013-5-23 09:40 static/image/common/back.gif
楼主很佩服你啊!学问做的这么好!还这么谦虚!

学问真不好,这点实在惭愧,现在快毕业了还差十来个学分没修完{:sad:} ,这个也不是用了什么高深的算法,推导过程不难,只是要花点时间写草稿演算

mcu5i51 发表于 2013-5-23 09:53:25

做个标记,有时间研究

biansf2001 发表于 2013-5-23 09:59:08

徐士良 常用算法程序集
里面都有

zouzhichao 发表于 2013-5-23 10:10:29

biansf2001 发表于 2013-5-23 09:59 static/image/common/back.gif
徐士良 常用算法程序集
里面都有

非常感谢

zouzhichao 发表于 2013-5-23 10:13:08

biansf2001 发表于 2013-5-23 09:59 static/image/common/back.gif
徐士良 常用算法程序集
里面都有

浏览了一下目录,好书!

ersha4877 发表于 2013-5-23 11:17:17

ma入库 好,有空看看{:lol:}

wangfei1956 发表于 2013-5-29 10:14:07

收下,不容错过!学习

zouzhichao 发表于 2013-5-29 10:22:03

wangfei1956 发表于 2013-5-29 10:14 static/image/common/back.gif
收下,不容错过!学习

谢谢支持

mypc16888 发表于 2013-5-29 11:01:46

不错的帖子,好好学习下

zouzhichao 发表于 2013-6-3 08:29:29

无级电工 发表于 2013-5-22 11:57 static/image/common/back.gif
LZ有没有指数曲线拟合的C程序?


拟合指数函数的示例,是这样子的吗?

无级电工 发表于 2013-6-3 08:33:27

zouzhichao 发表于 2013-6-3 08:29 static/image/common/back.gif
拟合指数函数的示例,是这样子的吗?

看起来效果不错啊。

zouzhichao 发表于 2013-6-3 08:46:27

无级电工 发表于 2013-6-3 08:33 static/image/common/back.gif
看起来效果不错啊。

用的迭代,有个问题,有时候不收敛,数据也不能太大

zouzhichao 发表于 2013-6-3 08:53:24

无级电工 发表于 2013-6-3 08:33 static/image/common/back.gif
看起来效果不错啊。

方法还是一阶偏导为零,但是没有一个确定的系数矩阵,不能直接求解行列式
解决办法是:增加一个二阶偏导(一个3*3的矩阵)
(1)给a,b,c一组初始值,求出一阶偏导的值,二阶偏导矩阵
(2)一阶偏导 / 二阶偏导矩阵,得到da,db,dc。(求解行列式)
(3)a = a - da; b = b - db; c = c - dc; 更新a,b,c的值,回到(1),进入迭代

无级电工 发表于 2013-6-3 08:54:07

zouzhichao 发表于 2013-6-3 08:46 static/image/common/back.gif
用的迭代,有个问题,有时候不收敛,数据也不能太大

这倒是个问题。不急,将来搞成熟了再发布也行啊。 网上指数函数拟合的C程序没见过,你这应该算头一份了。

zouzhichao 发表于 2013-6-3 09:09:18

本帖最后由 zouzhichao 于 2013-6-3 09:13 编辑

无级电工 发表于 2013-6-3 08:33 static/image/common/back.gif
看起来效果不错啊。

m代码分四个:
(1)求解一阶导数GetYijieDaoshu.m
%得到一阶导数的值,公式如下
%2 * exp(b * x) * (c - y + a * exp(b * x))
%2 * a * x * exp(b * x) * (c - y + a * exp(b * x))
%2 * c - 2 * y + 2 * a * exp(b * x)
function Res = GetYijieDaoshu(Xy, Para)
Res = zeros(3, 1);
for i = 1 : size(Xy, 1)
    Res(1, 1) = Res(1, 1) + 2 * exp(Para(2) * Xy(i, 1)) * (Para(3) - Xy(i, 2) + Para(1) * exp(Para(2) * Xy(i, 1)));
    Res(2, 1) = Res(2, 1) + 2 * Para(1) * Xy(i, 1) * exp(Para(2) * Xy(i, 1)) * (Para(3) - Xy(i, 2) + Para(1) * exp(Para(2) * Xy(i, 1)));
    Res(3, 1) = Res(3, 1) + 2 * Para(3) - 2 * Xy(i, 2) + 2 * Para(1) * exp(Para(2) * Xy(i, 1));
end
end

(2)求解二阶导数GetErjieDaoshu.m
%得到二阶导数的值,公式如下
%2 * exp(2 * b * x)
%2 * x * exp(b * x) * (c - y + a * exp(b * x)) + 2 * a * x * exp(2 * b * x)
%2 * exp(b * x)
%2 * x * exp(b * x) * (c - y + a * exp(b * x)) + 2 * a * x * exp(2 * b * x)
%2 * a^2 * x^2 * exp(2 * b * x) + 2 * a * x^2 * exp(b * x) * (c - y + a * exp(b * x))
%2 * a * x * exp(b * x)
%2 * exp(b * x)
%2 * a * x * exp(b * x)
%2
function Res = GetErjieDaoshu(Xy, Para)
Res = zeros(3, 3);
for i = 1 : size(Xy, 1)
    Res(1, 1) = Res(1, 1) + 2 * exp(2 * Para(2) * Xy(i, 1));
    Res(1, 2) = Res(1, 2) + 2 * Xy(i, 1) * exp(Para(2) * Xy(i, 1)) * (Para(3) - Xy(i, 2) + Para(1) * exp(Para(2) * Xy(i, 1))) + 2 * Para(1) * Xy(i, 1) * exp(2 * Para(2) * Xy(i, 1));
    Res(1, 3) = Res(1, 3) + 2 * exp(Para(2) * Xy(i, 1));
    Res(2, 1) = Res(2, 1) + 2 * Xy(i, 1) * exp(Para(2) * Xy(i, 1)) * (Para(3) - Xy(i, 2) + Para(1) * exp(Para(2) * Xy(i, 1))) + 2 * Para(1) * Xy(i, 1) * exp(2 * Para(2) * Xy(i, 1));
    Res(2, 2) = Res(2, 2) + 2 * Para(1)^2 * Xy(i, 1)^2 * exp(2 * Para(2) * Xy(i, 1)) + 2 * Para(1) * Xy(i, 1)^2 * exp(Para(2) * Xy(i, 1)) * (Para(3) - Xy(i, 2) + Para(1) * exp(Para(2) * Xy(i, 1)));
    Res(2, 3) = Res(2, 3) + 2 * Para(1) * Xy(i, 1) * exp(Para(2) * Xy(i, 1));
    Res(3, 1) = Res(3, 1) + 2 * exp(Para(2) * Xy(i, 1));
    Res(3, 2) = Res(3, 2) + 2 * Para(1) * Xy(i, 1) * exp(Para(2) * Xy(i, 1));
    Res(3, 3) = Res(3, 3) + 2;
end
end

(3)迭代计算CalPara.m
%计算系数的值
function Res = CalPara(Xy, ParaInit, N)
Res = ParaInit;
for i = 1 : N
    ParaN = GetYijieDaoshu(Xy, Res);
    ParaNN = GetErjieDaoshu(Xy, Res);
    Res = Res - 1 * (ParaNN \ ParaN);
end
end

(4)测试Test.m
clear;
clc;
%目标函数y = a * exp(b * x) + c
%预定参数
a = 1.7;
b = 0.2;
c = -4.6;
%待拟合的点,带上(rand(size(X)) - 0.5) * 0.2的随机偏差
X = -3 : 0.1 : 2.9;
Y = a * exp(b * X) + c + (rand(size(X)) - 0.5) * 0.2;
%拟合,注意需要给出a, b, c的初始值,不要离得太远,否则会发散导致迭代失败
%这里我们取',跟预定参数'距离不是很远,迭代200次
ParaABC = CalPara(, ', 200);
%绘图比较
X1 = min(X) : (max(X)- min(X))/9999 : max(X);
Z = ParaABC(1) * exp(ParaABC(2) * X1) + ParaABC(3);
plot(X, Y, '*', X1, Z);

测试图:


zouzhichao 发表于 2013-6-3 09:16:59

无级电工 发表于 2013-6-3 08:54 static/image/common/back.gif
这倒是个问题。不急,将来搞成熟了再发布也行啊。 网上指数函数拟合的C程序没见过,你这应该算头一份了 ...

目前还没有写C,但是没有用matlab现成的拟合函数,最高级的也只是用了一个矩阵除法(解行列式用),结合之前高阶多项式拟合的时候写的那个纯C的行列式求解代码,很容易纯C实现

gfy200866 发表于 2013-6-3 09:23:37

好帖子,顶一下!

zouzhichao 发表于 2013-6-3 09:24:26

gfy200866 发表于 2013-6-3 09:23 static/image/common/back.gif
好帖子,顶一下!

多谢帮顶

gfy200866 发表于 2013-6-3 09:25:05

zouzhichao 发表于 2013-6-3 08:29 static/image/common/back.gif
拟合指数函数的示例,是这样子的吗?

这样的,曲线该怎么实现算法。不知楼主有没有好的原创案例。

zouzhichao 发表于 2013-6-3 09:27:23

gfy200866 发表于 2013-6-3 09:25 static/image/common/back.gif
这样的,曲线该怎么实现算法。不知楼主有没有好的原创案例。

?没听明白

gfy200866 发表于 2013-6-3 09:27:33

如果不是一个标准的指数函数。但也接近。LZ有什么好的办法减少失真?

zouzhichao 发表于 2013-6-3 09:32:13

gfy200866 发表于 2013-6-3 09:27 static/image/common/back.gif
如果不是一个标准的指数函数。但也接近。LZ有什么好的办法减少失真?

*的点不是指数函数了,加了随机偏差,拟合的话就是找一条标准的指数函数曲线,使之跟那些点尽可能的近,那条实线就是拟合出来的

gfy200866 发表于 2013-6-3 09:34:27

zouzhichao 发表于 2013-6-3 09:32 static/image/common/back.gif
*的点不是指数函数了,加了随机偏差,拟合的话就是找一条标准的指数函数曲线,使之跟那些点尽可能的近, ...

LZ的意思是需要先找一个接近实际需要的标准函数,然后在想办法进行拟合吗?

gfy200866 发表于 2013-6-3 09:35:42

这个就是所说的最小二乘法吗?以前没接触过这个算法。有没有相关的资料分享一下吗?谢谢!

xtx8962 发表于 2013-6-3 09:36:05

{:smile:}{:smile:}牛....

zouzhichao 发表于 2013-6-3 09:39:21

本帖最后由 zouzhichao 于 2013-6-3 09:41 编辑

gfy200866 发表于 2013-6-3 09:35 static/image/common/back.gif
这个就是所说的最小二乘法吗?以前没接触过这个算法。有没有相关的资料分享一下吗?谢谢! ...

资料全部在上面,最后那一部分还是今天早上刚出炉的,别的资料我只能说大学课本+自己的思维了

x9fish 发表于 2013-6-3 09:46:01

支持一下~

zouzhichao 发表于 2013-6-3 09:47:39

x9fish 发表于 2013-6-3 09:46 static/image/common/back.gif
支持一下~

感谢支持

zouzhichao 发表于 2013-6-3 11:07:02

本帖最后由 zouzhichao 于 2013-6-3 11:08 编辑

gfy200866 发表于 2013-6-3 09:35 static/image/common/back.gif
这个就是所说的最小二乘法吗?以前没接触过这个算法。有没有相关的资料分享一下吗?谢谢! ...

楼上有网友推荐的那本算法书不错
33楼,徐士良 《常用算法程序集》

woshimajia222 发表于 2013-6-3 11:45:59

非常感谢楼主!入库收藏。

zouzhichao 发表于 2013-6-3 11:53:48

woshimajia222 发表于 2013-6-3 11:45 static/image/common/back.gif
非常感谢楼主!入库收藏。

还不太完善,有待改进

gfy200866 发表于 2013-6-3 15:06:11

zouzhichao 发表于 2013-6-3 11:07 static/image/common/back.gif
楼上有网友推荐的那本算法书不错
33楼,徐士良 《常用算法程序集》

非常感谢!{:handshake:}

jordonwu 发表于 2013-6-3 15:18:27

学习 谢谢

绿茶山人 发表于 2013-6-13 18:08:34

楼主太厉害了,虽然看不懂,先mark!

sunwu 发表于 2013-6-14 15:58:10

非常感谢。

liuanty7 发表于 2013-8-29 08:56:05

谢谢分享,下载看看

hephi 发表于 2013-8-31 13:39:55

MARK一下

armok 发表于 2013-9-19 15:09:45

COOL !

ljt80158015 发表于 2013-9-19 15:13:15

收藏!~                                 

szmini2006 发表于 2013-9-19 15:25:56

每次看算法都头疼。。。

凌滨 发表于 2013-9-19 15:40:45

这个帖子不错哦

blue.fox 发表于 2013-9-20 04:54:01

本帖最后由 blue.fox 于 2013-9-20 04:56 编辑

LZ代码好多啊。
LZ这个是在整非线性最小二乘法的叠代吧?有在DSP上实际测试过执行时间没?
另外,LZ矩阵求逆在DSP里怎么实现的?没看到相关的部分
谢谢

kneken 发表于 2013-9-20 08:03:23

mark{:smile:}

zouzhichao 发表于 2013-9-23 12:47:14

blue.fox 发表于 2013-9-20 04:54 static/image/common/back.gif
LZ代码好多啊。
LZ这个是在整非线性最小二乘法的叠代吧?有在DSP上实际测试过执行时间没?
另外,LZ矩阵求 ...

是的,没有在DSP中跑,都是在GCC和VC6.0,WinAVR GCC和GCC中编译测试的,前面的发帖是多项式拟合,里面有矩阵求逆的代码(解一个方程组),后面的发帖是为了解决其中有一个坛友提出的指数曲线拟合的问题,用的是迭代的方法,测试了指数函数拟合以及正弦函数拟合,都还行,只是对迭代初始化位置比较敏感,尤其是正弦函数的频率参数非常敏感,改进方法暂时没有新的思路

blue.fox 发表于 2013-9-23 12:55:49

zouzhichao 发表于 2013-9-23 12:47 static/image/common/back.gif
是的,没有在DSP中跑,都是在GCC和VC6.0,WinAVR GCC和GCC中编译测试的,前面的发帖是多项式拟合,里面有 ...

如果是多项式线性拟合的话,(A^T * A)^-1 * A^T部分是个常数矩阵,可以先用matlab算好,然后C语言就只需要做一次矩阵乘法了

指数曲线的迭代我也分析过。里面的求逆,判断奇异的运算量太大了。而且的确对初始位置比较敏感,DSP里面难以实现。最后换的别的方法。

追寻cheney 发表于 2013-10-2 15:47:31

好贴!!!{:victory:}

笑傲江湖 发表于 2013-10-16 19:54:01

matlab不行吗!                              

zouzhichao 发表于 2013-10-16 22:18:56

笑傲江湖 发表于 2013-10-16 19:54 static/image/common/back.gif
matlab不行吗!

从来就没说matlab不行啊

四轴飞行器 发表于 2013-10-17 07:38:36

很厉害的感觉

zhiqianlin 发表于 2014-1-13 12:07:12

强帖留名,刚好有需要,谢谢楼主分享.

ethan_free 发表于 2014-1-13 14:32:28

楼主厉害,赞一个!

moouse 发表于 2014-1-21 00:23:05

学习学习。。。

maimaige 发表于 2014-1-21 08:57:35

mark 一下

luhuizszw 发表于 2014-1-21 10:43:40

收藏成功

Robin_King 发表于 2014-1-21 12:10:29

先赞一个,有需要时再来下载。

棋间卒 发表于 2014-1-21 12:50:37

顶顶。。。。收藏了

cwei 发表于 2014-1-21 12:53:25

不错,看哈

zhudadragon 发表于 2014-1-21 13:52:07

mark下先,以后可能用到

不进则退 发表于 2014-1-21 14:34:30

小白学习中……{:smile:}

bisoo 发表于 2014-1-21 15:55:04

mark{:smile:}

active12 发表于 2014-1-21 21:38:23

晕死,matlab这些都还给老师了。

湛泸骏驰 发表于 2014-1-21 22:23:23

膜拜大神、、

孤独_求败 发表于 2014-1-22 08:37:58

mark,有空再看

zouzhichao 发表于 2014-1-22 12:36:01

湛泸骏驰 发表于 2014-1-21 22:23
膜拜大神、、

惭愧,非大神

xml2028 发表于 2014-1-22 12:43:33

绝对的好贴

HEYsir 发表于 2014-1-22 14:14:37

make下,回去学习{:lol:}

farmerzhangdl 发表于 2014-1-23 09:31:55

写得好,值得保留

redwolf310 发表于 2014-1-24 13:14:57

make下,回去学习

xz516@163.com 发表于 2014-1-24 15:06:54

mark,good

残忆视觉 发表于 2014-2-18 10:05:18

标记mark
页: [1] 2 3
查看完整版本: 最小二乘法C算法终极整理版本,绝对原创!