フィルターのクリア

lsqcurvefit not working properly

45 ビュー (過去 30 日間)
Jack
Jack 2024 年 7 月 4 日 15:49
コメント済み: Jack 2024 年 7 月 5 日 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')

採用された回答

Star Strider
Star Strider 2024 年 7 月 4 日 17:52
編集済み: Walter Roberson 2024 年 7 月 4 日 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 件のコメント
Jack
Jack 2024 年 7 月 5 日 13:42
Thank you Star!
Star Strider
Star Strider 2024 年 7 月 5 日 13:57
As always, my pleasure!

サインインしてコメントする。

その他の回答 (1 件)

William Rose
William Rose 2024 年 7 月 4 日 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 件のコメント
Jack
Jack 2024 年 7 月 5 日 13:47
Thank you William!
Jack
Jack 2024 年 7 月 5 日 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.

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeLinear and Nonlinear Regression についてさらに検索

タグ

製品


リリース

R2024a

Community Treasure Hunt

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

Start Hunting!

Translated by