Non linear fit of multiple data set

1 次查看（过去 30 天）
Daniele Sonaglioni
Daniele Sonaglioni 2021-5-15
回答： Matt J 2021-5-15
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;

请先登录，再回答此问题。

回答（1 个）

Matt J
Matt J 2021-5-15
Because of the constraints, I think you'll just have to do it by Monte Carlo simulation...

类别

在 Help Center 和 File Exchange 中查找有关 Material Sciences 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by