fitting implicit non-linear equation

Dear ALL, I have the data points below, defined by the variables X1,X2,X3,X4:

X1              X2      X3      X4
3.5E-4        0  0.99965  1
4.08935E-4  0.00138  0.99821  0.85588
0.00114          0.03473  0.96413  0.30777
6.43932E-4  0.00989  0.98947  0.54354
0.00145          0.04977  0.94878  0.24193
0.0019          0.08999  0.90811  0.18414
0.00268          0.14986  0.84745  0.13042
0.00268          0.14998  0.84734  0.13056
0.00465          0.2488  0.74655  0.07521

I need to fit these data points using a function of the type

a*b*ln(X4)=W1*X2*(1-X1)+W2*X3*(1-X1)-W3*X2*X3+W4*X2*X3*(1-2*X1)

where a,b are constants (8.31 and 973), and W1,W2,W3,W4 are adjustable parameters. Could you indicate me how to do that? Thank you in advance, Matthieu

Answer 1

[Torsten]

Use lsqnonlin with f(i) defined as

f(i) = a*b*ln(X4(i))-(W1*X2(i)*(1-X1(i))+W2*X3*(1-X1(i))-W3*X2(i)*X3(i)+W4*X2(i)*X3(i)*(1-2*X1(i)))

Best wishes

Torsten.

[MATT]

Hello, and thank. I am not sure I understand how to input the indices:

I define my function as:

function F = FIT(W, X1, X2, X3, X4)
i=1:9;
F = 8.31*973*ln(X4(i))-(W(1)*X2(i)*(1-X1(i))+W(2)*X3(i)*(1-X1(i))-W(3)*X2(i)*X3(i)+W(4)*X2(i)*X3(i)*(1-2*X1(i)));
end

Then I launch the following command:

W0=[0.10,1,40,1];
[W] = lsqnonlin(@FIT,W0, X1, X2, X3, X4);

This does not seem to work, I guess it inputs all of the X(i) rather than fitting the data. Any advice on how to debbug it? Best, Matthieu

[Torsten]

As I see now, your function is linear in the unknown Parameters.

So you can determine W simply by

A=zeros(9,4);

B=zeros(9);

for ii=1:9

A(ii,1)=X2(ii)*(1-X1(ii));

A(ii,2)=X3(ii)*(1-X1(ii));

A(ii,3)=-X2(ii)*X3(ii);

A(ii,4)=X2(ii)*X3(ii)*(1-2*X1(ii));

B(ii)=a*b*log(X4(ii));

end

W=A\B;

Best wishes

Torsten.

Answer 2

[arich82]

I question the physicality of your fitting equation:

Since you have a log on the left-hand side, and your data includes X4(1) == 1 and X2(1) == 0 exactly, it seems like W2 should necessarily be 0 for the equation to remain valid...

In your data, it seems like X3 is very nearly (1 - X2). If you use this substitution, you can plot your data (and your fit) with plot3.

Feel free to play with the code below. I used the backslash operator for fitting, which I believe results in a least-squares fit via QR (but don't hold me to that):

First, some preliminary definitions.

%%%%
%%data
%%%%
clear; close all; clc;
data = [...
  0.00035,        0,         0.99965,   1; ...
  0.000408935,    0.00138,   0.99821,   0.85588; ... % swapped rows
  0.000643932,    0.00989,   0.98947,   0.54354; ... % swapped rows
  0.00114,        0.03473,   0.96413,   0.30777; ...
  0.00145,        0.04977,   0.94878,   0.24193; ...
  0.0019,         0.08999,   0.90811,   0.18414; ...
  0.00268,        0.14986,   0.84745,   0.13042; ...
  0.00268,        0.14998,   0.84734,   0.13056; ...
  0.00465,        0.2488,   0.74655,   0.07521; ...
];
% Note: swapping rows 3 & 4 dowsn't affect the fitting, but makes
% ploting lines cleaner

X1 = data(:, 1);
X2 = data(:, 2);
X3 = data(:, 3);
X4 = data(:, 4);
a = 8.31;
b = 973;

% re-parameterize (unnecessary, but might simplify some later algebra)
T = a*b*log(X4);
X = X2;
Y = 1 - X1;
Z = X3;
% Note: Z = X3 ~~ (1 - X2) == (1 - X)

For the first fit, we'll assume X3 = 1 - X1

%%%%
%%w: fit assuming X3 = (1 - X2), using all four w's
%%%%
% T = w2*Y + (w1 - w2 + 2*w4)*X.*Y + (w3 + w4)*(X.^2 - X) - 2*w4*X.^2.*Y;

M = [Y, X.*Y, X.^2 - X, X.^2.*Y];
m = M\T;
Tm = M*m;
X4m = exp(Tm/a/b);

w(4) = m(4)/2;
w(3) = m(3) - w(4);
w(2) = m(1);
w(1) = m(2) + w(2) - 2*w(4);
w = w(:);
err_m = (X4 - X4m)./X4m;
w
err_m

Now let's see what happens if we enforce the (perhaps more physical) constraint W2 = 0

%%%%
%%w2: fit assuming X3 = (1 - X2), using w2 = 0
%%%%
% T = (w1 + 2*w4)*X.*Y + (w3 + w4)*(X.^2 - X) - 2*w4*X.^2.*Y;
M = [X.*Y, X.^2 - X, X.^2.*Y];
m2 = M\T;
Tm2 = M*m2;
X4m2 = exp(Tm2/a/b);

w2(4) = m(3)/2;
w2(3) = m(2) - w2(4);
w2(2) = 0;
w2(1) = m(2) - 2*w2(4);
w2 = w2(:);
err_m2 = (X4 - X4m2)./X4m2;
w2
err_m2

Now lets use the full equation (no assumption on X3)

%%%%
%%W: fit using all four w's (no assumptions)
%%%%
% 
% TW = ...
%       W(1)*X2.*(1 - X1) ...
%     + W(2)*X3.*(1 - X1) ...
%     - W(3)*X2.*X3 ...
%     + W(4)*X2.*X3.*(1 - 2*X1);

M = [X2.*(1 - X1), X3.*(1 - X1), -X2.*X3, X2.*X3.*(1 - 2*X1)];
W = M\(a*b*log(X4));
TW = M*W;
X4W = exp(TW/a/b);
err_W = (X4 - X4W)./X4W;
W
err_W

And finally, using W2 == 0 (no assumption on X3)

%%%%
%%W: fit using w2 = 0 (no additional assumptions)
%%%%
M = [X2.*(1 - X1), -X2.*X3, X2.*X3.*(1 - 2*X1)];
W2 = M\(a*b*log(X4));
TW2 = M*W2;
X4W2 = exp(TW2/a/b);
err_W2 = (X4 - X4W2)./X4W2;
W2 = [W2(1), 0, W2(2), W2(3)].';
W2
err_W2

results

%%%%
%%compare results
%%%%
compare_W = [w, w2, W, W2];
compare_W

compare_err = 100*[err_m, err_m2, err_W, err_W2];
compare_err

hf = figure('WindowStyle', 'docked');
ha = axes;
hp = plot3(X1, X2, X4, X1, X2, X4m,  X1, X2, X4m2,  X1, X2, X4W,  X1, X2, X4W2)
xlabel('X1')
ylabel('X2')
zlabel('X4')
legend('X4', 'X4m', 'X4m2', 'X4W', 'X4W2')
grid on;
title('compare fits')

You can use your judgment to see if this %error is reason able over the region of interest.

 >> compare_W =

   1.0e+08 *

   3.813510730697298  -0.000491562357117   0.025845576531701   0.032553401777266
  -0.000018277710575                   0  -0.000018555355233                   0
   2.867011340403814   0.959509492664202   1.017723304845835   1.273754444212918
  -0.947009230361176   0.960001055021319   0.991441774206119   1.240781494391939

 >> compare_err =

   25.3541         0   25.7751         0
    9.5807  -12.2780    9.8547  -12.3620
  -18.5631  -32.9458  -18.6725  -33.2918
  -19.8587  -26.6898  -20.2522  -27.2932
   -5.6106   -6.3826   -5.7995   -6.7560
   19.3484   26.9156   19.2132   26.7386
    0.7148    2.1690    0.9361    2.5826
    0.2742    1.5628    0.5416    2.0342
   -1.9625   -3.4741   -2.0987   -3.7206