confidence interval of non linear fit of multiple data set

Hi everyone,

I am trying to estimate the confidence interval of my fitting parameters. I have written the code with the constraints and it works properly but now I want to estimate the confidence intervals for the fitting parameters.

I decided to use bootci but I am not able to include this function in my code.

Do you have any suggestion?

The code is attached to the question:

clear
%fit per tln+glu under Tg
R=8.314; 
T1=220; %temperature reffered to y1
T2=300; %temperature reffered to y2

yfcn1 = @(b,x) b(1).*exp(-x(:,2).^2.*b(2)).*(1-2*(exp(-b(3)./(R.*x(:,1))+b(4)/R)./(1+exp(-b(3)./(R.*x(:,1))+b(4)/R)).^2).*(1-sin(x(:,2).*b(5))./(x(:,2)*b(5))));
yfcn2 = @(b,x) b(6).*exp(-x(:,2).^2.*b(7)).*(1-2*(exp(-b(3)./(R.*x(:,1))+b(4)/R)./(1+exp(-b(3)./(R.*x(:,1))+b(4)/R)).^2).*(1-sin(x(:,2).*b(5))./(x(:,2)*b(5))));
x=[0.5215    0.7756    1.2679    1.4701    1.6702    1.8680    2.0633    2.2693    2.4584    2.6442    2.8264    3.0046  3.0890    3.2611    3.4287    3.5917    3.7497    3.9309    4.0774    4.2183    4.3535    4.4827    4.5427    4.6628]; 
y1=[0.9936    0.9375    0.9081    0.8648    0.8568    0.8114    0.7711    0.8010    0.7884    0.7389    0.7901    0.7825  0.7903    0.7501    0.7070    0.7489    0.6441    0.7105    0.6735    0.6385    0.6357    0.6962    0.5946    0.6783]; 
y1_err= [ 0.0637    0.0526    0.0330    0.0235    0.0298    0.0223    0.0388    0.0223    0.0333    0.0326    0.0410    0.0282 0.0561    0.0235    0.0271    0.0218    0.0333    0.0252    0.0344    0.0261    0.0499    0.0396    0.0655    0.0901];
y2=[0.8748    0.8726    0.7922    0.7782  0.7396    0.6958    0.6603    0.6503    0.6556    0.6461    0.6021    0.5820 0.6220    0.5768    0.4950    0.5300    0.5234    0.5170    0.4369    0.4508    0.4409    0.4392    0.4100    0.6699];
y2_err=[ 0.0562    0.0480    0.0287    0.0211    0.0260    0.0194    0.0339    0.0188    0.0287    0.0289    0.0332    0.0225   0.0460    0.0191    0.0211    0.0169    0.0280    0.0198    0.0256    0.0204    0.0392    0.0283    0.0504    0.0856]; 


format long E
T1=220; %temperature reffered to y1
T2=300; %temperature reffered to y2
T1v = T1*ones(size(x));
T2v = T2*ones(size(x));
%T3v = T3*ones(size(x));
xm = x(:)*ones(1,2);
ym = [y1(:) y2(:)];%, y3(:)];
Tm = [T1v(:) T2v(:)];% T3v(:) ];
yerr=[y1_err(:) y2_err(:)];% y3_err(:)];
xv = xm(:);
yv = ym(:);
Tv = Tm(:);
yerrv=yerr(:);
weights=1./yerrv;
xTm = [Tv xv];

%B0 = randn(7,1)*0.1;
B0=[0.5 1e-3 0.5 1e-3 5e4 45 1.5]';

yfcn = @(b,x) yfcn1(b,x).*(x(:,1)==T1) + yfcn2(b,x).*(x(:,1)==T2);
fitfcnw = @(b) norm(weights.*(yv - yfcn(b,xTm)));
lb=[0,0,1e3,10,0,0,0];
ub=[1,1,4e5,1000,2,1,1];
A = @simple_constraint;
problem = createOptimProblem('fmincon', 'x0',B0,'nonlcon',A, 'objective',fitfcnw,'lb',lb,'ub',ub);%
gs = GlobalSearch('PlotFcns',@gsplotbestf);
[B,fval] = run(gs,problem)
B(:)
format short eng

figure(1)
for k = 1:2
    idx = (1:numel(x))+numel(x)*(k-1);
    subplot(2,1,k)
    errorbar(x.^2, ym(:,k),yerr(:,k), '.')
    hold on
    plot(x.^2, yfcn(B,xTm(idx,:)), '-r')
    hold off
    grid
    ylabel('Substance [Units]')
    title(sprintf('y_{%d}, T = %d', k,xTm(idx(1),1)))
    ylim([min(yv) max(yv)+0.2])
    if k == 1
        text(5, max(yv)+0.1, sprintf('$y = %.3f\\times e^{-x^2\\times %.3f}(1-2\\frac{e^{\\frac{%.3f}{%3d\\times R}-\\frac{%.3f}{R}}}{(1+e^{\\frac{%.3f}{%3d\\times R}-\\frac{%.3f}{R}})^2}(1-\\frac{sin(%.3fx}{%.3fx}))$',B(1:3),T1,B(4),B(3),T1,B(4:5),B(5)), 'Interpreter','latex', 'FontSize',12)
    elseif k == 2
        text(5, max(yv)+0.1, sprintf('$y = %.3f\\times e^{-x^2\\times %.3f}(1-2\\frac{e^{\\frac{%.3f}{%3d\\times R}-\\frac{%.3f}{R}}}{(1+e^{\\frac{%.3f}{%3d\\times R}-\\frac{%.3f}{R}})^2}(1-\\frac{sin(%.3fx}{%.3fx}))$',B(6:7),B(3),T2,B(4),B(3),T2,B(4:5),B(5)), 'Interpreter','latex', 'FontSize',12)
    end
end
xlabel('Q^2')

Answer 1

[Aditya]

To estimate the confidence intervals for your fitting parameters using the bootci function, you need to follow a few steps. Since bootci is typically used for bootstrap confidence intervals and requires resampling your data, it's not directly applicable to the output of optimization routines like fmincon used in your case. However, you can simulate a similar bootstrap process by resampling your residuals or performing the fitting process multiple times on subsets of your data.

Here's a simplified implementation outline:

% Assuming B contains your best-fit parameters
nBoot = 1000; % Number of bootstrap samples
bootParams = zeros(nBoot, length(B)); % Store bootstrapped parameters

% Residuals from the original fit
originalFit = yfcn(B, xTm);
residuals = yv - originalFit;

for i = 1:nBoot
    % Resample residuals and add to the original fit
    newResiduals = residuals(randsample(length(residuals), length(residuals), true));
    newY = originalFit + newResiduals;
    
    % Refit model to the new dataset (simplified, see note below)
    try
        newFit = fmincon(@(b) norm(weights.*(newY - yfcn(b,xTm))), B0, [], [], [], [], lb, ub, A);
        bootParams(i, :) = newFit;
    catch
        warning('Failed to converge for bootstrap sample %d', i);
    end
end

% Calculate confidence intervals for each parameter
confIntervals = prctile(bootParams, [2.5, 97.5]);

% Display confidence intervals
disp(confIntervals);

his code snippet omits some details for brevity and clarity. Specifically, you need to adjust the refitting process (fmincon call inside the loop) to match your full model fitting procedure, including handling any constraints or options you used in the original fit. The try-catch block is used to handle cases where fmincon might fail to converge for some bootstrap samples.