Improving chose of guesses in my fit

Hi everyone,

I do some fit on my data, and my guesses is no good enough.

Here an example of plotting my data:

And this is the code of the function that fit my data:

%fit function
function [Xmax,beta] = fit_function(x,y,peaks,t_peaks,peak_num,check_fit)

% ymax calculated analytically with Wolfram Mathematica
ymax = @(b) (2-b(2)/b(1))*(2*b(1)/b(2)-1)^(-b(2)/2/b(1));
modelfun = @(b,x) b(3)/ymax(b)*exp((x-b(4))*(b(1)-b(2))).*sech(b(1)*(x-b(4)));


bguess=[60, 15, peaks(peak_num), x(t_peaks(peak_num))];
%bguesses is [alpha, beta, amplitude, x offset]
%alpha and beta controlling of the slope of the peak, each one of them control in one slope side.

beta = nlinfit(x, y, modelfun, bguess);
Xmax=(log(-1+(2*beta(1)/beta(2)))/(2*beta(1)))+beta(4);



%check fit:
if check_fit==1
    fittedCurve = modelfun(beta, x);

    hold on;
    plot(x, fittedCurve, 'm');
    legend('Original Data','Fitted Curve');
else
end

end

In this fit, I have 4 guesses. The first two of them(call alpha and beta), chosing as a constant number, according to my attempts to see what would be suitable. The other two of them is depend on the data and according to the data, the code chose the appropriate guesses. (actually is the amplitude and x offset)

So, if you look on the code, you see that 60 and 15 is the constant guesses, and the "peaks(peak_num)" and "x(t_peaks(peak_num))", is the guesses that variable according to the maximum point of the data.

The problem is that in some cases,alpha and beta, the two that are constant, sometimes suitable and sometime unsuitable.

I want to insert guesses, that will be change according to the data, exactlly like the last two guesses of the amplitude and x offset. I thinking about derivative of the peaks or something like this, but I managed to mess with it..

do you have idae for good variable guesses or how to do something that will be appropriate to any try of fitting my data?

thank you all (:

Answer 1

[Matt J]

1. It appears from the plot that you are giving nlinfit the data for all peaks simultaneously. You should not be doing that. Rather, you should be extracting the x,y data local to the peak you are currently fitting. Doing so would likely make the fit more robust to the choice of initial guess
2. The model appears to try fitting both the left and right side of the peak. It seems more prudent to fit the right side of the peak only, which looks like it would follow a much simpler model, e.g., a two-parameter exponential. The left side of the peak looks like a simple jump discontinuity or a linear slope, which needn't be fit explicitly.
3. If you insist on using your current model, an exhaustive search over alpha and beta should be very tractable. It's only two unknown parameters, and your model function is abundantly vectorizable in alpha and beta.

[Image Analyst]

Excellent suggestions, Matt. @Yohay are you sure you really want a formula to fit the whole signal, which may have a variable number of peaks? Maybe you really just want to smooth your data and use the smoothed data in subsequent processing. Explain why you think you need a fit, and what you will do with it assuming you could get it. Also, attach your data in a .mat or text file with the paperclip icon.

[Yohay]

@Image Analyst

Ok, the data comes from a Fabry-Parot scanner of two lasers, where the goal is to make a correction to one of the lasers all the time.

How will it work? I check the difference between two peaks (the first and the third) every few minutes, at the same time I check the difference between two other peaks (the second and the fourth). One difference is my reference, and the other will be corrected according to the firts difference, so that the difference remains the same. So basically what I need is the exact location of each peak, so I thought I'd make a fit.

I Attached one file with the scan data.

[Yohay]

@Matt J

Thanks for the advice!

I'm actually begginer in MATLAB, so I'd appreciate some more detailed guidance. I have no problem changing the model of the fit I did, it's just what I was able to do.

In the comment above I explained why I even need the fit, maybe it will help you understand more what I'm trying to do.

1. Are you actually saying cut my data to only the segment of the specific pik at a time? I'm a bit confused about how to do this in the sense of writing the code...

2. I haven't really been able to find a model function to fit these peaks other than this one, if you have an idea how to write it better with fewer variables that would be great.

finally, i need the locations of the picks so I can do error correction with their data, that's all (:

Thank you!

[Matt J]

You can use findpeaks for peak detection. Then you can extract the i-th inter-peak data by doing,

I=peaks(i):peaks(i+1)-5;

xi=x(I);
yi=y(I);

Answer 2

[Sarthak]

You can try to determine the values of the constant guesses based on the data and use a peak detection algorithm to identify the peaks in the data and then calculate the derivatives at those peak points. The values of alpha and beta can then be set based on the properties of the peaks and their derivatives.

Here's an example of how you can modify your code:

function [Xmax, beta] = fit_function(x, y, peaks, t_peaks, peak_num, check_fit)
    % Calculate the derivatives of the data
    dy = diff(y);
    dx = diff(x);
    derivatives = dy ./ dx;

    % Find the peaks in the data
    [~, peak_indices] = findpeaks(y);

    % Calculate the average derivative at the peak points
    avg_derivative = mean(derivatives(peak_indices));

    % Calculate the values of alpha and beta based on the average derivative
    alpha = avg_derivative / 2;
    beta = 2 * alpha;

    % Use the calculated values of alpha and beta in the model function
    ymax = @(b) (2 - b(2) / b(1)) * (2 * b(1) / b(2) - 1)^(-b(2) / 2 / b(1));
    modelfun = @(b, x) b(3) / ymax(b) * exp((x - b(4)) * (b(1) - b(2))) .* sech(b(1) * (x - b(4)));
    bguess = [alpha, beta, peaks(peak_num), x(t_peaks(peak_num))];

    % Fit the data using the modified guesses
    beta = nlinfit(x, y, modelfun, bguess);
    Xmax = (log(-1 + (2 * beta(1) / beta(2))) / (2 * beta(1))) + beta(4);

    % Plot the fitted curve if check_fit is true
    if check_fit == 1
        fittedCurve = modelfun(beta, x);
        hold on;
        plot(x, fittedCurve, 'm');
        legend('Original Data', 'Fitted Curve');
    end
end

See if this works.