Confidence interval calculation. Cannot compute confidence interval for imaginary values.
5 views (last 30 days)
Show older comments
clear; clc
load Steady_State_Data.mat % this contains the wavelength of light and absorbance of substrate and sample
%% Preamble
% Fundamental constants
h = 4.0135667696*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units
% Clean up of data to select range of values
Absorbance = log10(T_substrate./T_sample);
E = (h*c./(Lambda*10^-9));
e = E >= 0 & E <= 2.0; % returns boolean value of the indices that satisfies the logical condition defined above
A = find(e); % gives indices that are non-zero
N = length(A); % no. of elements that are non-zero
n = length(E) + 1;
% Data for fitting
E_p = E(n-N:167);
Abs = Absorbance(n-N:167) - min(Absorbance);
% Abs = Absorbance(n-N:167);
% Fitting function
function F = EM_SS(p, e_p)
for i = 1:numel(e_p)
E_p = e_p(i);
F(i) = p(1)*(2*pi*sqrt(p(4))/E_p)*(1/p(6))*...
(integral(@(E)sech(((E_p - E)./p(6)))*(1 + 10*p(5)*(E - p(3)) + ...
126*(p(5))^2*(E - p(3))^2)/(1 - exp(-2*pi*sqrt(p(4)/(E - p(3))))), p(3), Inf, 'ArrayValued', 1)) + ...
p(2)*(2*pi*p(4)^3/2)*1/p(6)*(...
(1/1^3)*sech((E_p - p(3) + p(4)/1^2)./p(6)) + ...
(1/2^3)*sech((E_p - p(3) + p(4)/2^2)./p(6)) + ...
(1/3^3)*sech((E_p - p(3) + p(4)/3^2)./p(6)) + ...
(1/4^3)*sech((E_p - p(3) + p(4)/4^2)./p(6)) + ...
(1/5^3)*sech((E_p - p(3) + p(4)/5^2)./p(6)) + ...
(1/6^3)*sech((E_p - p(3) + p(4)/6^2)./p(6)) + ...
(1/7^3)*sech((E_p - p(3) + p(4)/7^2)./p(6)));
end
F = F(:);
end
%% Curve fitting options
% Initial parameter guess and bounds
lb = [0,0,0,0,0,0]; ub = [Inf, Inf, 1.65, 0.03, Inf, Inf];
p0 = [0.13, 0.11, 1.4, 0.025, 0.12, 0.035]; % refer to the next line for their order
% p0 = [A1, A2, Eg, Eb, R, g]
% 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',1000, 'MaxIterations', 1000, 'FunctionTolerance',10^-9, 'StepTolerance', 10^-9);
% optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'trust-region-reflective', 'MaxFunctionEvaluations',1000, 'MaxIterations', 'FunctionTolerance',10^-9, 'StepTolerance', 10^-9);
% optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'interior-point', 'MaxFunctionEvaluations',1000, 'MaxIterations', 1000, 'FunctionTolerance',10^-9, 'StepTolerance', 10^-9);
% Solver for lsqcurvefit
[p, residual, exitflag, output, jacobian] = lsqcurvefit(@EM_SS, p0, E_p, Abs, lb, ub);
%% Standard error calculation
ci = nlparci(p, residual, 'Jacobian', jacobian);
lowerbar = p - ci(:,1);
uppperbar = ci(:,2) - p;
if p == real(p)
disp('Real')
else
disp('Complex with non-zero imaginary')
end
if residual == real(residual)
disp('Real')
else
disp('Complex with non-zero imaginary')
end
if J == real(J)
disp('Real')
else
disp('Complex with non-zero imaginary')
end
%% Plot command
plot(E_p, Abs, 'o')
hold on
plot(E_p, EM_SS(p, E_p))
xlabel('Probe energy (eV)')
ylabel('Absorbance.O.D')
legend('Experimental Data', 'Fitted Curve')
hold off
%% Parameter values (refer to command window)
p1 = p(1,1);
p2 = p(1,2);
p3 = p(1,3);
p4 = p(1,4);
p5 = p(1,5);
p6 = p(1,6);
X1 = [' A1 = ', num2str(p1)];
X2 = [' A2 = ', num2str(p2)];
X3 = [' Eg = ', num2str(p3)];
X4 = [' Eb = ', num2str(p4)];
X5 = [' R = ', num2str(p5)];
X6 = [' g = ', num2str(p6)];
disp(X1);
disp(X2);
disp(X3);
disp(X4);
disp(X5);
disp(X6);
Hi all, I am doing a fitting to my experimental data and I am now trying to calculate the confidence interval. However, it seems like there is something wrong with my code and the computer can't really calculate it. May I ask where I am wrong? I tried following the instructions on this link but I still can't get any results.
0 Comments
Accepted Answer
Star Strider
on 12 Jun 2024
Edited: Star Strider
on 12 Jun 2024
The parameters aren’t complex. All the results are real.
The problem with your code is that in MATLAB, function arguments and outputs are positional (except for the name-value pair arguments), and the output of your lsqcurvefit call did not account for all the outputs in order, so you weren’t getting what you thought you were getting.
I provided variable names for all of them here, for the ones you don’t need in your code, you can use the tilde (~) instead. An alternative to my changes would be:
[p, ~, residual, exitflag, output, ~, jacobian] = lsqcurvefit(@EM_SS, p0, E_p, Abs, lb, ub);
With that change, it now works —
clear; clc
load Steady_State_Data.mat % this contains the wavelength of light and absorbance of substrate and sample
% whos('-file','Steady_State_Data')
%% Preamble
% Fundamental constants
h = 4.0135667696*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units
% Clean up of data to select range of values
Absorbance = log10(T_substrate./T_sample);
E = (h*c./(Lambda*10^-9));
e = E >= 0 & E <= 2.0; % returns boolean value of the indices that satisfies the logical condition defined above
A = find(e); % gives indices that are non-zero
N = length(A); % no. of elements that are non-zero
n = length(E) + 1;
% Data for fitting
E_p = E(n-N:167);
Abs = Absorbance(n-N:167) - min(Absorbance);
% Abs = Absorbance(n-N:167);
% Fitting function
function F = EM_SS(p, e_p)
for i = 1:numel(e_p)
E_p = e_p(i);
F(i) = p(1)*(2*pi*sqrt(p(4))/E_p)*(1/p(6))*...
(integral(@(E)sech(((E_p - E)./p(6)))*(1 + 10*p(5)*(E - p(3)) + ...
126*(p(5))^2*(E - p(3))^2)/(1 - exp(-2*pi*sqrt(p(4)/(E - p(3))))), p(3), Inf, 'ArrayValued', 1)) + ...
p(2)*(2*pi*p(4)^3/2)*1/p(6)*(...
(1/1^3)*sech((E_p - p(3) + p(4)/1^2)./p(6)) + ...
(1/2^3)*sech((E_p - p(3) + p(4)/2^2)./p(6)) + ...
(1/3^3)*sech((E_p - p(3) + p(4)/3^2)./p(6)) + ...
(1/4^3)*sech((E_p - p(3) + p(4)/4^2)./p(6)) + ...
(1/5^3)*sech((E_p - p(3) + p(4)/5^2)./p(6)) + ...
(1/6^3)*sech((E_p - p(3) + p(4)/6^2)./p(6)) + ...
(1/7^3)*sech((E_p - p(3) + p(4)/7^2)./p(6)));
end
F = F(:);
end
%% Curve fitting options
% Initial parameter guess and bounds
lb = [0,0,0,0,0,0]; ub = [Inf, Inf, 1.65, 0.03, Inf, Inf];
p0 = [0.13, 0.11, 1.4, 0.025, 0.12, 0.035]; % refer to the next line for their order
% p0 = [A1, A2, Eg, Eb, R, g]
% 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',1000, 'MaxIterations', 1000, 'FunctionTolerance',10^-9, 'StepTolerance', 10^-9);
% optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'trust-region-reflective', 'MaxFunctionEvaluations',1000, 'MaxIterations', 'FunctionTolerance',10^-9, 'StepTolerance', 10^-9);
% optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'interior-point', 'MaxFunctionEvaluations',1000, 'MaxIterations', 1000, 'FunctionTolerance',10^-9, 'StepTolerance', 10^-9);
% Solver for lsqcurvefit
[p, residualnorm, residual, exitflag, output, lambda, jacobian] = lsqcurvefit(@EM_SS, p0, E_p, Abs, lb, ub);
% p
% residual
% jacobian
%% Standard error calculation
ci = nlparci(p, residual, 'Jacobian', jacobian);
% lowerbar = p - ci(:,1);
% uppperbar = ci(:,2) - p;
Parameter_Confidence_Intervals = table(ci(:,1), p(:), ci(:,2),'VariableNames',{'CI 0.025','p','CI 0.975'})
[Ypred,delta] = nlpredci(@EM_SS,E_p,p,residual,'Jacobian',jacobian); % Prediction Confidence Limits
% if p == real(p)
% disp('Real')
% else
% disp('Complex with non-zero imaginary')
% end
%
% if residual == real(residual)
% disp('Real')
% else
% disp('Complex with non-zero imaginary')
% end
%
% if J == real(J)
% disp('Real')
% else
% disp('Complex with non-zero imaginary')
% end
%% Plot command
plot(E_p, Abs, '.')
hold on
plot(E_p, Ypred, '-r')
errorbar(E_p, Ypred, delta, 'r')
xlabel('Probe energy (eV)')
ylabel('Absorbance.O.D')
legend('Experimental Data', 'Fitted Curve', 'Prediction Confidence Intervals', 'Location','best')
hold off
%% Parameter values (refer to command window)
p1 = p(1,1);
p2 = p(1,2);
p3 = p(1,3);
p4 = p(1,4);
p5 = p(1,5);
p6 = p(1,6);
X1 = [' A1 = ', num2str(p1)];
X2 = [' A2 = ', num2str(p2)];
X3 = [' Eg = ', num2str(p3)];
X4 = [' Eb = ', num2str(p4)];
X5 = [' R = ', num2str(p5)];
X6 = [' g = ', num2str(p6)];
disp(X1);
disp(X2);
disp(X3);
disp(X4);
disp(X5);
disp(X6);
The second parameter is not significantly different from zero (the confidence limits are of opposite signs, so include zero), so it may not be necessary in your model. An alternative approach, if you have the Global Optimization Toolbox, would be to use the ga (genetic algorithm) function (my favourite, or one of the other multistart functions) first to search for the global minimum, and then use lsqcurvefit with those initial parameter estimates to refine the parameter estimates.
EDIT — (12 Jun 2024 at 05:47)
I just remembered that you want error bars, and the parameter confidence limits, while important, will not produce them. I added the nlpredci call to calculate the confidence intervals on the predictions (here, the fitted curve), and then used the errorbar function to plot them. (The error bars are vanishingly small in the plot, so you may need to zoom in to see them distinctly.)
.
2 Comments
Star Strider
on 12 Jun 2024
As always, my pleasure!
No worries! I could have labeled them a bit differently, however I generally strive for accuracy.
More Answers (0)
See Also
Categories
Find more on Annotations 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!