最小二乘法拟合求教

zouzhichao · 发表于 2013-6-3 00:48:25

本帖最后由 zouzhichao 于 2013-6-3 00:56 编辑

这几天在弄最小二乘法曲线拟合（目标函数不是多项式，而是正弦函数），遇到一个问题。
假设要拟合的目标函数为y = a * sin(b * x) + c，有没有什么好的思路？
（目前我用的迭代，但是很容易震荡，对参数的初始值有很大依赖性，尤其是参数b）

zouzhichao · 发表于 2013-6-3 00:58:51

%得到一阶导数的值，公式如下
%2 * sin(b * x) * (c - y + a *sin(b * x))
%2 * a * x * cos(b * x) * (c - y + a * sin(b * x))
%2 * c - 2 * y + 2 * a * sin(b * x)
function Res = GetYijieDaoshu(Xy, Para)
Res = zeros(3, 1);
for i = 1 : size(Xy, 1)
Res(1, 1) = Res(1, 1) + 2 * sin(Para(2) * Xy(i, 1)) * (Para(3) - Xy(i, 2) + Para(1) * sin(Para(2) * Xy(i, 1)));
Res(2, 1) = Res(2, 1) + 2 * Para(1) * Xy(i, 1) * cos(Para(2) * Xy(i, 1)) * (Para(3) - Xy(i, 2) + Para(1) * sin(Para(2) * Xy(i, 1)));
Res(3, 1) = Res(3, 1) + 2 * Para(3) - 2 * Xy(i, 2) + 2 * Para(1) * sin(Para(2) * Xy(i, 1));
end
end

%得到二阶导数的值，公式如下
%2 * sin(b * x)^2
%2 * x * cos(b * x) * (c - y + a * sin(b * x)) + 2 * a * x * cos(b * x) * sin(b * x)
%2 * sin(b * x)

%2 * x * cos(b * x) * (c - y + a * sin(b * x)) + 2 * a * x * cos(b * x) * sin(b * x)
%2 * a^2 * x^2 * cos(b * x)^2 - 2 * a * x^2 * sin(b * x) * (c - y + a * sin(b * x))
%2 * a * x * cos(b * x)

%2 * sin(b * x)
%2 * a * x * cos(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 * sin(Para(2) * Xy(i, 1))^2;
Res(1, 2) = Res(1, 2) + 2 * Xy(i, 1) * cos(Para(2) * Xy(i, 1)) * (Para(3) - Xy(i, 2) + Para(1) * sin(Para(2) * Xy(i, 1))) + 2 * Para(1) * Xy(i, 1) * cos(Para(2) * Xy(i, 1)) * sin(Para(2) * Xy(i, 1));
Res(1, 3) = Res(1, 3) + 2 * sin(Para(2) * Xy(i, 1));

Res(2, 1) = Res(2, 1) + 2 * Xy(i, 1) * cos(Para(2) * Xy(i, 1)) * (Para(3) - Xy(i, 2) + Para(1) * sin(Para(2) * Xy(i, 1))) + 2 * Para(1) * Xy(i, 1) * cos(Para(2) * Xy(i, 1)) * sin(Para(2) * Xy(i, 1));
Res(2, 2) = Res(2, 2) + 2 * Para(1)^2 * Xy(i, 1)^2 * cos(Para(2) * Xy(i, 1))^2 - 2 * Para(1) * Xy(i, 1)^2 * sin(Para(2) * Xy(i, 1)) * (Para(3) - Xy(i, 2) + Para(1) * sin(Para(2) * Xy(i, 1)));
Res(2, 3) = Res(2, 3) + 2 * Para(1) * Xy(i, 1) * cos(Para(2) * Xy(i, 1));

Res(3, 1) = Res(3, 1) + 2 * sin(Para(2) * Xy(i, 1));
Res(3, 2) = Res(3, 2) + 2 * Para(1) * Xy(i, 1) * cos(Para(2) * Xy(i, 1));
Res(3, 3) = Res(3, 3) + 2;
end
end

%计算系数的值
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

clear;
clc;
a = 1.7;
b = 3.3;
c = -4.6;
X = 200 : 0.01 : 203;
Y = a * sin(b * X) + c + (rand(size(X)) - 0.5) / 1.0;
ParaABC = CalPara([X', Y'], [20.91, 3.3, -34.6]', 40);
Z = ParaABC(1) * sin(ParaABC(2) * X) + ParaABC(3);
plot(X, Y, X, Z);

asma · 发表于 2013-6-4 09:40:40

拿低通滤波得了。

笑笑我笑了 · 发表于 2013-6-4 09:57:19

最小二乘貌似有两种推导方式,一种是通过偏导数的,一种是通过矩阵的.

zouzhichao · 发表于 2013-6-4 10:00:07

笑笑我笑了发表于 2013-6-4 09:57
最小二乘貌似有两种推导方式,一种是通过偏导数的,一种是通过矩阵的.

是的，我用的偏导，线性代数没学好，矩阵的推导还没完全搞通

zouzhichao · 发表于 2013-6-4 10:07:08

笑笑我笑了发表于 2013-6-4 09:57
最小二乘貌似有两种推导方式,一种是通过偏导数的,一种是通过矩阵的.

目标函数为多项式的时候容易一些，参数之间没有耦合。目标函数为别的的话，参数之间有耦合，不好弄，我现在用的迭代，测试了正（余）弦函数，指数函数，正（余）弦函数y = a * sin(b * x) + c（还没有把相位弄进去）的参数b的初始值很敏感。指数函数y = a * exp(b * x) + c稍微好一点，初始值的影响不是那么大，一般都能收敛。
不知道有没有好点的思路可以借鉴

笑笑我笑了 · 发表于 2013-6-4 11:56:39

zouzhichao 发表于 2013-6-4 10:07
目标函数为多项式的时候容易一些，参数之间没有耦合。目标函数为别的的话，参数之间有耦合，不好弄，我现 ...

觉得还是用矩阵好,线性代数有一套很完善的数值计算理论.

zouzhichao · 发表于 2013-6-4 14:06:41

笑笑我笑了发表于 2013-6-4 11:56
觉得还是用矩阵好,线性代数有一套很完善的数值计算理论.

求提点。。。

quackonchen · 发表于 2013-6-19 22:58:43

这个波形，不禁想起卡尔曼滤波。。。

zouzhichao · 发表于 2013-6-21 00:46:21

quackonchen 发表于 2013-6-19 22:58
这个波形，不禁想起卡尔曼滤波。。。

这个确实是采用“预测—实测—修正”的思路迭代的

最小二乘法拟合求教

本帖子中包含更多资源

阿莫论坛20周年了！感谢大家的支持与爱护！！