Determining the uncertainty of the derivative of a fit

Hi, I have a dataset of values x, for which I am computing an exponential fit exp_fit = fit(x, y, 'exp2'). This return a 1x1 cfit, of which I am able to compute the confidence intervals using predint(exp_fit, x,0.95,'functional','on').

Now want to differentiate the fit, for which I am using differentiate(exp_fit, x)'. This returns a double, for which I am unable to extarct the confidence inervals.

Does anyone know how to obtain the confidence intervals of the derivative of a fit?

I have included an example with the census data:

load census

exp_fit = fit(cdate,pop,'exp1');
p12 = predint(exp_fit,cdate,0.95,'functional','on');
dy_exp_fit = differentiate(exp_fit, cdate)';


fig = figure('Visible', 'off');
fig.Position = [360.3333   50.3333  650.6667  670.6667];

ax1 = subplot(2,1,1); % Define ax1 for the first subplot
yyaxis(ax1, 'left');
plot(ax1, exp_fit, cdate, pop)
hold on;
plot(ax1, cdate, p12)
hold on;
plot(ax1, cdate, dy_exp_fit)
legend


fig = figure('Visible', 'on');

Thank you in advance!

Answer 1

[Hitesh]

Hi Elena,

You need to consider how uncertainty in the fit parameters affects the derivative for obtaining the confidence intervals for the derivative of a fitted curve. However, you can approximate these intervals using a Monte Carlo simulation approach.Kindly follow the following steps to obtain the confidence intervals of the derivative :

• Parameter Uncertainty: You need to obtain the standard errors of the fit parameters. This can be done using the "confint" function, which provides confidence intervals for the parameters.
• Monte Carlo Simulation: Then generate a large number of random samples of the fit parameters, assuming they are normally distributed around the estimated values with standard deviations equal to their standard errors.
• Differentiate Each Sample: After that for each sample of parameters, compute the derivative of the fitted function.
• Compute Confidence Intervals: Finally, calculate the confidence intervals from the distribution of the derivatives obtained in the previous step.

Refer the below code as an example:

load census
exp_fit = fit(cdate, pop, 'exp1');

% Get parameter estimates and their standard errors
coeffs = coeffvalues(exp_fit);
conf = confint(exp_fit);
param_std_err = (conf(2,:) - conf(1,:)) / (2 * 1.96); % Assuming 95% CI

% Monte Carlo simulation
num_samples = 1000;
param_samples = zeros(num_samples, length(coeffs));
for i = 1:length(coeffs)
    param_samples(:, i) = normrnd(coeffs(i), param_std_err(i), num_samples, 1);
end

% Compute derivatives for each sample
dy_samples = zeros(num_samples, length(cdate));
for i = 1:num_samples
    % Calculate the derivative for each sample
    dy_samples(i, :) = param_samples(i, 1) * param_samples(i, 2) * exp(param_samples(i, 2) * cdate);
end

% Compute confidence intervals for derivatives
dy_lower = prctile(dy_samples, 2.5);
dy_upper = prctile(dy_samples, 97.5);

% Compute the derivative of the original fit
dy_exp_fit = differentiate(exp_fit, cdate)';

% Plot results
figure;
subplot(2,1,1);
yyaxis left;
plot(cdate, pop, 'o');
hold on;
plot(cdate, exp_fit(cdate), '-');
plot(cdate, dy_exp_fit, '-r', 'DisplayName', 'Derivative');
plot(cdate, dy_lower, '--r', 'DisplayName', 'Derivative Lower CI');
plot(cdate, dy_upper, '--r', 'DisplayName', 'Derivative Upper CI');
legend('Data', 'Fit', 'Derivative', 'Derivative Lower CI', 'Derivative Upper CI');
title('Exponential Fit and Derivative with Confidence Intervals');

For more information on "Monte Carlo Simulation", kindly refer to the following MATLAB documentation:

https://www.mathworks.com/discovery/monte-carlo-simulation.html

https://www.mathworks.com/matlabcentral/fileexchange/157501-monte-carlo-simulation

[Elena]

Hi Hitesh, thank you so much for your detailed answer - as I am still going through, may I ask where the '1.96' comes from in the param_std_error?

[Hitesh]

Hi Elena,

The factor 1.96 is used because of the properties of the normal distribution.The critical value for a standard normal distribution is approximately 1.96 for a 95% confidence interval. This value represents the number of standard deviations from the mean required to encompass 95% of the probability mass of a standard normal distribution.

param_std_err = (conf(2,:) - conf(1,:)) / (2 * 1.96);

In the above code snippet, the confidence intervals for the parameter estimates are extracted as follows:

• conf(2,:) and conf(1,:) are the upper and lower bounds of the 95% confidence interval for the parameters, respectively.
• (conf(2,:) - conf(1,:)) calculates the width of the confidence interval.
• Dividing by 2 * 1.96 gives the standard error of the parameter estimates. The factor of 2 accounts for the fact that the confidence interval is symmetric around the estimate, and dividing by 1.96 converts the interval width to a standard error assuming a normal distribution.