Borrar filtros
Borrar filtros

lsqcurvefit not working properly

45 visualizaciones (últimos 30 días)
Jack
Jack el 4 de Jul. de 2024 a las 15:49
Comentada: Jack el 5 de Jul. de 2024 a las 14:08
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')

Respuesta aceptada

Star Strider
Star Strider el 4 de Jul. de 2024 a las 17:52
Editada: Walter Roberson el 4 de Jul. de 2024 a las 18:59
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 comentarios
Jack
Jack el 5 de Jul. de 2024 a las 13:42
Thank you Star!
Star Strider
Star Strider el 5 de Jul. de 2024 a las 13:57
As always, my pleasure!

Iniciar sesión para comentar.

Más respuestas (1)

William Rose
William Rose el 4 de Jul. de 2024 a las 22:03
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))
Size(data)=430 x 247.
%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))
Size(data_new)=401 x 3.
fprintf('Size(E_Fit)=%d x %d.\n',size(E_fit))
Size(E_Fit)=340 x 1.
% 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
Local minimum possible. lsqcurvefit stopped because the relative size of the current step is less than the value of the step size tolerance.
n=1: p0=[-5.79e+00 3481 0.0014]; p=[-7.52e+09 710 0.0007]
Local minimum possible. lsqcurvefit stopped because the relative size of the current step is less than the value of the step size tolerance.
n=2: p0=[-6.58e+00 3481 0.0019]; p=[-1.14e+10 703 0.0009]
Local minimum possible. lsqcurvefit stopped because the relative size of the current step is less than the value of the step size tolerance.
n=3: p0=[-7.34e+00 3481 0.0022]; p=[-3.54e+09 739 0.0013]
% 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.
  2 comentarios
Jack
Jack el 5 de Jul. de 2024 a las 13:47
Thank you William!
Jack
Jack el 5 de Jul. de 2024 a las 14:08
To add on, I think I should have been mroe specific when I say that my fit is not good. Basically, I tried to fit all the 186 delay times, the resulting fitting is far off from the data points. That's what I meant. Hope this is clear.

Iniciar sesión para comentar.

Categorías

Más información sobre Linear and Nonlinear Regression en Help Center y File Exchange.

Etiquetas

Productos


Versión

R2024a

Community Treasure Hunt

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

Start Hunting!

Translated by