最小二乘法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;
}
自己顶一下 上传比较完整的分析文档
低阶拟合算法的优化:
当阶数比较低时,如直线或者抛物线,解行列式可以优化,不需要做成递推模式,直接强行解就是,而且一般不会出现溢出问题。
还有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;
}
上传拟合高阶曲线的VC6.0工程,经测试在20个数据,6阶的拟合中与matlab拟合结果基本吻合,没有出现数据溢出,由于取消了限幅函数,数据量大(线性增加),阶数高(觅增长),数据大(指数增长)的时候估计还是会溢出,欢迎测试 LZ有没有指数曲线拟合的C程序? 无级电工 发表于 2013-5-22 11:57 static/image/common/back.gif
LZ有没有指数曲线拟合的C程序?
目前没有,你要的话可以试着帮你推一个,这个多项式拟合也就前天才开始推导,但指数比这个难一点,需要比较严谨的数学方法证明才行 无级电工 发表于 2013-5-22 11:57 static/image/common/back.gif
LZ有没有指数曲线拟合的C程序?
如果你仔细看那个一阶或二阶的推导过程,应该你也能把指数的推出来,只不过把多项式换成指数函数,(最好经过严格的数学证明,或者有参考工具可以校验,比如matlab) zouzhichao 发表于 2013-5-22 12:09 static/image/common/back.gif
如果你仔细看那个一阶或二阶的推导过程,应该你也能把指数的推出来,只不过把多项式换成指数函数,(最好 ...
说错了,指数跟这个还是很不相同的,第一步无法得到一个确定的系数矩阵,所以得换思路才行 太感谢楼主了,雪中送炭。
不过楼主你能否弄个 递推抛物拟合。就是不断的有新数据进入,老数据推出;或者对老数据遗忘处理?
kmani 发表于 2013-5-22 12:20 static/image/common/back.gif
太感谢楼主了,雪中送炭。
不过楼主你能否弄个 递推抛物拟合。就是不断的有新数据进入,老数据推出;或者对 ...
不太明白你的意思,是动态的拟合吗?每次的XY数据有类似FIFO的更新? kmani 发表于 2013-5-22 12:20 static/image/common/back.gif
太感谢楼主了,雪中送炭。
不过楼主你能否弄个 递推抛物拟合。就是不断的有新数据进入,老数据推出;或者对 ...
抛物线拟合的代码和数学推导都比较简单,你试着改一改吧 zouzhichao 发表于 2013-5-22 12:54 static/image/common/back.gif
不太明白你的意思,是动态的拟合吗?每次的XY数据有类似FIFO的更新?
嗯,是的。
因为我想递推拟合的话,计算量可能会小一些。
直接拟合每次都重新计算了。 zouzhichao 发表于 2013-5-22 12:55 static/image/common/back.gif
抛物线拟合的代码和数学推导都比较简单,你试着改一改吧
直接拟合的我会,就是递推的不行,不太会推公式... 顶一个先 kmani 发表于 2013-5-22 12:56 static/image/common/back.gif
嗯,是的。
因为我想递推拟合的话,计算量可能会小一些。
直接拟合每次都重新计算了。 ...
你的意思是,假设这里有1000个数据,但是这1000个不是立马就有的,而是一边运行,一边产生的?
拟合1——100个数
拟合2——101个数
拟合3——102个数
拟合4——103个数
是这样吗?那么当前的拟合结果和上一次拟合结果需要关联吗? kmani 发表于 2013-5-22 12:58 static/image/common/back.gif
直接拟合的我会,就是递推的不行,不太会推公式...
明白了,你要拟合1000个数据,但一次只拟合100个,但是之前的拟合结果会影响这次的拟合结果,最后的效果是可以做到实时拟合,并且不需要计算1000个数的拟合,只需要计算100个数的拟合,是这样吧? 那这个的逻辑关系就复杂了,暂时没有好的思路 至少先数学推导才敢下笔 本帖最后由 kmani 于 2013-5-22 18:02 编辑
zouzhichao 发表于 2013-5-22 13:05 static/image/common/back.gif
明白了,你要拟合1000个数据,但一次只拟合100个,但是之前的拟合结果会影响这次的拟合结果,最后的效果 ...
嗯,就是这样的,每次我只用100个数据拟合。下一次计算时进来一个新的数据,同时去掉一个最旧的数据,再次拟合。
最终的目的就是实时拟合,新数据不断地进来,旧数据不断地出去,保证每次拟合都是100个数据。
我为啥想用递推算法呢,因为假如使用楼主的算法,每次都计算了100个点的拟合,下一次拟合还是计算100个点,就没有使用上一次计算的部分结果,这样增大了计算量。
就跟拉格朗日插值一样,只要新增加点之后,所有的项都要重新计算。但是牛顿插值就只用计算新的项就可以了,不用全部从头算。但是我的想法是有了新的数据后,还有自动去掉旧的数据。
假如用递推方式计算,我估计只能采用加遗忘因子的方法了,不能准确的使用100个数据。 加精吧。很好{:lol:} kmani 发表于 2013-5-22 17:45 static/image/common/back.gif
嗯,就是这样的,每次我只用100个数据拟合。下一次计算时进来一个新的数据,同时去掉一个最旧的数据,再 ...
要实现精确的数学模型来实现有点难度,可能模糊地实现(不保证数学上的最优点,但是更接近最优点)可能好实现一些 jacobson 发表于 2013-5-22 18:21 static/image/common/back.gif
加精吧。很好
多谢夸奖 kmani 发表于 2013-5-22 17:45 static/image/common/back.gif
嗯,就是这样的,每次我只用100个数据拟合。下一次计算时进来一个新的数据,同时去掉一个最旧的数据,再 ...
多交流哈,递推的拟合你要是有什么好的思路教教我,突然对这个很有兴趣 目前临近毕业了,杂七杂八的事情多了,沉不下心看书,等忙完了再思考思考上面提到的指数曲线拟合和递推拟合 厉害,向楼主学习 xiaowenshao 发表于 2013-5-22 23:22 static/image/common/back.gif
厉害,向楼主学习
{:smile:} 厉害,向楼主学习{:lol:}最佩服有知识的人了 收下 收藏 谢谢。 楼主很佩服你啊!学问做的这么好!还这么谦虚! jsszdfdn 发表于 2013-5-23 09:40 static/image/common/back.gif
楼主很佩服你啊!学问做的这么好!还这么谦虚!
学问真不好,这点实在惭愧,现在快毕业了还差十来个学分没修完{:sad:} ,这个也不是用了什么高深的算法,推导过程不难,只是要花点时间写草稿演算 做个标记,有时间研究 徐士良 常用算法程序集
里面都有
biansf2001 发表于 2013-5-23 09:59 static/image/common/back.gif
徐士良 常用算法程序集
里面都有
非常感谢 biansf2001 发表于 2013-5-23 09:59 static/image/common/back.gif
徐士良 常用算法程序集
里面都有
浏览了一下目录,好书! ma入库 好,有空看看{:lol:} 收下,不容错过!学习 wangfei1956 发表于 2013-5-29 10:14 static/image/common/back.gif
收下,不容错过!学习
谢谢支持 不错的帖子,好好学习下 无级电工 发表于 2013-5-22 11:57 static/image/common/back.gif
LZ有没有指数曲线拟合的C程序?
拟合指数函数的示例,是这样子的吗? zouzhichao 发表于 2013-6-3 08:29 static/image/common/back.gif
拟合指数函数的示例,是这样子的吗?
看起来效果不错啊。 无级电工 发表于 2013-6-3 08:33 static/image/common/back.gif
看起来效果不错啊。
用的迭代,有个问题,有时候不收敛,数据也不能太大 无级电工 发表于 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),进入迭代 zouzhichao 发表于 2013-6-3 08:46 static/image/common/back.gif
用的迭代,有个问题,有时候不收敛,数据也不能太大
这倒是个问题。不急,将来搞成熟了再发布也行啊。 网上指数函数拟合的C程序没见过,你这应该算头一份了。 本帖最后由 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);
测试图:
无级电工 发表于 2013-6-3 08:54 static/image/common/back.gif
这倒是个问题。不急,将来搞成熟了再发布也行啊。 网上指数函数拟合的C程序没见过,你这应该算头一份了 ...
目前还没有写C,但是没有用matlab现成的拟合函数,最高级的也只是用了一个矩阵除法(解行列式用),结合之前高阶多项式拟合的时候写的那个纯C的行列式求解代码,很容易纯C实现 好帖子,顶一下! gfy200866 发表于 2013-6-3 09:23 static/image/common/back.gif
好帖子,顶一下!
多谢帮顶 zouzhichao 发表于 2013-6-3 08:29 static/image/common/back.gif
拟合指数函数的示例,是这样子的吗?
这样的,曲线该怎么实现算法。不知楼主有没有好的原创案例。 gfy200866 发表于 2013-6-3 09:25 static/image/common/back.gif
这样的,曲线该怎么实现算法。不知楼主有没有好的原创案例。
?没听明白 如果不是一个标准的指数函数。但也接近。LZ有什么好的办法减少失真? gfy200866 发表于 2013-6-3 09:27 static/image/common/back.gif
如果不是一个标准的指数函数。但也接近。LZ有什么好的办法减少失真?
*的点不是指数函数了,加了随机偏差,拟合的话就是找一条标准的指数函数曲线,使之跟那些点尽可能的近,那条实线就是拟合出来的 zouzhichao 发表于 2013-6-3 09:32 static/image/common/back.gif
*的点不是指数函数了,加了随机偏差,拟合的话就是找一条标准的指数函数曲线,使之跟那些点尽可能的近, ...
LZ的意思是需要先找一个接近实际需要的标准函数,然后在想办法进行拟合吗? 这个就是所说的最小二乘法吗?以前没接触过这个算法。有没有相关的资料分享一下吗?谢谢! {:smile:}{:smile:}牛.... 本帖最后由 zouzhichao 于 2013-6-3 09:41 编辑
gfy200866 发表于 2013-6-3 09:35 static/image/common/back.gif
这个就是所说的最小二乘法吗?以前没接触过这个算法。有没有相关的资料分享一下吗?谢谢! ...
资料全部在上面,最后那一部分还是今天早上刚出炉的,别的资料我只能说大学课本+自己的思维了 支持一下~ x9fish 发表于 2013-6-3 09:46 static/image/common/back.gif
支持一下~
感谢支持 本帖最后由 zouzhichao 于 2013-6-3 11:08 编辑
gfy200866 发表于 2013-6-3 09:35 static/image/common/back.gif
这个就是所说的最小二乘法吗?以前没接触过这个算法。有没有相关的资料分享一下吗?谢谢! ...
楼上有网友推荐的那本算法书不错
33楼,徐士良 《常用算法程序集》 非常感谢楼主!入库收藏。 woshimajia222 发表于 2013-6-3 11:45 static/image/common/back.gif
非常感谢楼主!入库收藏。
还不太完善,有待改进 zouzhichao 发表于 2013-6-3 11:07 static/image/common/back.gif
楼上有网友推荐的那本算法书不错
33楼,徐士良 《常用算法程序集》
非常感谢!{:handshake:} 学习 谢谢
楼主太厉害了,虽然看不懂,先mark! 非常感谢。 谢谢分享,下载看看 MARK一下 COOL ! 收藏!~ 每次看算法都头疼。。。 这个帖子不错哦 本帖最后由 blue.fox 于 2013-9-20 04:56 编辑
LZ代码好多啊。
LZ这个是在整非线性最小二乘法的叠代吧?有在DSP上实际测试过执行时间没?
另外,LZ矩阵求逆在DSP里怎么实现的?没看到相关的部分
谢谢 mark{:smile:} blue.fox 发表于 2013-9-20 04:54 static/image/common/back.gif
LZ代码好多啊。
LZ这个是在整非线性最小二乘法的叠代吧?有在DSP上实际测试过执行时间没?
另外,LZ矩阵求 ...
是的,没有在DSP中跑,都是在GCC和VC6.0,WinAVR GCC和GCC中编译测试的,前面的发帖是多项式拟合,里面有矩阵求逆的代码(解一个方程组),后面的发帖是为了解决其中有一个坛友提出的指数曲线拟合的问题,用的是迭代的方法,测试了指数函数拟合以及正弦函数拟合,都还行,只是对迭代初始化位置比较敏感,尤其是正弦函数的频率参数非常敏感,改进方法暂时没有新的思路 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里面难以实现。最后换的别的方法。 好贴!!!{:victory:} matlab不行吗! 笑傲江湖 发表于 2013-10-16 19:54 static/image/common/back.gif
matlab不行吗!
从来就没说matlab不行啊 很厉害的感觉 强帖留名,刚好有需要,谢谢楼主分享. 楼主厉害,赞一个! 学习学习。。。 mark 一下 收藏成功 先赞一个,有需要时再来下载。 顶顶。。。。收藏了 不错,看哈 mark下先,以后可能用到 小白学习中……{:smile:} mark{:smile:} 晕死,matlab这些都还给老师了。 膜拜大神、、 mark,有空再看 湛泸骏驰 发表于 2014-1-21 22:23
膜拜大神、、
惭愧,非大神 绝对的好贴 make下,回去学习{:lol:} 写得好,值得保留 make下,回去学习 mark,good
标记mark