lsqcurvefit not working properly
7 views (last 30 days)
Show older comments
Hi all, I have multiple sets of data that are to be fitted by my fitting function. Since these data sets are kind of related, what I did was to use an inital guess to find the parameter values for the first dataset, and then use the parameter values as the guess for the second dataset so on and forth. However, the fitting was really bad. May I ask if there is another way around this kind of problem involving fitting of multiple datasets?
tic
%% Preparation
clear; clc
data = importdata("Experimental data\Transient Absorption\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
% Data
wavelength = data(1:end, 1); % contains all of the wavelengths
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 of data needed
Range_E = E>=1.5 & E<=2.2;
Range_W = wavelength >= (h*c)/(2.2*10^-9) & wavelength <= (h*c)/(1.5*10^-9);
Range_T = delay_t>=0.5 & delay_t<=1000;
E_p = E(Range_E); % probe energies for scatterplot
w_p = wavelength(Range_W);
delay_T = delay_t(Range_T);
Range_Efit = E_p>=1.62 & E_p<=max(E_p);
Range_Wfit = w_p >= (h*c)/(2.2*10^-9) & w_p <= (h*c)/(1.62*10^-9);
E_fit = E_p(Range_Efit); % probe energies for fitting
W_fit = w_p(Range_Wfit);
t_min = min(delay_T);
t_max = max(delay_T);
w_min = min(w_p);
w_max = max(w_p);
[row1, col1] = find(data == w_min);
[row2, col2] = find(data == w_max);
[row3, col3] = find(data == t_min);
[row4, col4] = find(data == t_max);
% New cleaned up data
data_new = data(row1:row2, col3:col4);% new data containing required delta A
% for n = 1:length(delay_T)
% delta_Abs(:,n) = -1*data_new(:,n);
% delta_Abs_norm(:,n) = delta_Abs(:,n)/max(delta_Abs(:,n));
% delta_Abs_fit(:,n) = data_new(Range_Wfit,n);
% delta_Abs_norm_fit(:,n) = delta_Abs_fit(:,n)/max(delta_Abs_fit(:, n));
% % plot command
% plot(E_p, delta_Abs_norm)
% xlabel('Probe Energy (eV)')
% ylabel('Normalised \Delta A (a.u.)')
% legend('Experimental Data')
% end
% Fitting function: Maxwell-Boltzmann distribution
function F = MB(y, e_fit)
data = importdata("Experimental data\Transient Absorption\FCPIB-293K-2.5mW-400nm-Jan072021 -ibg -bg -chirp.csv");
h = 4.0135667696*10^-15;
c = 3*10^8;
kB = 8.617333268*10^-5;
delay_t = data(1, :);
Range_T = delay_t>=0.5 & delay_t<=1000;
delay_T = delay_t(Range_T);
wavelength = data(1:end, 1);
E = (h*c)./(wavelength*10^-9);
Range_E = E>=1.5 & E<=2.2;
E_p = E(Range_E);
Range_Efit = E_p>=1.62 & E_p<=max(E_p);
E_fit = E_p(Range_Efit);
for i = 1:length(E_fit)
E_fit = e_fit(i);
F(i) = y(1).*exp(-(E_fit./(kB.*y(2)))) + y(3);
end
F = F(:);
end
% For loop to create required datasets for plotting
for n = 2:length(delay_T)
delta_Abs(:,1) = -1*data_new(:,1);
delta_Abs_norm(:,1) = delta_Abs(:,1)./max(abs(delta_Abs(:,1)));
delta_Abs_fit(:,1) = data_new(Range_Wfit,1);
delta_Abs_norm_fit(:,1) = delta_Abs_norm(Range_Wfit);
delta_Abs(:,n) = -1*data_new(:,n);
delta_Abs_norm(:,n) = delta_Abs(:,n)./max(abs(delta_Abs(:,n)));
delta_Abs_fit(:,n) = data_new(Range_Wfit,n);
delta_Abs_norm_fit(:,n) = delta_Abs_norm(Range_Wfit);
% lsqcurvefit
optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'levenberg-marquardt', 'MaxFunctionEvaluations',10^20, 'MaxIterations', 10^20, 'FunctionTolerance',10^-20, 'StepTolerance', 10^-20);
lb = [0, 293, -1]; ub = [Inf, 6000, 1];
lb1 = [0, 293, -1]; ub1 = [Inf, 2000, 1];
lb2 = [0, 293, -1]; ub2 = [Inf, 1200, 1];
lb3 = [0, 293, -1]; ub3 = [Inf, 1000, 1];
lb4 = [0, 293, -1]; ub4 = [Inf, 750, 1];
y1 = [2*10^7, 1000, 0.5];
y2 = [2*10^7, 800, 0.5];
y3 = [2*10^10, 800, 0.5];
y4 = [2*10^11, 600, 0.5];
y5 = [2*10^14, 600, 0.5];
y(1,:) = lsqcurvefit(@MB, y1, E_fit, delta_Abs_norm_fit(:, 1), lb1, ub1, optim_lsq);
y(n,:) = lsqcurvefit(@MB, y(n-1, :), E_fit, delta_Abs_norm_fit(:, n), lb, ub, optim_lsq);
% plot command
plot(E_p, delta_Abs_norm)
hold on
plot(E_fit, MB(y(1,:), E_fit), 'red')
plot(E_fit, MB(y(n,:), E_fit), 'red')
xlabel('Probe Energy (eV)')
ylabel('Normalised \Delta A (a.u.)')
legend('Fitted curve')
end
carrier_T = y(:,2);
disp(carrier_T)
toc
%% Miscellaneous
% For a rectangular pic:
% drag horizontally until it just covers the letter 'c' in the word 'col'
% For a square pic:
% drag horizontall until it just covers the letter 'o' in the word 'Contribution' in the next line)
% legend('Experimental Data', 'Fitted Curve', 'Carrier Contribution','Excitonic Contribution')
0 Comments
Accepted Answer
Star Strider
on 4 Jul 2024
Edited: Walter Roberson
on 4 Jul 2024
It is difficult for me to understand what you are doing in this version of your code.
You have a curve that is in part an exponential decay, so one option to estimate tthe initial parameters would be to linearise it, do a linear regression on it, re-transform the estimated parameters to use with your nonlinear regression function.
Lv = any(delta_Abs_norm_fit > 0,2);
B = [log(E_fit(Lv)) ones(size(E_fit(Lv)))] \ log(delta_Abs_norm_fit(Lv,:));
figure
plot(E_fit, delta_Abs_norm_fit)
hold on
plot(E_fit(Lv), exp([log(E_fit(Lv)) ones(size(E_fit(Lv)))]*B(:,1)))
hold off
title('Linear Regression Parameter Estimate Results')
The transformation would then be:
y1 = [exp(B(2,1)) B(1,1) rand]
The purpose of the ‘Lv’ calculation is to avoid taking the logarithms of negative numbers.
When I experimented with that approach with your other version of your code, it worked. The upside is that it’s reasonably efficient.
Experiment with it to see if it works in this version of your code.
.
8 Comments
More Answers (1)
William Rose
on 4 Jul 2024
It would be helpful to run your script in the window to show that the "really bad" fit looks like. Post only the minimum code and data necessary to illustrate the issue. Modify your script so it runs in the Matlab answers window.
The value of Planck constant in eV/Hz is incorrect by a few percent; fixed below.
I had an error when I tried to run. The error went away when I changed the filename to something more simple.
Then the code timed out ("Your code took longer than 235 seconds to run. Simplify code and then run again.").
Your data file has >300 columns and 1024 rows of data. Plus a header line (the delay times) and lines of metadata at the end. Much of the data is NaNs.
You select data at a range of wavelengths (i.e. a range of rows) for plotting, from 547.8 to 802.5 nm. You select a subset of this data for fitting. The dashed line in plot below shows the begigning of the range for fitting. I edited the data file to eliminate rows and columns outside this range of plotting interest. And I adjusted Range_T so that only 3 spectra (3 columns) would be fitted.
The fitting function has 3 parameters: p(1:3). The fitting function is a decaying exponential function of energy, with parameters p(1)=amplitude of exponential, p(2)=temperature [K], which determines the energy constant for decay of the exponential, and p(3)=vertical offset of the baseline.
The absorbance data, in the energy range of interest, looks like an upside-down decaying exponential. The absorbance data wavelength increases with row number, so the energy decreases with row number. So the highest energies are the lowest number rows, and the lowest energy is the last row.
For the three columns plotted below, it looks like a good initial guess for the parameters is p0(2)=0.3 eV/kB~=3500 K, p0(3)=+0.002, and p0(1)=(ydata(end)-p0(3))/exp(-E(end)/(kB*p0(2))). A better initial guess for p0(3) is the mean value of the absorbance at the highest energy range being fitted.
data=importdata('spectrum1a.csv');
fprintf('Size(data)=%d x %.d.\n',size(data))
%data = importdata("FCPIB-293K-2.5mW-400nm-Jan072021 -ibg -bg -chirp.csv");
%% Preamble
% Fundamental constants
h = 4.135667696*10^-15; % units: eV/ Hz
%h = 4.0135667696*10^-15;
c = 3*10^8; % SI units
kB = 8.617333268*10^-5; % units: eV/ K
% Data
wavelength = data(1:end, 1); % wavelengths [nm]
delay_t = data(1, 1:end); % delay times
E = (h*c)./(wavelength*10^-9); % probe energies [eV]
% Range of data needed
Range_E = E>=1.5 & E<=2.2; % [eV]
%Range_T = delay_t>=0.5 & delay_t<=1000;
Range_T = delay_t>=0.5 & delay_t<=0.7;
E_p = E(Range_E); % probe energies for scatterplot
w_p = wavelength(Range_E);
% w_p = wavelength(Range_W);
delay_T = delay_t(Range_T);
%Range_Efit = E_p>=1.62 & E_p<=max(E_p);
Emin=1.62; % minimum energy for fitting [eV]
Range_Efit = E_p>=Emin;
E_fit = E_p(Range_Efit); % probe energies for fitting
t_min = min(delay_T);
t_max = max(delay_T);
w_min = min(w_p);
w_max = max(w_p);
[row1, ~] = min(find(data(:,1) == w_min));
[row2, ~] = max(find(data(:,1) == w_max));
[~, col1] = find(data == t_min);
[~, col2] = find(data == t_max);
% New cleaned up data
data_new = data(row1:row2, col1:col2);% new data containing required delta A
fprintf('Size(data_new)=%d x %.d.\n',size(data_new))
fprintf('Size(E_Fit)=%d x %d.\n',size(E_fit))
% Plot columns 1-3
figure
plot(E_p,data_new(:,1),'-r',E_p,data_new(:,2),'-g',E_p,data_new(:,3),'-b')
hold on; xline(Emin,'--k');
grid on; xlabel('Energy (eV)'); ylabel('\Delta Absorbance')
legend(['delay=',num2str(delay_T(1))],['delay=',num2str(delay_T(2))],['delay=',num2str(delay_T(3))],'Location','southeast')
% Fitting function: Maxwell-Boltzmann distribution
function F = MB(p, e)
kB = 8.617333268*10^-5; % [eV/K]
F = p(1).*exp(-(e./(kB.*p(2)))) + p(3);
end
% For loop to fit the data
figure
p=zeros(length(delay_T),3); % allocate array for fitted parameters
for n = 1:length(delay_T)
ydata=data_new(Range_Efit,n);
p0(3)=mean(ydata(1:10)); % initial guess
p0(2)=0.3/kB; % initial guess [deg K]
p0(1)=(ydata(end)-p0(3))/exp(-E_fit(end)/(kB*p0(2))); % initial guess
optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'levenberg-marquardt', 'MaxFunctionEvaluations',10^20, 'MaxIterations', 10^20, 'FunctionTolerance',10^-20, 'StepTolerance', 10^-20);
lb = [-Inf, 293, -1]; ub = [Inf, 6000, 1];
p(n,:) = lsqcurvefit(@MB, p0, E_fit, ydata, lb, ub, optim_lsq);
fprintf('n=%d: p0=[%.2e %.0f %.4f]; p=[%.2e %.0f %.4f]\n',n,p0,p(n,:))
end
% Plot data and fitted curves
figure
plot(E_fit,data_new(Range_Efit,1),'-r',E_fit,data_new(Range_Efit,2),'-g',E_fit,data_new(Range_Efit,3),'-b')
hold on
plot(E_fit,MB(p(1,:),E_fit),'--r',E_fit,MB(p(2,:),E_fit),'--g',E_fit,MB(p(3,:),E_fit),'--b')
grid on; xlabel('Energy (eV)'); ylabel('\Delta Absorbance')
legend(['delay=',num2str(delay_T(1))],['delay=',num2str(delay_T(2))],['delay=',num2str(delay_T(3))],...
'fit 1','fit 2','fit 3','Location','southeast')
The three fits look reasonable. The fitted p(1) and p(2) are not close to the initial gueeses, but it seems to work, i.e. the fitting routine finds a reasonable fit.
See Also
Categories
Find more on 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!