I understand that you are trying to fit a known analytical function to three datasets, each corresponding to a different temperature, and you want to ensure the fit shares a common set of parameters across all three. While your current approach fits the data well, it does not take into account the experimental uncertainties in the 'y' values. I went ahead and updated your fitting code to account for the experimental errors in your ‘y’ data by using weighted least squares. This gives more weight to data points with less uncertainty, which seems to be exactly what you need.
Instead of minimizing the plain norm, I have modified the cost function to include weights based on the inverse of the squared errors.
Please try out the below steps:
- Use weighted least squares: Since you have errors associated with your ‘y’ data, you can assign weights as the inverse of the squared errors.
weights = 1 ./ (y_err_all.^2);
- Update the objective function: Instead of using the ‘norm’ method, define a custom function that calculates the weighted sum of squared residuals. You can replace this line:
B = fminsearch(@(b) norm(yv - yfcn(b,xTm)), B0);
With this line:
weighted_residual = @(b) sum(((y_all - yfcn(b,xT)).^2) .* weights);
[B, resnorm] = fminsearch(weighted_residual, B0);
- Parameter uncertainties calculation using Jacobian: After finding the best-fit parameters, estimate the Jacobian matrix using finite differences, and calculate the covariance matrix. This will help us get the final error. Add the below code, after the fitting code:
delta = 1e-6;
J = zeros(length(y_all), length(B));
for i = 1:length(B)
Bp = B;
Bm = B;
Bp(i) = B(i) + delta;
Bm(i) = B(i) - delta;
J(:,i) = (yfcn(Bp,xT) - yfcn(Bm,xT)) / (2*delta);
end
W = diag(weights);
CovB = inv(J' * W * J);
param_errors = sqrt(diag(CovB));
- Finally plot the errors: Replace your plot logic with the below code:
errorbar(x, y_all(idx), y_err_all(idx), 'k.', 'MarkerSize', 10)
hold on
plot(x, yfcn(B,xT(idx,:)), '-r', 'LineWidth', 1.5)
With this, you now get a fit that shows the experimental uncertainties, and the output parameters come with estimated confidence intervals. I tried this approach on MATLAB R2024b, and I got the following plot as well as the experimental errors for each of the fitted parameters.
Hope this will help you.
