Non linear fit of multiple data set

Hi everyone,

I am trying to extract the confidence interval of my fitted parameters. I have used a GlobalSearch routine with fmincon as solver. Do you have any suggestion?

Below i attach the code:

clear
%fit per tln+glu under Tg
R=8.314; 
T1=220; %temperature reffered to y1
T2=300; %temperature reffered to y2

yfcn1 = @(b,x) b(1).*exp(-x(:,2).^2.*b(2)).*(1-2*(exp(b(3)./(R.*x(:,1))-b(4)/R)./(1+exp(b(3)./(R.*x(:,1))-b(4)/R)).^2).*(1-sin(x(:,2).*b(5))./(x(:,2)*b(5))));
yfcn2 = @(b,x) b(6).*exp(-x(:,2).^2.*b(7)).*(1-2*(exp(b(3)./(R.*x(:,1))-b(4)/R)./(1+exp(b(3)./(R.*x(:,1))-b(4)/R)).^2).*(1-sin(x(:,2).*b(5))./(x(:,2)*b(5))));
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=[0.9936    0.9375    0.9081    0.8648    0.8568    0.8114    0.7711    0.8010    0.7884    0.7389    0.7901    0.7825  0.7903    0.7501    0.7070    0.7489    0.6441    0.7105    0.6735    0.6385    0.6357    0.6962    0.5946    0.6783]; 
y1_err= [ 0.0637    0.0526    0.0330    0.0235    0.0298    0.0223    0.0388    0.0223    0.0333    0.0326    0.0410    0.0282 0.0561    0.0235    0.0271    0.0218    0.0333    0.0252    0.0344    0.0261    0.0499    0.0396    0.0655    0.0901];
y2=[0.8748    0.8726    0.7922    0.7782  0.7396    0.6958    0.6603    0.6503    0.6556    0.6461    0.6021    0.5820 0.6220    0.5768    0.4950    0.5300    0.5234    0.5170    0.4369    0.4508    0.4409    0.4392    0.4100    0.6699];
y2_err=[ 0.0562    0.0480    0.0287    0.0211    0.0260    0.0194    0.0339    0.0188    0.0287    0.0289    0.0332    0.0225   0.0460    0.0191    0.0211    0.0169    0.0280    0.0198    0.0256    0.0204    0.0392    0.0283    0.0504    0.0856]; 


format long E
T1=220; %temperature reffered to y1
T2=300; %temperature reffered to y2
T1v = T1*ones(size(x));
T2v = T2*ones(size(x));
%T3v = T3*ones(size(x));
xm = x(:)*ones(1,2);
ym = [y1(:) y2(:)];%, y3(:)];
Tm = [T1v(:) T2v(:)];% T3v(:) ];
yerr=[y1_err(:) y2_err(:)];% y3_err(:)];
xv = xm(:);
yv = ym(:);
Tv = Tm(:);
yerrv=yerr(:);
weights=1./yerrv;
xTm = [Tv xv];

%B0 = randn(7,1)*0.1;
B0=[0.5 1e-3 0.5 1e-3 12e4 45 1.5]';

yfcn = @(b,x) yfcn1(b,x).*(x(:,1)==T1) + yfcn2(b,x).*(x(:,1)==T2);
fitfcnw = @(b) norm(weights.*(yv - yfcn(b,xTm)));
lb=[0,0,10e3,0,0,0,0];
ub=[1,1,4e5,1000,2,1,1];
A = @simple_constraint;
problem = createOptimProblem('fmincon', 'x0',B0,'nonlcon',A, 'objective',fitfcnw,'lb',lb,'ub',ub);%
gs = GlobalSearch('PlotFcns',@gsplotbestf);%,'Aineq',ConstraintFunction
[B,fval] = run(gs,problem)
B(:)
% ci = nlparci(B,resid,'jacobian',J1);
format short eng
mdl = fitnlm(xTm,yv,yfcn,B,'Weights',weights)
figure(1)
for k = 1:2
    idx = (1:numel(x))+numel(x)*(k-1);
    subplot(2,1,k)
    errorbar(x.^2, ym(:,k),yerr(:,k), '.')
    hold on
    plot(x.^2, 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)+0.2])
    if k == 1
        text(5, max(yv)+0.1, sprintf('$y = %.3f\\times e^{-x^2\\times %.3f}(1-2\\frac{e^{\\frac{%.3f}{%3d\\times R}-\\frac{%.3f}{R}}}{(1+e^{\\frac{%.3f}{%3d\\times R}-\\frac{%.3f}{R}})^2}(1-\\frac{sin(%.3fx}{%.3fx}))$',B(1:3),T1,B(4),B(3),T1,B(4:5),B(5)), 'Interpreter','latex', 'FontSize',12)
    elseif k == 2
        text(5, max(yv)+0.1, sprintf('$y = %.3f\\times e^{-x^2\\times %.3f}(1-2\\frac{e^{\\frac{%.3f}{%3d\\times R}-\\frac{%.3f}{R}}}{(1+e^{\\frac{%.3f}{%3d\\times R}-\\frac{%.3f}{R}})^2}(1-\\frac{sin(%.3fx}{%.3fx}))$',B(6:7),B(3),T2,B(4),B(3),T2,B(4:5),B(5)), 'Interpreter','latex', 'FontSize',12)
    end
end
xlabel('Q^2')
title('GlobalSearch')

% Tm=[B];
% dlmwrite('C:\Users\Utente\Desktop\neutron_fit\doppia_buca_analisi\genetic_algorithm\prova_params_gs2.txt',Tm,'precision','%.6f');

    



pos = get(gcf, 'Position');
set(gcf, 'Position',pos+[-150 -150 +150 +150])
[ynew,ynewci] = predict(mdl,xTm);
figure(2)
for k = 1:2
    idx = (1:numel(x))+numel(x)*(k-1);
    subplot(2,1,k)
    errorbar(x.^2, ym(:,k),yerr(:,k), '.')
    hold on
    plot(x.^2, ynew(idx), '-r')
    plot(x.^2, ynewci(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('Q^2')
title('GlobalSearch bis')
 

    

and the subroutine used:

function [c, ceq] = simple_constraint(b)
c=[b(4)/8.314-15;
   b(4)-50e4;
   1-b(5)];
ceq=[];

%1e4-b(3)/8.314;

Answer 1

[Matt J]

Because of the constraints, I think you'll just have to do it by Monte Carlo simulation...