# 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

### More Answers (1)

Image Analyst
on 29 Jun 2024

##### 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

### Community Treasure Hunt

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

Start Hunting!