Fit of multiple data sets with error

Hi everybody,

i want to make a fit in Matlab of my experimental data with a known function. I have a code that is able to perform the fit of the data but it does not account for the error on my y-data. This fit is a bit particular because i want to fit three data set with the same function and thus the output fitted parameters are unique for the three cases. Each y-data set corresponds to a different temperature and my fitting function takes into account this fact.

My fitting function is this one:

where . The fitting parameters are .

Now my problem is: how can i make a fit of my data considering the experiemental error too?

Below i show the code with the data and the error associated. If possible, i want that the output parameters have their error.

Kb=8.26e-23;
m=17246;
yfcn = @(b,x) exp(-(x(:,2).*3.*Kb.*x(:,1)/m).^2) .* (1-2*b(1).*b(2).*(1 - sin(x(:,2).*b(3))./(x(:,2).*b(3))));
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=[1.0290    1.0025    1.0158    1.0068    0.9705    0.9646    0.9596    0.9499    0.9811    0.9669    0.9519    0.9573 0.8989    0.9315    0.9618    0.9260    1.0481    0.9245    0.9733    0.8830    1.0203    0.9851    0.9314    0.9204];
y1_err=[0.0637    0.0520    0.0322    0.0239    0.0291    0.0232    0.0405    0.0228    0.0347    0.0364    0.0424    0.0298  0.0531    0.0252    0.0313    0.0230    0.0428    0.0255    0.0375    0.0278    0.0602    0.0467    0.0736    0.0836];
y2=[0.9773    1.0156    0.9362    1.0103    0.9414    0.9167    0.9257    0.9225    0.9127    0.9230    0.9380    0.8899 0.8047    0.9339    0.9012    0.9079    0.9497    0.8449    0.8837    0.8250    0.8738    0.8602    0.9106    0.8656];
y2_err=[0.0616    0.0523    0.0303    0.0239    0.0284    0.0223    0.0394    0.0223    0.0329    0.0351    0.0418    0.0283 0.0493    0.0252    0.0298    0.0226    0.0398    0.0239    0.0350    0.0265    0.0541    0.0425    0.0722    0.0801];
y3=[1.0433    0.9711    0.9913    0.9902    0.9427    0.9146    0.8849    0.9010    0.8876    0.9175    0.9329    0.8639  0.6970    0.8260    0.8675    0.8568    0.8748    0.8156    0.8041    0.7679    0.8333    0.8443    0.8336    0.8344];
y3_err=[0.0638    0.0504    0.0311    0.0232    0.0280    0.0219    0.0375    0.0215    0.0317    0.0343    0.0410    0.0272 0.0445    0.0227    0.0283    0.0213    0.0369    0.0228    0.0323    0.0248    0.0512    0.0405    0.0664    0.0766];
T1=120; %temperature reffered to y1
T2=140; %temperature reffered to y2
T3=160; %temperature reffered to y3
T1v = T1*ones(size(x));
T2v = T2*ones(size(x));
T3v = T3*ones(size(x));
xm = x(:)*ones(1,3);
ym = [y1(:) y2(:), y3(:)];
Tm = [T1v(:) T2v(:) T3v(:)];
xv = xm(:);
yv = ym(:);
Tv = Tm(:);
xTm = [Tv xv];
B0 = rand(3,1);                                                                     % Use Appropriate Initial Estimates
B = fminsearch(@(b) norm(yv - yfcn(b,xTm)), B0);                                    % Estimate Parameters
figure
for k = 1:3
    idx = (1:numel(x))+numel(x)*(k-1);
    subplot(3,1,k)
    plot(x, ym(:,k), '.')
    hold on
    plot(x, 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)])
end
xlabel('x')
sgtitle(sprintf('$y=e^{-(\\frac{3xK_BT}{m})^2} (1-2%.3f\\ %.3f (1-\\frac{sin(x\\ %.3f)}{x\\ %.3f}))$',B,B(3)), 'Interpreter','latex')

Answer 1

[Abhishek]

Hi @Daniele Sonaglioni,

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.