Why doesn't my fit converge?

15 Ansichten (letzte 30 Tage)
Corey Crandall
Corey Crandall am 24 Aug. 2019
Kommentiert: Corey Crandall am 24 Aug. 2019
Hello,
I have some data that I am simulating and would like to fit a function to. The function I am trying to fit to the data is given below, where I am trying to do a monte carlo type simulation. The parameters are varied by a small amount, and if the sum of the residuals squared is less than before, the new guesses are kept. If it is greater, the old parameters are kept and then are varied by a different amount. With a large number of runs, one would expect the parameters to converge to the correct ones. But for some reason, the parameters just don't seem to converge (even for as many as 2000000 runs). I've been banging my head against the wall for a while and was wondering if anyone can see something I'm missing.
clearvars;
clc;
runs = 200000;
x = -15:0.1:15;
%The domain of the data.
y = 3*(sinc(0.2*x + 0.5)).^2 +0.2;
%The original function--the parameters are arbitrary, it's just something I
%want the fit to converge to.
a=zeros(size(runs+1)); b=zeros(size(runs+1)); c=zeros(size(runs+1)); d=zeros(size(runs+1)); fun = zeros(runs+1,301);
%Initialize variables.
a(1) = 3.4; b(1) = 0.25;c(1) = 0.54; d(1) = 0.15;
fun(1,:) = a(1)*(sinc(b(1).*x + c(1))).^2 +d(1);
%Initial guesses
residuals = (fun(1,:) - y).^2;
chi_squared(1) = sum(residuals);
%Calculate chi squared for funtion with initial guesses for parameters.
tic;
for i = 1:runs
a(i+1) = a(i) + 0.001*rand(1); a_temp = a(i+1);
b(i+1) = b(i) + 0.001*rand(1); b_temp = b(i+1);
c(i+1) = c(i) + 0.001*rand(1); c_temp = c(i+1);
d(i+1) = d(i) + 0.001*rand(1); d_temp = d(i+1);
%Adds a random number to vary the parameters slightly.
fun(i+1,:) = a(1)*(sinc(b(1).*x + c(i+1))).^2 +d(1);
residuals = (fun(i+1,:) - y).^2;
chi_squared(i+1) = sum(residuals);
%Calculates the chi squared for the new guesses.
[~,I] = min(chi_squared);
%Finds the where the minimum value is for the chi-squared. If the
%new guesses for the parameters gives a minimum for the chi
%squared (meaning: it is less than the previous run) then the new
%guesses for the parameters are kept; if the chi squuared is larger
%for this run than the previous run, the parameters are not updated.
if chi_squared(i+1) > min(chi_squared)
a(i+1) = a(i);
b(i+1) = b(i);
c(i+1) = c(i);
d(i+1) = d(i);
else
a(i+1) = a_temp;
b(i+1) = b_temp;
c(i+1) = c_temp;
d(i+1) = d_temp;
end
end
toc;
figure
hold on
y_fit = a(runs+1)*(sinc(b(runs+1).*x + c(1))).^2 +d(runs+1);
plot(x,y,'.')
plot(x,y_fit)
figure
plot(chi_squared)
  2 Kommentare
John D'Errico
John D'Errico am 24 Aug. 2019
Is there a good reason why you are writing code to do (apparently poorly) what high quality professionally written codes are available to do well?
Corey Crandall
Corey Crandall am 24 Aug. 2019
Not really I suppose, I'm just trying to learn through practice.

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

John D'Errico
John D'Errico am 24 Aug. 2019
Bearbeitet: John D'Errico am 24 Aug. 2019
Assuming you have some good reason to want to do this (something I seriously doubt is a good reason)...
Why is it not converging? READ THE HELP FOR THE TOOLS YOU ARE USING.
What are the properties of rand? (Hint: uniform random numbers, between 0 and 1.)
What is the actual parameter set here?
y = 3*(sinc(0.2*x + 0.5)).^2 +0.2;
Answer: [3, 0.2, 0.5, 0.2]
What is your starting guess?
a(1) = 3.4; b(1) = 0.25;c(1) = 0.54; d(1) = 0.15;
Answer: [3.4, 0.25, 0.54, 0.15]
What direction do most of those parameters need to move?
Answer: Lower, for 3 of the 4 parameters.
Can this ever happen if your random steps are always greater than zero?
Answer: No.
What random steps are you taking?
a(i+1) = a(i) + 0.001*rand(1); a_temp = a(i+1);
b(i+1) = b(i) + 0.001*rand(1); b_temp = b(i+1);
c(i+1) = c(i) + 0.001*rand(1); c_temp = c(i+1);
d(i+1) = d(i) + 0.001*rand(1); d_temp = d(i+1);
Answer: uniformly distributed numbers in the interval (0,0.001), so ALWAYS positive
The solution is to not write code to do what pre-existing, well written code is already available to solve. That is a recipe for failure.
I'm sorry, but the code you have written is a pretty poor optimization algorithm. Learn to use existing tools, perhaps the optimization toolbox, the curve fitting toolbox, the stats toolbox, the global optimization toolbox. If you have none of them, then learn to use tools like fminsearch, which is entirely capable of solving your problem. Also you should learn about things like preallocation, so as to not dynamically grow vectors, but that is just good programming technique.
  1 Kommentar
Corey Crandall
Corey Crandall am 24 Aug. 2019
Thanks! I don't know how I didn't see that. As I said it's just for practice, and I'm still pretty new to optimization in general. I'll look into those things you suggested.

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by