Why doesn't my fit converge?

15 ビュー (過去 30 日間)
Corey Crandall
Corey Crandall 2019 年 8 月 24 日
コメント済み: Corey Crandall 2019 年 8 月 24 日
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 件のコメント
John D'Errico
John D'Errico 2019 年 8 月 24 日
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 2019 年 8 月 24 日
Not really I suppose, I'm just trying to learn through practice.

サインインしてコメントする。

採用された回答

John D'Errico
John D'Errico 2019 年 8 月 24 日
編集済み: John D'Errico 2019 年 8 月 24 日
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 件のコメント
Corey Crandall
Corey Crandall 2019 年 8 月 24 日
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.

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeEigenvalue Problems についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by