Curve fitting with loop

Hi all, I have multiple sets of data to be fitted using a custom fitting function. I would like to do it such that I will only have to provide one initial guess for the first set of data, and after the computer managed to get the solution of the parameters for the first set of data, it will use the solution of the first set of data as the initial guess for the second set of data, so on and forth.

I have attached my data and code for the fitting:

%% Preparation 
clear;clc 

%data = importdata("Experimental data\Transient Absorption\FCPIB-293K-2.5mW-400nm-Jan072021 -ibg -bg -chirp.csv"); % insert file path within parenthesis 
data = importdata("FCPIB-293K-2.5mW-400nm-Jan072021 -ibg -bg -chirp.csv"); % insert file path within parenthesis 


%% Preamble 

% Fundamental constants
h = 4.0135667696*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units
kB = 8.617333268*10^-5; % units: eV/ K

% Clean up of data to select range of values 
wavelength = data(1:end, 1);
delay_t = data(1, 1:end); % conatains all of the delay times
E = (h*c)./(wavelength*10^-9); % contains all of the probe energies
Range_E = E>=1.5 & E<=2.2;
Range_T = delay_t>=0.5 & delay_t<=1000;

% for one delay time
T = find(Range_T);
T_min = min(T);
T_max = max(T);
t = min(T); % choose an integer b/w T_min and T_max
delaytime = delay_t(1, t);


% Data for fitting 
E_p = E(Range_E); % selected probe energies
delta_Abs = -1*data(Range_E,t); 
delta_Abs_norm = delta_Abs./max(delta_Abs); % normalised delta_Abs
Range_Efit = E_p>=1.65 & E_p<=max(E_p);
E_fit = E_p(Range_Efit);
delta_Abs_norm_fit = delta_Abs_norm(Range_Efit);

% Fitting function
function F = MB(y, E_fit)
    F = y(1).*exp(-(E_fit./(8.617333268*10^-5.*y(2)))) + y(3);
end 


%% Curve fitting options 

% Initial parameter guess and bounds 
lb = [0, 293, -Inf]; ub = [Inf, Inf, Inf];
y0 = [4*10^7, 900, 2];


% lsqcurvefit and choose between different algorithm that lsqcurvefit employs (3C1, comment those lines that are not choosen and uncomment the line that is choosen, if not, matlab will take the last line of "optim_lsq" by default)
% optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'levenberg-marquardt', 'MaxFunctionEvaluations',10^5, 'MaxIterations', 10^5, 'FunctionTolerance',10^-10, 'StepTolerance', 10^-10);
optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'trust-region-reflective', 'MaxFunctionEvaluations',10^5, 'MaxIterations',10^5, 'FunctionTolerance',10^-20, 'StepTolerance', 10^-20);
% optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'interior-point', 'MaxFunctionEvaluations',1000, 'MaxIterations', 1000, 'FunctionTolerance',10^-20, 'StepTolerance', 10^-20);

% Solver for lsqcurvefit
[y, residualnorm, residual, exitflag, output, lambda, jacobian] = lsqcurvefit(@MB, y0, E_fit, delta_Abs_norm_fit, lb, ub);

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

%% Error bars calculation 

ci = nlparci(y, residual, 'Jacobian', jacobian);
PCI = table(ci(:,1), y(:), ci(:,2),'VariableNames',{'CI 0.025','y','CI 0.975'});
Parameter_CI = table2array(PCI);
upper_bar = (Parameter_CI(:,3) - Parameter_CI(:,2))./2;
lower_bar = (Parameter_CI(:,2) - Parameter_CI(:,1))./2;

%% Plot command 

plot(E_p, delta_Abs_norm,'Black')
hold on 
plot(E_fit, MB(y, E_fit), 'LineWidth', 1.0, 'Color', 'red')
xlabel('Probe Photon Energy (eV)')
ylabel('Normalised \Delta A (a.u.)')
legend('Experimental Data', 'Fitted Curve')

Some background info regarding the data:

1. There's a range of values in the first column and these corresponds to this quantity called wavelength. After I have extracted the wavelength from there, I'll process it into this quantitiy call E. However, E takes a range of values and I'll only need the values ranging from E = 1.5 to E = 2.2 (this is the data used for the scatterplot, it is named as E_p in my code)
2. After the first step is done, I'll now need to choose the range of E_p (E_p = 1.65 to E_p = 2.2) for the fitting of my data (this is named as E_fit in my code) note: the range of values for fitting ≠ range of values for scatter plot
3. Subsequently, I'll need to look at the first row of my data and select some of them (they need to be positive). Once I have selected them, I'll need to select the data in the same column as them and this needs to correspond to both E_p and E_fit.

In conclusion, the data selected in 1. and 2. goes to the horizontal axis and the data selected in 3. goes to the vertical axis.

Thank you.

Answer 1

[Star Strider]

I would do something like this, puttting the entire code in a loop and saving the fitted parameters (and possibly the plots as well) in each iteration —

csvs = dir('*-chirp.csv');

for k = 1:numel(csvs)
    data = readmatrix(csvs(k).name); % insert file path within parenthesis 
    %% Preamble 
    % Fundamental constants
    h = 4.0135667696*10^-15; % units: eV/ Hz
    c = 3*10^8; % SI units
    kB = 8.617333268*10^-5; % units: eV/ K
    % Clean up of data to select range of values 
    wavelength = data(1:end, 1);
    delay_t = data(1, 1:end); % conatains all of the delay times
    E = (h*c)./(wavelength*10^-9); % contains all of the probe energies
    Range_E = E>=1.5 & E<=2.2;
    Range_T = delay_t>=0.5 & delay_t<=1000;
    % for one delay time
    T = find(Range_T);
    T_min = min(T);
    T_max = max(T);
    t = min(T); % choose an integer b/w T_min and T_max
    delaytime = delay_t(1, t);
    % Data for fitting 
    E_p = E(Range_E); % selected probe energies
    delta_Abs = -1*data(Range_E,t); 
    delta_Abs_norm = delta_Abs./max(delta_Abs); % normalised delta_Abs
    Range_Efit = E_p>=1.65 & E_p<=max(E_p);
    E_fit = E_p(Range_Efit);
    delta_Abs_norm_fit = delta_Abs_norm(Range_Efit);
    % Fitting function
    MB = @(y, E_fit) y(1).*exp(-(E_fit./(8.617333268*10^-5.*y(2)))) + y(3);
    %% Curve fitting options 
    % Initial parameter guess and bounds 
    lb = [0, 293, -Inf]; ub = [Inf, Inf, Inf];
    if k == 1
        y0 = [4*10^7, 900, 2];
    else 
        y0 = yv(k-1,:)
    end
    % lsqcurvefit and choose between different algorithm that lsqcurvefit employs (3C1, comment those lines that are not choosen and uncomment the line that is choosen, if not, matlab will take the last line of "optim_lsq" by default)
    % optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'levenberg-marquardt', 'MaxFunctionEvaluations',10^5, 'MaxIterations', 10^5, 'FunctionTolerance',10^-10, 'StepTolerance', 10^-10);
    optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'trust-region-reflective', 'MaxFunctionEvaluations',10^5, 'MaxIterations',10^5, 'FunctionTolerance',10^-20, 'StepTolerance', 10^-20);
    % optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'interior-point', 'MaxFunctionEvaluations',1000, 'MaxIterations', 1000, 'FunctionTolerance',10^-20, 'StepTolerance', 10^-20);
    % Solver for lsqcurvefit
    [y, residualnorm, residual, exitflag, output, lambda, jacobian] = lsqcurvefit(MB, y0, E_fit, delta_Abs_norm_fit, lb, ub);
    yv(k,:) =  y;
    % % Local minimum possible.
    % % 
    % % lsqcurvefit stopped because the final change in the sum of squares relative to 
    % % its initial value is less than the value of the function tolerance.
    
    %% Error bars calculation 
    ci = nlparci(y, residual, 'Jacobian', jacobian);
    PCI{k} = table(ci(:,1), y(:), ci(:,2),'VariableNames',{'CI 0.025','y','CI 0.975'})
    Parameter_CI = table2array(PCI{k});
    upper_bar = (Parameter_CI(:,3) - Parameter_CI(:,2))./2;
    lower_bar = (Parameter_CI(:,2) - Parameter_CI(:,1))./2;
    
    %% Plot command 
    plot(E_p, delta_Abs_norm,'Black')
    hold on 
    plot(E_fit, MB(y, E_fit), 'LineWidth', 1.0, 'Color', 'red')
    xlabel('Probe Photon Energy (eV)')
    ylabel('Normalised \Delta A (a.u.)')
    legend('Experimental Data', 'Fitted Curve')
end

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

PCI = 1x1 cell array
    {3x3 table}

It would be helpful to have at least one additional .csv file to test this with, so lacking them, I have to leave that to you. Make appropriate changes to the dir call to import only the files you want.

Consider using the saveas function to save the plots (preferably as .png files unless you want to save all the data in them as well, so save them as .fig files in that instance).

.

[Jack]

thank you so much!

[Star Strider]

As always, my pleasure!

Note thtat tthe ‘MB’ anonymous function and several of the constants can be defined outside (before) the loop. I kept them in first because of simplicity, and second, with only three files, the efficiency improvement is likely not significant.

[Image Analyst]

@Jack I don't think this gets around your objections (in your comment to me) of calling lsqcurvefit multiple times and of your data being too huge. But I doubt those were actually valid objections anyway.

Answer 2

[Image Analyst]

So why can't you just do what you said? Just call lsqcurvefit three times, updating variables before the second and third calls?

[Jack]

That's a possible option, but it's inefficient because the data size is too large and I don't know how to code it to make it efficient i.e. not calling lsqcurvefit multiple times

[Image Analyst]

I don't know how you can update the initial parameters and call it again if you don't call it again to see what those parameters are. Even @Star Strider calls it multiple times in his answer. I don't know of any way around it.

MATLAB does have ways of dealing with extremely large data sets, for example memmapfile, see

https://www.mathworks.com/help/matlab/large-files-and-big-data.html

However if it's too large to call lsqcurvefit three times in a row, then it would be too large to call it even once. If it didn't fail the first time then I don't know why it would fail the subsequent times. If you do create temporary variables that are very large, you can call clear to release them from memory.

clear('myHugeVariable');

Worst case, you might be able to subsample your data. Chances are that if you have tens of millions of data points and are trying to fit a curve, you'd get essentially the same curve (same coefficients) if you just used a few thousand data points.