Global fitting of two different function with shared parameters on two different data sets

I am trying to generate a script to fit two different (tediously long) functions that describes two different properties of a certain experimental object. The two functions shares six fitting paremeters and is a function of "T". I have two data sets for each of those respective functions and I am trying to global fit them.

I searched around Matlab answers and found a code that seemingly did what I wanted to do (https://www.mathworks.com/matlabcentral/answers/496168-fitting-2-data-sets-simultaneously-using-two-different-equations-with-some-shared-fit-parameters). And then I made the following code:

syms kB hhat AT1 dG

kB = 0.695034800;

hhat = 5.308959927e-12;

AT1 = 50;

dG = 800;

%Data sets

Adata = [0.48; 0.50; 0.52; 0.7; 0.75; 0.81; 0.82];

Bdata = [8.9e-4; 8.9e-4; 8.8e-4; 8.75e-4; 8.6e-4; 8.5e-4; 7.9e-4];

GFD = [Adata, Bdata];

T = (100:50:400); %T range

%Symbols

%X(1) parameter 1

%X(2) parameter 2

%X(3) parameter 3

%X(4) parameter 4

%X(5) parameter 5

%X(6) parameter 6

%Fitting function

GFF = @(X,T) [X(3)*(X(1)*((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)-dG)^2/(4*X(6)*kB*T)))+X(2)*(AT1+((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)-dG)^2/(4*X(6)*kB*T)))))/((AT1+((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)-dG)^2/(4*X(6)*kB*T))))*(X(3)+X(4)+((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)+dG)^2/(4*X(6)*kB*T))))-((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)+dG)^2/(4*X(6)*kB*T)))*((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)-dG)^2/(4*X(6)*kB*T)))), 2/(AT1+((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)-dG)^2/(4*X(6)*kB*T)))+X(3)+X(4)+((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)+dG)^2/(4*X(6)*kB*T)))-((AT1+((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)-dG)^2/(4*X(6)*kB*T)))-X(3)-X(4)-((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)+dG)^2/(4*X(6)*kB*T))))^2+4*((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)+dG)^2/(4*X(6)*kB*T)))*((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)-dG)^2/(4*X(6)*kB*T))))^(1/2))];

%Root Mean Squared

RMS = @(X) rms(GFD - GFF(X,T));

options = optimset('MaxFunEvals', 1000000, 'MaxIter',1000000, 'Display', 'off', 'TolX', 1e-5);

FIT = fminsearch(RMS,[0.9 0.5 400 400 1e-6 1500],options);

I have two problem here.

1) The stopping criteria (TolX) is different for two functions, but I don't know how to specify that in my code.

2) The bigger problem, is that the it returns an error, "Incorrect dimensions for raising a matrix to a power. Check that the matrix is square and the power is a scalar. To perform elementwise matrix powers, use '.^'."

I am assuming there's a problem in the last fminsearch command with its compatibility with what I'm trying to do, but I have no idea what I need to do.

Could someone help, please? Thank you.

Answer 1

[Walter Roberson]

your T is a vector. Your function involves an expression of T, then ^(1/2) . With T being nonscalar, the ^ operator is the Matrix Power operator, so you are asking for the matrix square root. But matrix power only works for square matrices, not for vectors. You need the .^ operation

[Hayao]

My gosh, how could I forget that. You're right. Thank you. I would also need to do the same about * and /.

But now I got "Arrays have incompatible sizes for this operation". For the line:

RMS = @(X) rms(GFD - GFF(X,T));

I'm assuming this is because I put GFF(X,T), in which T has already been implemented. So I instead tried:

RMS = @(X) rms(GFD - GFF(X));

Then I got "not enough input arguments".

[Hayao]

After twirking things around I got the arrays to match by transposing the T array. It was a row vector, so I needed to change into a column vector. I also removed T off of @(X,T) and did this:

syms kB hhat AT1 dG

kB = 0.695034800;       
hhat = 5.308959927e-12; 
AT1 = 50;             
dG = 800;               

%Data sets
Adata = [0.48; 0.50; 0.52; 0.7; 0.75; 0.81; 0.82];
Bdata = [8.9e-4; 8.9e-4; 8.8e-4; 8.75e-4; 8.6e-4; 8.5e-4; 7.9e-4];
GFD = [Adata, Bdata];
    
T = (100:50:400).';       %T range

%Symbols
%X(1) parameter 1
%X(2) parameter 2
%X(3) parameter 3
%X(4) parameter 4
%X(5) parameter 5
%X(6) parameter 6
%Fitting function
GFF = @(X) [X(3)*(X(1)*((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)-dG)^2/(4*X(6)*kB*T)))+X(2)*(AT1+((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)-dG)^2/(4*X(6)*kB*T)))))/((AT1+((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)-dG)^2/(4*X(6)*kB*T))))*(X(3)+X(4)+((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)+dG)^2/(4*X(6)*kB*T))))-((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)+dG)^2/(4*X(6)*kB*T)))*((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)-dG)^2/(4*X(6)*kB*T)))), 2/(AT1+((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)-dG)^2/(4*X(6)*kB*T)))+X(3)+X(4)+((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)+dG)^2/(4*X(6)*kB*T)))-((AT1+((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)-dG)^2/(4*X(6)*kB*T)))-X(3)-X(4)-((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)+dG)^2/(4*X(6)*kB*T))))^2+4*((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)+dG)^2/(4*X(6)*kB*T)))*((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)-dG)^2/(4*X(6)*kB*T))))^(1/2))];

%Root Mean Squared
RMS = @(X) rms(GFD - GFF(X));
RMS([0.9 0.5 400 400 1e-6 1500])

Error using ^
Incorrect dimensions for raising a matrix to a power. Check that the matrix is square and the power is a scalar. To operate on each element of the matrix individually, use POWER (.^) for elementwise
power.

Error in solution (line 21)
GFF = @(X) [X(3)*(X(1)*((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)-dG)^2/(4*X(6)*kB*T)))+X(2)*(AT1+((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)-dG)^2/(4*X(6)*kB*T)))))/((AT1+((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)-dG)^2/(4*X(6)*kB*T))))*(X(3)+X(4)+((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)+dG)^2/(4*X(6)*kB*T))))-((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)+dG)^2/(4*X(6)*kB*T)))*((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)-dG)^2/(4*X(6)*kB*T)))), 2/(AT1+((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)-dG)^2/(4*X(6)*kB*T)))+X(3)+X(4)+((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)+dG)^2/(4*X(6)*kB*T)))-((AT1+((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)-dG)^2/(4*X(6)*kB*T)))-X(3)-X(4)-((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)+dG)^2/(4*X(6)*kB*T))))^2+4*((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)+dG)^2/(4*X(6)*kB*T)))*((2*pi/hhat)*X(5)*(4*pi*X(6)*kB*T)^(-1/2)*exp(-(X(6)-dG)^2/(4*X(6)*kB*T))))^(1/2))];

Error in solution (line 24)
RMS = @(X) rms(GFD - GFF(X));

options = optimset('MaxFunEvals', 1000000, 'MaxIter',1000000, 'Display', 'off', 'TolX', 1e-5);
FIT = fminsearch(RMS,[0.9 0.5 400 400 1e-6 1500],options);

Now the error is "Unable to perform assignment because the size of the left side is 1-by-1 and the size of the right side is 1-by-2." in the last line.

I don't understand.

[Torsten]

I called your objective function with your vector of initial values for the parameters and it could not be evaluated (see above).

RMS must return a scalar value that usually is the sum of squared differences between your measurement data and the fitted data,

[Hayao]

I see, thanks. I instead used 2-norm. I also need to use constraints for physically valid reason so I used fmincon instead. This is the code:

syms kB hhat AT1 dG

kB = 0.695034800;

hhat = 5.308959927e-12;

AT1 = 50;

dG = 800;

%Data sets

Adata = [0.48; 0.50; 0.52; 0.7; 0.75; 0.81; 0.82];

Bdata = [8.9e-4; 8.9e-4; 8.8e-4; 8.75e-4; 8.6e-4; 8.5e-4; 7.9e-4];

GFD = [Adata, Bdata];

T = (100:50:400).'

%Symbols

%X(1) Parameter 1

%X(2) Parameter 2

%X(3) Parameter 3

%X(4) Parameter 4

%X(5) Parameter 5

%X(6) Parameter 6

%Fitting function

GFF = @(X,T) [X(3).*(X(1).*((2.*pi/hhat).*X(5).*(4.*pi.*X(6).*kB.*T).^(-1/2).*exp(-(X(6)-dG).^2./(4.*X(6).*kB.*T)))+X(2).*(AT1+((2.*pi/hhat).*X(5).*(4.*pi.*X(6).*kB.*T).^(-1/2).*exp(-(X(6)-dG).^2./(4.*X(6).*kB.*T)))))./((AT1+((2.*pi/hhat).*X(5).*(4.*pi.*X(6).*kB.*T).^(-1/2).*exp(-(X(6)-dG).^2./(4.*X(6).*kB.*T)))).*(X(3)+X(4)+((2.*pi/hhat).*X(5).*(4.*pi.*X(6).*kB.*T).^(-1/2).*exp(-(X(6)+dG).^2./(4.*X(6).*kB.*T))))-((2.*pi/hhat).*X(5).*(4.*pi.*X(6).*kB.*T).^(-1/2).*exp(-(X(6)+dG).^2./(4.*X(6).*kB.*T))).*((2.*pi/hhat).*X(5).*(4.*pi.*X(6).*kB.*T).^(-1/2).*exp(-(X(6)-dG).^2./(4.*X(6).*kB.*T)))), 2./(AT1+((2.*pi/hhat).*X(5).*(4.*pi.*X(6).*kB.*T).^(-1/2).*exp(-(X(6)-dG).^2./(4.*X(6).*kB.*T)))+X(3)+X(4)+((2.*pi/hhat).*X(5).*(4.*pi.*X(6).*kB.*T).^(-1/2).*exp(-(X(6)+dG).^2./(4.*X(6).*kB.*T)))-((AT1+((2.*pi/hhat).*X(5).*(4.*pi.*X(6).*kB.*T).^(-1/2).*exp(-(X(6)-dG).^2./(4.*X(6).*kB.*T)))-X(3)-X(4)-((2.*pi/hhat).*X(5).*(4.*pi.*X(6).*kB.*T).^(-1/2).*exp(-(X(6)+dG).^2./(4.*X(6).*kB.*T)))).^2+4.*((2.*pi/hhat).*X(5).*(4.*pi.*X(6).*kB.*T).^(-1/2).*exp(-(X(6)+dG).^2./(4.*X(6).*kB.*T))).*((2.*pi/hhat).*X(5).*(4.*pi.*X(6).*kB.*T).^(-1/2).*exp(-(X(6)-dG).^2./(4.*X(6).*kB.*T)))).^(1/2))];

%Root Mean Squared

SE = @(X) norm(GFD - GFF(X,T));

options = optimset('MaxFunEvals', 1000000, 'MaxIter',1000000, 'Display', 'off', 'TolX', 1e-5);

FIT = fmincon(SE,[0.9 0.1 400 400 1e-6 1500],[],[],[],[],[0 0 0 0 0 0], [1,1, 2000, 2000, 10000, 100000], [], options);

And now it worked. I now have a separate question, so I will post them as a new question.