Curve fitting with loop
4 views (last 30 days)
Show older comments
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);
%% 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:
- 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)
- 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
- 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.
0 Comments
Accepted Answer
Star Strider
on 29 Jun 2024
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
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).
.
3 Comments
Star Strider
on 1 Jul 2024
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
on 2 Jul 2024
@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.
More Answers (1)
Image Analyst
on 29 Jun 2024
So why can't you just do what you said? Just call lsqcurvefit three times, updating variables before the second and third calls?
2 Comments
Image Analyst
on 30 Jun 2024
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.
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.
See Also
Categories
Find more on Get Started with Curve Fitting Toolbox in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!