How to calculate Bootstrap confidence interval

Hi all. I wonder how to calculate bootstrapped confidence interval for following data sent. The conventiona confidence interval is quite large and I wonder if it may make more sense with bootstrapped CIs.

x = [6 10 14 20 26 34 38];

y = [122 107 119 120 105 95 92];

I am using a*exp(-x/b) for fitting which gives a = 127, b = 130. I want to calculated the CIs for bootstrapped sample.

Answer 1

[the cyclist]

As you may have noticed already, the data have a huge variance around that fit:

x = [6 10 14 20 26 34 38]';
y = [122 107 119 120 105 95 92]';

% Tabulate the data
tbl = table(x,y);

% %  Fit the model

% Define function that will be used to fit data
% (F is a vector of fitting parameters)
f = @(F,x) F(1).*exp(-x/F(2));

beta0 = [1 1];
mdl = fitnlm(tbl,f,beta0)

mdl = 
Nonlinear regression model:
    y ~ F1*exp( - x/F2)

Estimated Coefficients:
          Estimate      SE      tStat       pValue  
          ________    ______    ______    __________

    F1     127.28     6.6636    19.101    7.2494e-06
    F2     130.02     39.984    3.2518      0.022651


Number of observations: 7, Error degrees of freedom: 5
Root Mean Squared Error: 7.44
R-Squared: 0.69,  Adjusted R-Squared 0.628
F-statistic vs. zero model: 750, p-value = 6.36e-07

    
% Calculate the model values at the empirical x
y_predicted = predict(mdl,x);

% Plot the data and fit
figure
plot(x,y,'*',x,y_predicted,'g');
legend('data','fit')

But, here is the code to build the bootstrap intervals. b in particular has a huge variance, and the fit to the resampled data off throws off warnings. I turned off warnings here, so you don't have to scroll through all of them to see the output, but I recommend leaving them on locally, so that you can investigate them. Caveat emptor.

rng default

x = [6 10 14 20 26 34 38]';
y = [122 107 119 120 105 95 92]';

numberObservations = numel(x);

NUMBER_RESAMPLE = 5000;
bootstrappedCoefficients = zeros(2,NUMBER_RESAMPLE);

warning("off")
for nr = 1:NUMBER_RESAMPLE
    
    % if mod(nr,100)==0
    %     fprintf("Iteration %d of %d",nr,NUMBER_RESAMPLE)
    % end

    resampleIdx = randsample(numberObservations,numberObservations,true);
    
    xr = x(resampleIdx);
    yr = y(resampleIdx);

    % Tabulate the data
    tbl = table(xr,yr);
    
    % %  Fit the model
    
    % Define function that will be used to fit data
    % (F is a vector of fitting parameters)
    f = @(F,x) F(1).*exp(-x/F(2));
    
    beta0 = [1 1];
    mdl = fitnlm(tbl,f,beta0);
    
    bootstrappedCoefficients(:,nr) = mdl.Coefficients.Estimate;
    
end
warning("on")

figure
tiledlayout(2,1)

nexttile
histogram(bootstrappedCoefficients(1,:))

nexttile
histogram(bootstrappedCoefficients(2,:))

coef1      =    median(bootstrappedCoefficients(1,:))

coef1 = 126.8955

coef1_lo95 = prctile(bootstrappedCoefficients(1,:), 2.5)

coef1_lo95 = 113.1878

coef1_hi95 = prctile(bootstrappedCoefficients(1,:), 97.5)

coef1_hi95 = 145.1811

coef2      =    median(bootstrappedCoefficients(2,:))

coef2 = 130.2354

coef2_lo95 = prctile(bootstrappedCoefficients(2,:), 2.5)

coef2_lo95 = 78.1614

coef2_hi95 = prctile(bootstrappedCoefficients(2,:), 97.5)

coef2_hi95 = 710.0117

[Abhishek K Singh]

Thanks! This is very helpful.