I am having trouble using lsqcurvefit with an function

Here is the code I have:

clear all, clc
data=load('OneLorentzian.txt');
t = data(:,1);
y = data(:,2);
plot(t,y,'r.')
F = @(y, ydata) (1/pi*ydata).*( ydata.^2./((t-0).^2 + ydata.^2));
x0 = [1 1 1 0];
[x,resnorm,~,exitflag,output] = lsqcurvefit(F,x0,t,y)
hold on
plot(t,F(x,t))
hold off

When I run the code my function does not fit the spectra. I get this message at the bottom:

Initial point is a local minimum.

Optimization completed because the size of the gradient at the initial point

is less than the value of the optimality tolerance.

<stopping criteria details>

x =

1 1 1 0

resnorm =

1.9478e+07

exitflag =

1

output =

struct with fields:

firstorderopt: 0

iterations: 0

funcCount: 5

cgiterations: 0

algorithm: 'trust-region-reflective'

stepsize: 1

message: '↵Initial

Answer 1

[Star Strider]

Experiment with different values for ‘x0’ since nonlinear parameter estimation routines are very sensitive to the initial parameter estimates. Initial estimates that are far from the global minimum will generally not be changed, as they are not here.

[Russ Rohloff]

I have tried a large range of numbers in all the fields with the same result.

Answer 2

[Walter Roberson]

Initial point is a local minimum.

That is the normal output when the fitting function does the best that can possibly be done.

The function you are using does not use calculus to analyze the function to determine what the global minima of the function is -- and even if it did, there are a lot of functions where the calculus approach is too difficult to navigate, realistically.

Because of this, the fitting function does not assume that it has definitely found the global minima: it says that it found a local minima because it does not have the ability to prove that there is absolutely no possible set of parameters that could do better. It would be overpromising to say otherwise.

Consider the simple function

f = @(x) (x-5).^2 - 1000*(x == 8.2511135319564 )

This has a global minima at 8.2511135319564 but the function cannot analyze the code to determine that it has a special case to analyze, and instead just analyzes the hessians and so on. The hessians converge at the point x = 5 so the function would report x = 5 but would not claim that it is a global minima, just that it is a local minima, since it knows that you might have had something weird in the function that was beyond the ability of the numeric estimation routines to notice.

[Walter Roberson]

F = @(y, ydata) (1/pi*ydata).*( ydata.^2./((t-0).^2 + ydata.^2));

Bad function for fitting purposes. You ignore y and you use the external value t

When the fitting functions run, they pass in a vector of proposed parameters for the function, and a vector of values of the indepenent variable. The function is responsible for calculating the values the set of parameters would predict at the given locations of the independent variable. lsqcurvefit() then subtracts the values of the dependent variable, calculates the residue, figures out how to adjust the parameter.

When the fitting starts out, it may probe at a subset of the full data in order to determine whether the fitting function is behaving properly.

Under the assumption that y is your set of parameters, then your function should probably be

F = @(y, t) (1/pi*t).*( t.^2./((y-0).^2 + t.^2));

but I would not recommend using y as the names of the parameters, as that conflicts with using y as the independent variable.. You are calling the results x and using x0 so you should probably

F = @(x, t) (1/pi*t).*( t.^2./((x-0).^2 + t.^2));

I get the impression that you might have reversed the parameters and the independent variable in your fitting function.

[Walter Roberson]

You are passing in a 1 x 4 vector for your parameters (x0). You are receiving a column vector of values. You are doing vectorized calculations. The result is going to be an N x 4 array that the N x 1 dependent values (y) are subtracted from. This is unlikely to be what you want. Please be explicit about how you want the vector of parameters to affect the output. If that means sum() along the second dimension, then sum() along the second dimension.

I would point out that at the moment, you do not distinguish the meaning of the four different parameters, so there is nothing driving a particular ordering of them.

[Russ Rohloff]

Thanks,

I tried both sets of code but still no luck. here is the curve I am trying to fit. The function is one that was given to me for the spectra I am suppused to fit.

[Russ Rohloff]

This is the first time I have tried to work with fitting a complex curve.

[Walter Roberson]

If you are fitting a lorenzian then you would typically need to fit the location, magnitude, width (or half-width), and phase. If you look at your original function,

F = @(y, ydata) (1/pi*ydata).*( ydata.^2./((t-0).^2 + ydata.^2));

then the t-0 would appear to be an attempt to deal with the location, but it is assuming that the location is 0. If we let be the location of the peak then y would expect a term

Your function is not suitable for fitting a single lorentzian, at least not without making assumptions about some of the parameters... and if you make assumptions about those parameters then you would reduce the number of parameters to be fitted.

[Walter Roberson]

https://www.mathworks.com/matlabcentral/fileexchange/33775-lorentzfit-x-y-varargin

[Russ Rohloff]

Walter,

Thats part of the problem. I am really new at working in matlab and not really sure what I need. I am starting with an arbitary function and the spectra data file. The function was given to me and I am asked to fit the data to the arbitary function and thats where I get stuck. I can get the spectra plotted fine but don't know how to fit the function. I was told to use the lsqcurvefit for the curve. I have had a total of just 6 hours instruction in matlab, so I am not even a rookie. Attached is a screen shot of the function.

[Walter Roberson]

Perhaps your x0 should not be a vector of length 4.

It is not immediately obvious to me what parameter is being fit.

[Torsten]

So for the bell-shaped curve, I suggest using

t = data(:,1);
y = data(:,2);
plot(t,y,'r.')
F = @(x,t) x(1)*exp(-((t-x(2))/x(3)).^2);
x0 = [0.035;1920;5];
[x,resnorm,~,exitflag,output] = lsqcurvefit(F,x0,t,y)
hold on
plot(t,F(x,t))
hold off