I am confused about why my code doesn't lead to a fitting result.
2 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Jinglei
el 2 de Feb. de 2023
Editada: Walter Roberson
el 6 de Feb. de 2023
x=[0 1 2 3 4 5 6 7 8 9 10]
xs = 0:0.01:120;
y=[0 4 9 16 25 36 49 64 81 100 121]
myfunc=@fitting;
rng('shuffle');
best_r = inf;
best_p = [-inf, -inf];
best_c = [-inf];
a=0.00001;
b=3;
d=2;
for iters = 1 : 100
p0 = a+(b-a).*rand(1,2);
c0 = a+(d-a).*rand(1,1);
ytotal0 = fitting(p0,c0,xs);
try
[p, c, r] =nlinfit(x,y,myfunc,p0,c0);
if any(p < 0), or (c < 0)
fprintf('rejected negative output on iteration %d\n', iters);
continue;
end
if any(norm(r)==best_r)
continue;
end
nr = norm(r);
if nr < best_r
fprintf('improved on iteration %d\n', iters);
best_r = nr;
best_p = p;
best_c = c;
end
catch ME
fprintf('failed iteration %d\n', iters);
end
end
c = best_c
p = best_p
%c = lsqcurvefit(myfunc,c0,x,y)
ytotal = fitting(p,c,x);
%for adding points
hold on
plot(x,y,'o', 'displayname', 'original')
hold off
xlim auto
ylim auto
legend show
%for adding points
hold on
xs = 0:0.1:150;
ytotal = fitting(p,c,xs);
plot(xs,ytotal)
hold off
R2=1-(sum((fitting(p,c,x)-y).^2)/sum((y-mean(y)).^2));
function ytotal = fitting(p,c,x)
xspan = x;
y0 = zeros(2,1);
[t,ye] = ode15s(@(x,y)eq2(x,y,p,c), xspan, y0);
%protect against failure of integration
if t(end) ~= xspan(end)
ytotal = inf(numel(xspan),1);
else
ytotal = ye(:,1)+ye(:,2)+c(1);
end
end
function dy=eq2(x,y,p,c)
%y(1)=CE,y(2)=LE
dy=zeros(2,1);
dy(1) = p(1) .* x;
dy(2) = p(2);
end
0 comentarios
Respuesta aceptada
Walter Roberson
el 5 de Feb. de 2023
Editada: Walter Roberson
el 5 de Feb. de 2023
I had to fix every different things to get it to work.
The improvement to summarize failed iterations instead of listing each one was cosmetic, not strictly required. It took a while before the iterations started working properly and I was tired of the long stream of fail messages.
x=[0 1 2 3 4 5 6 7 8 9 10]
xs = 0:0.01:120;
y=[0 4 9 16 25 36 49 64 81 100 121]
rng('shuffle');
best_r = inf;
best_p = [-inf, -inf];
best_c = [-inf];
a=0.00001;
b=3;
d=2;
arefailing = false;
firstfail = -inf;
currentfail = -inf;
for iters = 1 : 100
p0 = a+(b-a).*rand(1,2);
c0 = a+(d-a).*rand(1,1);
ytotal0 = fitting(p0,c0,xs);
try
[pc, r] = nlinfit(x,y(:),@(pc,x)fitting(pc(1:2),pc(3),x),[p0,c0]);
%with the try/catch we never reach this statement if the nlfit
%succeeded, so the fact we got here means that an iteration worked
%(but might have returned values we do not want.)
if arefailing
fprintf('failed iterations %d to %d\n', firstfail, currentfail);
arefailing = false;
end
p = pc(1:2); c = pc(3);
if any(p < 0), or (c < 0)
fprintf('rejected negative output on iteration %d\n', iters);
continue;
else
end
if any(norm(r)==best_r)
continue;
end
nr = norm(r);
if nr < best_r
fprintf('improved on iteration %d\n', iters);
best_r = nr;
best_p = p;
best_c = c;
end
catch ME
if arefailing
currentfail = iters;
else
arefailing = true;
firstfail = iters;
currentfail = iters;
disp(ME)
if ~isempty(ME.cause); disp(ME.cause{1}); end
end
end
end
if arefailing
fprintf('failed iterations %d to %d\n', firstfail, currentfail);
if firstfail == 1 %everything failed
fprintf('failed all iterations!!! Not doing any plotting\n');
else
arefailing = false; %we got at least one valid iteration
end
end
if ~arefailing
c = best_c
p = best_p
%c = lsqcurvefit(myfunc,c0,x,y)
ytotal = fitting(p,c,x);
%for adding points
hold on
plot(x,y,'o', 'displayname', 'original')
hold off
xlim auto
ylim auto
legend show
%for adding points
hold on
xs = 0:0.1:150;
ytotal = fitting(p,c,xs);
plot(xs,ytotal)
hold off
R2=1-(sum((fitting(p,c,x)-y).^2)/sum((y-mean(y)).^2));
end
function ytotal = fitting(p,c,x)
xspan = x;
y0 = zeros(2,1);
[t,ye] = ode15s(@(x,y)eq2(x,y,p,c), xspan, y0);
%protect against failure of integration
if t(end) ~= xspan(end)
ytotal = inf(numel(xspan),1);
else
ytotal = ye(:,1)+ye(:,2)+c;
end
end
function dy=eq2(x,y,p,c)
%y(1)=CE,y(2)=LE
dy=zeros(2,1);
dy(1) = p(1) .* x;
dy(2) = p(2);
end
6 comentarios
Walter Roberson
el 6 de Feb. de 2023
In fact, it got stuck in the middle this time, running forever without stopping manually. Do you know what was wrong?
No. I do not know what is wrong with code that I have not seen.
Más respuestas (1)
Bora Eryilmaz
el 2 de Feb. de 2023
If you print out the exception inside the catch block, it will tell you what the error is.
catch ME
ME
fprintf('failed iteration %d\n', iters);
To start with, you are not passing the right number of arguments to the fitting() function: it would be passed two arguments from nlinfit, but your code expects three input arguments.
6 comentarios
Ver también
Categorías
Más información sobre Systems of Nonlinear Equations en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!