Fit of multiple data sets

Hello everyone,

I have a problem in making a fit of some set of data. I have 4 arrays of data: the first is my x-axis while the other three are my y-axes. I want to fit the y-arrays with the same function but i want only one set of parameters to fit the three arrays.

My fitting function is the subsequent:

where . My fitting parameters are and each y-array is referred to a different value of T.

How can i write the code for this kind of fit?

Answer 1

[Star Strider]

That is relatively straightforward. If your objective function outputs only one vector (as opposed to three), concatenate the dependent variables in one vector, create a matching independent variable vector, then use whatever nonlinear paramerer estimation routine you want to estimate the parameters.

It would help to have your data and the and m values.

[Daniele Sonaglioni]

is the Blotzmann constant and equal to while .

If i understand well: you propose to merge all my y-arrays in one macro-array and fit the whole array with an x-array of the same length?

[Star Strider]

Yes.

Concatenate all the ‘y’ values in one vector and all the ‘T’ values in a matching vector (The individual ‘T’ and ‘y’ vectors do not have to be the same lengths as the others, so long as the individual ‘T’ and ‘y’ vector dimensions match).

Then do the regressions. To plot them, simulate the ‘y’ vectors with the individual ‘T’ vectors in each plot.

I did the same sort of approach recently fitting the results to differential equations in Parameter estimation for a system of differential equations with multiple time spans

[Daniele Sonaglioni]

I have tried to follow the procedure suggested by @William Rose but i'm not sure that the code works correctly.

In fact i have defined a matrix in which the second colums is an array of the temperature but i am not able to plot the fitting function for a single temperature.

Can you give me some suggestion? Below i report my code and my x,y1,y2,y3 arrays.

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];
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];
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];
T1=120; %temperature reffered to y1
T2=140; %temperature reffered to y2
T3=160; %temperature reffered to y3

A=[x' y1' y2' y3'];
N=length(x);
B=[A(:,1),ones(N,1)*T1,A(:,2); A(:,1),ones(N,1)*T2,A(:,3); A(:,1),ones(N,1)*T3,A(:,4)];
K=8.26e-23;
m=17246;

x1=B(:,1);
y=B(:,3);
t=B(:,2);

fun=@(r,x1,t) exp(-3*K*x1.*t./m).*(1-2*r(2)*r(3)*(1-sin(x1*r(4)))./(x1*r(4)));
r0=[1,1,1,1];
r=lsqcurvefit(fun,r0,x1,y);

[William Rose]

@Daniele Sonaglioni

In a meeting now. I will look at it tonight.

[Star Strider]

Try this:

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];
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];
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];
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')

producing:

EDIT — (31 Apr 2021 at 1:07)

Updated sgtitle and resulting plot. Code otherwise unchanged.

[William Rose]

@Daniele Sonaglioni,

In your original post you wrote

where u=2*kB*T/m

but your code has

fun=@(r,x1,t) exp(-3*K*x1.*t./m).*(...)

in other words, your code does not square the stuff inside the exp(). Which is correct?

[Star Strider]

My ‘yfcn’ squares everything inside the exponent, although the sgtitle should be corrected to:

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')

to emphasize that.

[William Rose]

@Daniele Sonaglioni,

My code is very similar to @Star Strider, but mine is less elegant.

I did not square the quantity xu=3x*kB*T/m inside the exp(), since @Daniele Sonaglioni did not, in his code, and because I have seen in physical models, but I have not seen . But my experience is limited.

The quantity -x*u, or -(x*u)^2, inside the exp() is the only part of the predicted y that depends on temperature. The magnitude of x*u is on the order of 10^-24, or 10^-48 if you square it, at all three temperatures and for all values of x in this dataset. Therefore exp(-x*u)=1, and exp(-(x*u)^2)=1. Therefore it doesn't matter if you square x*u, and the predictions are independent of temperature. Thus you can replace exp(-x*u), or exp(-(x*u)^2), with 1, and you can discard the temperature information, because it is irrelevant.

I suspect that the temperature information is not actually irrelevant. I suspect that the model equation or the value of one of the constants is incorrect. I would check the value of m.

I think there is another problem. It looks like the predicted values depend only the product of p1 and p2, and are unaffected by replacing (p1,p2) with (p1/c, p2*c), where c is any non-zero constant. Therefore their individual values cannot be estimated with this model. Only their product can be estimated.

[William Rose]

@Daniele Sonaglioni, when I say "It looks like the predicted values depend only the product of p1 and p2", I mean that that is what your original formula y=... shows, and that is how the function in your code works.

[Star Strider]

I used the expression provided in the Question, and did not look further.

If that is not the correct expression, it would likely be straightforward to change it.

Note the my code incorporates ‘T’ as the first column of the ‘xTm’ matrix, rather than as a separate vector, and then refers to it and ‘x’ as the appropriate columns of ‘xTm’ in ‘yfun’.

[Daniele Sonaglioni]

@Star Strider @William Rose your code is wonderfuld and it works greatly! In the code i have showed to you i have forgotten to square the exponent but, as @William Rose noted, it is quasi irrelevant.

Anyway i have only more question: how can i account fot the y-error in my fit? Because the best thing is to make a weighted fit.

[Star Strider]

Thank you!

I do not understand ‘account fot the y-error in my fit’. Do you want confidence intervals on the fit, or something else? That would likely require changing from fminsearch to the Statistics and Machine Learning Toolbox fitnlm function. That is relatively straightforward so not a problem, however you will have to have that Toolbox to run my code.

Also, what weighting values do you want to use for the weighted fit? That is also relatively straightforward to do, however I need the weighting vector.

[Daniele Sonaglioni]

I have solved on my own this problem by switching from fminsearch to nlinfit. What i mean with y-error is the experimental uncertainty of the mesurements and what i need is the confidence interval for the output parameters.

With nlinfit i can accomplish to this aim by passing to it the weigths that are equal to 1/(error)^2 and the confidence interval is equal to sqrt(diag(Covariance_matrix)).

Anyway, thank you!

[Daniele Sonaglioni]

I think that the two methods (mine and your) are equivalent.

[Star Strider]

With nlinfit, use nlparci to get the confidence intervals on the parameters and nlpredci to get the confidence intervals on the fit. The fitnlm function calculates the parameter confidence intervals automatically, and the predict function calculates the confidence intervals on the fit.

Answer 2

[William Rose]

@Daniele Sonaglioni,

@Star Strider is right. Specifically, if your data is in array A, and if A has N rows x 4 columns, and the columns are x, y1, y2, y3, and the 3 temps are T1, T2, T3, then do

B=[A(:,1),ones(N,1)*T1,A(:,2); A(:,1),ones(N,1)*T2,A(:,3); A(:,1),ones(N,1)*T3,A(:,4)]; 

The above command turns your Nx4 array into 3*N by 3 array. Note the use of commas and semicolons in the line above. I inserted a column for temperature. The first N rows have temp T1, the next N rows have temp T2, etc. Then you fit it all at once as @Star Strider said.