Add constraints to nonlinear multiple data sets fit

Hello everybody,

I am trying to make a fit of several data sets with constraints. I have the code for the fit without the constraints but I am not able to add constraints in the code.

The fitting function is the following:

, where are my fitting parameters. In particular, remains the same for all the data set while are specific for each data set.

Now i want to put some constraints on the values of the fitting parameters and in particular to the global ones.

Attached to the question is reported the code:

clear

x=[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.3318    0.2682    0.2101    0.1831    0.1686    0.1326    0.1333    0.1026    0.1060    0.0787    0.0941    0.0587    0.0415 0.0431    0.0469    0.0455    0.0358    0.0400    0.0336    0.0636    0.0773    0.1207];
y1_err=[0.0127    0.0083    0.0096    0.0067    0.0123    0.0061    0.0101    0.0097    0.0124    0.0074    0.0153    0.0055    0.0069 0.0049    0.0088    0.0062    0.0094    0.0073    0.0163    0.0103    0.0222    0.0355];
y2=[0.1874    0.1506    0.1131    0.0896    0.0765    0.0619    0.0772    0.0474    0.0472    0.0420    0.0533    0.0228    0.0123 0.0195    0.0176    0.0167    0.0085    0.0310    0.0120    0.0480    0.0454    0.0805];
y2_err=[0.0093    0.0063    0.0076    0.0051    0.0098    0.0050    0.0089    0.0086    0.0112    0.0068    0.0140    0.0050    0.0064 0.0045    0.0082    0.0058    0.0089    0.0071    0.0158    0.0099    0.0211    0.0338];

dsid=[ones(size(x)); 2.*ones(size(x))];

T1=400;
T2=420;
T1v = T1*ones(size(x));
T2v = T2*ones(size(x));
xm = x(:)*ones(1,2);
ym = [y1(:) y2(:)];
Tm = [T1v(:) T2v(:) ];
yerr=[y1_err(:) y2_err(:)];
xv = xm(:);
yv = ym(:);
Tv = Tm(:);
yerrv=yerr(:);
weights=1./yerrv;
xTm = [Tv xv dsid];
R=8.314;

B0=[0.2610    0.0523    0.2809    0.7660    0.3668    0.1720    0.4964    0.6377    0.3386    0.9974    0.8835    0.9455  0.9514    0.0231    0.4623]';
%B0=rand(7,1);

[B,Res,J,CovB]=nlinfit(xTm,yv,@subfun2,B0,'Weights',weights);

groups=[1, 2, 3, 4, 5, 6];


figure(1)
gscatter(xv,yv,dsid)
% for i = 1:length(groups)
%     errorbar(xv(dsid==groups(i)),yv(dsid==groups(i)),yerrv(dsid==groups(i)),'LineStyle','None')
% end
line(x,B(1)*exp(-x.^2*B(3)).*( 1-2*exp(B(5)./(R.*T1v)-B(6)./R )./((1+exp(B(5)./(R.*T1v)-B(6)./R )).^2).*(1-sin(x*B(7))./(x*B(7)))),'color','r')
line(x,B(2)*exp(-x.^2*B(4)).*( 1-2*exp(B(5)./(R.*T1v)-B(6)./R )./((1+exp(B(5)./(R.*T1v)-B(6)./R )).^2).*(1-sin(x*B(7))./(x*B(7)))),'color','g')

and the subroutine used:

function yfit = subfun2(params,xTm)
x1=xTm(:,1);
x2=xTm(:,2);
dsid = xTm(:,3);
R=8.314;

A0 = params(1:2); % different A2 for each dataset
A1= params(3:4);   % different A1 for each dataset
A3=params(5);
A4=params(6);
A5=params(7);


yfit=A0(dsid).*exp(-x2.^2.*A1(dsid)).*( 1-2*(exp(A3./(R*x1)-A4/R)./(1+exp(A3./(R*x1)-A4/R)).^2).*(1-sin(x2.*A5)./(x2.*A5)));           


end

[Star Strider]

The nlinfit function will not allow parameter constraints. Use lsqcurvefit for that instead.

Note that it is then possible to use the lsqcurvefit outputs to create statistics on the parameters and the fit with nlparci and nlpredci, although calculating the correct confidence limits and other statistics on a fit using constrained parameters may not be possible with these functions.

I do not understand the ‘global’ reference since I do not see that you are using global variables (please never use them). If you are instead referring to the GlobalSearch function (with respect to your previous Question Fit of multiple data set with variable parameters), it is possible to constrain the parameters there as well.

[Daniele Sonaglioni]

I have tried to put constraints on GlobalSearch function but the fit is ambiguous in the sense that, depending on the input parameters (generated with rand function), the fitted values change. The same holds also for lsqcurvefit, after the imposition of constraints.

Do you have any suggestion? If no, i will accept your answer the same.

[Daniele Sonaglioni]

@Star Strider I am not able to accept your answer. Please, copy the text and answer again to the question.

[Star Strider]

@Daniele Sonaglioni — Done! Thank you!

Answer 1

[Star Strider]

The nlinfit function will not allow parameter constraints. Use lsqcurvefit for that instead.

Note that it is then possible to use the lsqcurvefit outputs to create statistics on the parameters and the fit with nlparci and nlpredci, although calculating the correct confidence limits and other statistics on a fit using constrained parameters may not be possible with these functions.

I do not understand the ‘global’ reference since I do not see that you are using global variables (please never use them). If you are instead referring to the GlobalSearch function (with respect to your previous Question Fit of multiple data set with variable parameters), it is possible to constrain the parameters there as well.

EDIT —

‘I have tried to put constraints on GlobalSearch function but the fit is ambiguous in the sense that, depending on the input parameters (generated with rand function), the fitted values change. The same holds also for lsqcurvefit, after the imposition of constraints.

Do you have any suggestion?’

Nonlinear parameter estimation functions are sensitive to the initial conditions, and can be extremely sensitive to them. That is the reason GlobalSearch worked in your previous Question. I would not use random estimates, instead using a ‘best guess’ approach as to what they would lilkely be. The current problem appeasrs to be similar to the previous one, in that some parameters are stable while others vary with respect to different data sets, and are weighted with respect to both data sets.

To use GlobalSearch with this, use 'fmincon' as the optimiser, since it allows constraints, and the appropriate other arguments to the createOptimProblem call. Also, for constrained problems, the statistics returned by nlparci (and nlpredci) may not be appropriate or accurate.

I do not see a GlobalSearch call in the code example you posted, so I assume you are doing something similar to what we did in the previous problem. My approach would be to use the same techniques on the current problem as I did in the previous problem, with changes required by whatever the differences are between the two problems.

[Daniele Sonaglioni]

Hi,

in the question above I was using a different approach but now I turned back to the one you have suggested me and I am currenctly using 'fmincon' as optimiser, as you seggested.

I noted too that the fit is very sensitive to the starting parameter: i try to use a 'best guess', as you suggested. I hope that following your suggestions the code will work.

Thanks!

[Star Strider]

As always, my pleasure!

I will continue to help as I can.