我在谭里看了高斯牛顿法,也编了代码,发现两个问题1、数据过多,超过一个周期,就会发散(没想到解决办法)。2、初值要求较高,需要在精确值附近(可能可以用随机初值解决?)。
可是我需要只靠数据就得到a,b,c参数,我该怎么改,或着重新找什么方法去拟合。我在cftool上拟合sum of sine方法表现很好,我需要原理,参考代码,网上实在没找到。有路过大神还请赐教
高斯牛顿代码:
x1 = 1:7;
x1 = x1';
TargetA = 2pi/7;
TargetB = -pi/8;
y1 = sin(TargetAx1+TargetB);
plot(x1,y1);
lamda = 0.3;
a = 0.5;
b = 0.5;
disp(sprintf("original function f(x) = sin((%.2f)*x+(%.2f)) ",TargetA,TargetB));
disp("data:")
disp(x1');
disp(y1');
for k = 1:1:1000
fy = [];
JacobiMatrix = ones(length(x1),2);
w=[a,b];
[res,R,fit] = evaluateFit(y1,x1,w);
if R>0.9998
return;
end
for i = 1:1:length(x1)
fy(i,1)= sin(ax1(i)+b);
JacobiMatrix(i,1) = calc_pA(x1(i),a,b) ;
...
JacobiMatrix(i,2) = calc_pB(x1(i),a,b);
end
delta_w = inv((JacobiMatrix)'(JacobiMatrix))JacobiMatrix'(y1-fy);
%update a b
a = a +lamdadelta_w(1) ;
b = b+lamdadelta_w(2);
disp(sprintf("Simulation %d fit => sin((%.2f)*x+(%.2f)), R:%.2f",k,a,b,R));
end
%
function [res,R2,fit] = evaluateFit(y,x,w)
fit = getFittingValue(x,w);
res = norm(y-fit)/sqrt(length(fit));
yu = mean(y);
R2 = 1-norm(y-fit)^2/norm(y-yu)^2;
end
...
function fit = getFittingValue(x,w)
len = length(x);
fit = ones(len,1);
for i = 1:1:len
fit(i) = hypothesis(x(i),w);
end
end
function val= hypothesis(x,w)
a = w(1);
b= w(2);
val = sin(ax+b);
end
function val = calc_pA(x,A,B)
val = xcos(Ax+B);
end
function val = calc_pB(x,A,B)
val = cos(Ax+B);
end