How to fit two equations to two data set simultaneously and with multiple fitting parameters?

Hi everyone,

I want to fit a set of two specific equations (f1 and f2, in the following) with respect to a set of multiple parameters (g1, g2, t1 and t2, in the following). The dataset is given by vectors x, G1 and G2.

Actually, my code works when I consider only g1 and g2 in the formulation for f1 and f2:

x = [100 63.1 39.8 25.1 15.8 10 6.31 3.98 2.51 1.58 1 0.631 0.398 0.251 0.158 0.1];
G1 = [18 15.397 15.014 11.07 7.6413 4.8958 3.1734 2.1804 1.5285 1.0407 0.73853 0.5509 0.39496 0.24756 0.17789 0.11671];
G2 = [914.23 599.24 388.97 251.14 161.92 103.29 66.23 42.402 27.141 17.48 11.303 7.3379 4.8356 3.1736 2.1078 1.4052];

syms g1 t1 g2 t2
f1=((g1.*(t1.*x).^2)./(1+(t1.*x).^2))+0*((g2.*(t2.*x).^2)./(1+(t2.*x).^2));
f2=((g1.*(t1.*x).^1)./(1+(t1.*x).^2))+0*((g2.*(t2.*x).^1)./(1+(t2.*x).^2));
residue1 = sum((f1 - G1).^2);
residue2 = sum((f2 - G2).^2);
residue = residue1 + residue2;
g1_partial = solve(diff(residue,g1),g1)

eqn = subs(residue,g1,g1_partial);
t1_partial = vpasolve(diff(eqn,t1),t1)

t1_partial(abs(imag(t1_partial)) > 1e-30) = [];
t1_partial(((t1_partial)) < 0) = [];
t1_partial = real(t1_partial);
a_full = double(subs(g1_partial, t1_partial))
b_full = double(t1_partial)

G1_pred=subs(f1,g1,a_full);
G1_pred=double(subs(G1_pred,t1,b_full));
G2_pred=subs(f2,g1,a_full);
G2_pred=double(subs(G2_pred,t1,b_full));

figure
hold on
plot(x,G1,'bo-');
plot(x,G1_pred,'ro-');
legend('G1 exp','G1 fit');
hold off

figure
hold on
plot(x,G2,'bo-');
plot(x,G2_pred,'ro-');
legend('G2 exp','G2 fit');
hold off

As You can see, in f1 and f2 I have not taken into account g2 and t2, and the correspondence between fit and expected data is quite low especially on G1.

It is expected that by considering also g2 and t2 the fit would be better.

How can I generalize my code to consider also these two additional parameters?

Thanks in advance,

Alessio

[Mathieu NOE]

sorry to ask you this way, but is there any mistake in the f1 and f2 formulas ?

[Star Strider]

In my version of the code, they both appear to fit their respective data well, using the common parameters.

[Alex Sha]

@Alessio Pricci, there seems to be something wrong with the fitting functions:

f1=((g1.*(t1.*x).^2)./(1+(t1.*x).^2))+0*((g2.*(t2.*x).^2)./(1+(t2.*x).^2));

f2=((g1.*(t1.*x).^1)./(1+(t1.*x).^2))+0*((g2.*(t2.*x).^1)./(1+(t2.*x).^2));

Note the part "0*((g2.*(t2.*x)...", since this part will always be zero due to multiply by 0, so the parameters of "g2" and "t2" will have no effect and can be any values.

Answer 1

[Star Strider]

The Symbolic MathToolbox is not the best option here.

Try a numeric approach, such as this one —

x = [100 63.1 39.8 25.1 15.8 10 6.31 3.98 2.51 1.58 1 0.631 0.398 0.251 0.158 0.1];
G1 = [18 15.397 15.014 11.07 7.6413 4.8958 3.1734 2.1804 1.5285 1.0407 0.73853 0.5509 0.39496 0.24756 0.17789 0.11671];
G2 = [914.23 599.24 388.97 251.14 161.92 103.29 66.23 42.402 27.141 17.48 11.303 7.3379 4.8356 3.1736 2.1078 1.4052];

f1 = @(g1,g2,t1,t2,x) ((g1.*(t1.*x).^2)./(1+(t1.*x).^2))+0*((g2.*(t2.*x).^2)./(1+(t2.*x).^2));
f2 = @(g1,g2,t1,t2,x) ((g1.*(t1.*x).^1)./(1+(t1.*x).^2))+0*((g2.*(t2.*x).^1)./(1+(t2.*x).^2));
f1_fcn = @(b,x) f1(b(1),b(2),b(3),b(4),x);
f2_fcn = @(b,x) f2(b(1),b(2),b(3),b(4),x);
f12_fcn = @(b,x) [f1_fcn(b,x); f2_fcn(b,x)];

format longE
opts = optimset('MaxFunEvals', 1E+6, 'MaxIter',1E+6);
[B,resnorm] = fminsearch(@(b) norm([G1; G2] - f12_fcn(b,x)), randn(4,1), opts)

B = 4×1

1.0e+00 * 3.692495678803561e+12 8.092025209238474e+11 2.549446176516681e-12 -1.848717108306850e+12

resnorm = 
     3.986769725404514e+01

figure
hp1 = plot(x, G1, '.g', 'DisplayName','G1');
hold on
hp2 = plot(x, G2, '.r', 'DisplayName','G2');
hp3 = plot(x, f1_fcn(B,x), '-g', 'DisplayName','G1 Fit');
hp4 = plot(x, f2_fcn(B,x), '-r', 'DisplayName','G2 Fit');
hold off
grid
xlabel('x')
ylabel('G1, G2 Values')
legend('Location','best')

EDIT — (25 Jan 2024 at 12:58)

Added:

format longE

Code unchanged.

EDIT — (25 Jan 2024 at 15:29)

Changing the 0 multipliers to 1 in both and produces this —

f1 = @(g1,g2,t1,t2,x) ((g1.*(t1.*x).^2)./(1+(t1.*x).^2))+1*((g2.*(t2.*x).^2)./(1+(t2.*x).^2));
f2 = @(g1,g2,t1,t2,x) ((g1.*(t1.*x).^1)./(1+(t1.*x).^2))+1*((g2.*(t2.*x).^1)./(1+(t2.*x).^2));
f1_fcn = @(b,x) f1(b(1),b(2),b(3),b(4),x);
f2_fcn = @(b,x) f2(b(1),b(2),b(3),b(4),x);
f12_fcn = @(b,x) [f1_fcn(b,x); f2_fcn(b,x)];

format longE
opts = optimset('MaxFunEvals', 1E+6, 'MaxIter',1E+6);
[B,resnorm] = fminsearch(@(b) norm([G1; G2] - f12_fcn(b,x)), randn(4,1), opts)

B = 4×1

1.0e+00 * -3.924249086830523e+03 5.520395924084962e+02 -1.769994832668654e-03 5.739114582695544e-03

resnorm = 
     1.346916291020836e+01

figure
hp1 = plot(x, G1, '.g', 'DisplayName','G1');
hold on
hp2 = plot(x, G2, '.r', 'DisplayName','G2');
hp3 = plot(x, f1_fcn(B,x), '-g', 'DisplayName','G1 Fit');
hp4 = plot(x, f2_fcn(B,x), '-r', 'DisplayName','G2 Fit');
hold off
grid
xlabel('x')
ylabel('G1, G2 Values')
legend('Location','best')

The fit appears to be reliable, with different estimated parameters. The norm of the residuals is actually less (better fit) when the previously omitted term is included, 13.5 with the entire equaton and 39.9 if the second term is omitted.

.