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

9 次查看（过去 30 天）
Hayao
Hayao 2023-2-1
评论： Hayao 2023-2-2
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.

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

采纳的回答

Walter Roberson
Walter Roberson 2023-2-1
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
  4 个评论
显示 2更早的评论
Torsten
Torsten 2023-2-1
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
Hayao 2023-2-2
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.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Conversion Between Symbolic and Numeric 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by