MATLAB Answers

Fit of multiple data set with variable parameters

2 views (last 30 days)
Hello everyone,
I have a problem in fitting my experimental data. I have a theoretical function of this kind: where A, b, c, d, e are my fitting parameters and R is the gas constant and T is the temperature value. The peculiarity of this function is that A and b varies with temperature.
Now i want ot fit this function to several set of data taken at different temperatures.
The fit i want to make has to be this output: A and b for each temperature and c,d,e global.
I have a code that make the global fit, namely A,b,c,d,e unique fot each temperature and i want to improve it.
Do you have any suggestion?
Below i report the code i am using:
clear
R=8.314;
yfcn= @(b,x) b(1)*exp(-x(:,2).^2.*b(2)).*(1-2*exp(b(3)./(R.*x(:,1))-b(5)/R).*(1-sin(x(:,2).*b(4))./(x(:,2)*b(4))));
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];
T1=220; %temperature reffered to y1
T2=300; %temperature reffered to y2
%T3=320; %temperature reffered to y3
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 =[0.2941 0.5306 0.8324 0.5975 0.3353];%rand(5,1); % Use Appropriate Initial Estimates 0.6596
%B1 = fminsearch(@(b) norm(sqrt(weights).*(yv - yfcn(b,xTm))), B0); % Estimate Parameters
[B,R,J,CovB] = nlinfit(xTm,yv,yfcn,B0,'Weights',weights);
% X=B(1);
% Y=B(2);
% Z=B(3); 0.7640 0.8182 0.1002 0.1781 0.3596
% H=B(4); 0.0567 0.5219 0.3358 0.1757 0.2089
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)])
end
xlabel('Q^2')
%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')

Accepted Answer

Star Strider
Star Strider on 1 Apr 2021
(I was in the middle of posting an Answer here when my desktop crashed. Again. I really hate Windows!)
That aside, the problem needs a better set of initial parameter estimates, and the GlobalSearch funciton (in the Global Optimization Toolbox) provides it:
fitfcnw = @(b) norm(weights.*(yv - yfcn(b,xTm)));
problem = createOptimProblem('fmincon', 'x0',B0, 'objective',fitfcnw);
gs = GlobalSearch('PlotFcns',@gsplotbestf);
[B,fval] = run(gs,problem)
with the result:
GlobalSearch stopped because it analyzed all the trial points.
All 26 local solver runs converged with a positive local solver exit flag.
B =
8.862485547533210e-01
1.509391962458558e-02
-1.870848919343233e+04
1.088251234643560e+00
-4.504651478455145e+01
fval =
7.223670210423392e+00
(The code before and after this remains the same.)
that with:
mdl = fitnlm(xTm,yv,yfcn,B,'Weights',weights)
(that I added simply to get some statistics), produces:
Nonlinear regression model:
y ~ b1*exp( - x2^2*b2)*(1 - 2*exp(b3/(R*x1) - b5/R)*(1 - sin(x2*b4)/(x2*b4)))
Estimated Coefficients:
Estimate SE tStat pValue
________ _________ _______ __________
b1 0.90447 0.018777 48.169 4.7384e-39
b2 0.013486 0.0039251 3.436 0.0013206
b3 -11214 5454.7 -2.0558 0.045905
b4 1.1185 0.094728 11.807 4.4142e-15
b5 -20.751 17.026 -1.2188 0.22956
Number of observations: 48, Error degrees of freedom: 43
Root Mean Squared Error: 0.215
R-Squared: 0.932, Adjusted R-Squared 0.926
F-statistic vs. zero model: 3.28e+03, p-value = 2.38e-54
demonstrating an excellent fit, and this plot:
This appears to be an improvement over the previous result using only nonlinfit.
Plot the results of the fitnlm with:
[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')
if you so desire.
  6 Comments
Star Strider
Star Strider on 13 Jul 2021
My pleasure!
I am quite pleased that my code is useful in similar situations, and is apparently robust!
For 20 temperatures, my initial suggestion would be to set parts of it up in a loop, however the way I wrote this does not easily lend itself to that approach, since there are only two sets of data vectors and hard-coding it was appropriate.
A for loop might be appropriate for the initial set-up, for example creating ‘T1v’, ‘T2v’, ‘ym’, ‘Tm’ and ‘y_err’ (that could be made even easier if all of them were matrices), and the other loops could be expanded to accommodate the additional temperatures, instead of only two. Otherwise, the only option is to hard-code the initial set-up steps.
I hope it will be sufficiently robust to deal with 20 temperatures.
If you encounter problems with those changes, please post them as a new Question (with your code and all the necessary data, preferably as text or spreadsheet files) and provide a link to it here. I will follow up as best I can to help you get it to work.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!