Non linear least squares for a system of equations

Hi, I want to estimate 3 parameters using non linear least squares (lsqnonlin) from a system of 3 equations. I know how to write the code for one equation but not sure how to specify for 3 equations?

This is what I am trying to do:

arg Min Σ(Y - X)^2 for 3 Ys and 3 Xs.

Y1 = exp(-x(2)*0.5)*Ydata+(1-exp(-x(2)))*x(1)+(x(3)^2

Y2 = exp(-x(2)*0.1)*Ydata+(1-exp(-x(2)))*x(1)+(x(3)^2

Y3 = exp(-x(2)*0.15)*Ydata+(1-exp(-x(2)))*x(1)+(x(3)^2

X1data = ...

[1,1.5,1.6 etc]

X2data = ...

[2.5, 2,7,2.8 etc]

X3data = ...

[2.7,2.8,3.0 etc]

Ydata = ...

[1.5,1.6,1.7 etc]

For 1 equation (Y1 & X1) I would have written the following code:

fun = @(x)exp(-x(2)*0.5)*Ydata+(1-exp(-x(2)))*x(1)+(x(3)^2-X1data;

x0 = [1.08, 0.4, 0.1];

x = lsqnonlin(fun,x0)

For 3 equations would this be the correct specification?

fun = @(x)exp(-x(2)*0.5)*Ydata+(1-exp(-x(2)))*x(1)+(x(3)^2-X1data+exp(-x(2)*0.1)*Ydata+(1-exp(-x(2)))*x(1)+(x(3)^2-X2data+exp(-x(2)*0.15)*Ydata+(1-exp(-x(2)))*x(1)+(x(3)^2-X3data;

x0 = [1.08, 0.4, 0.1];

x = lsqnonlin(fun,x0)

Many thanks in advance

Answer 1

[Torsten]

fun = @(z)[z(1)^0.5*Ydata+z(2)-X1data,...
           z(1)^0.1*Ydata+z(2)-X2data,...
           z(1)^0.15*Ydata+z(2)-X3data];
z0 = [1 1];
z = lsqnonlin(fun,z0)

where I set

z(1) = exp(-x(2))
z(2) = (1-exp(-x(2)))*x(1)+x(3)^2

So you need to fit in 2 instead of 3 unknown parameters.

[Sofia Philippou]

Thank you! If I change my code to the following (based on your specification above) would it still be correct?

fun = @(x)[exp(-x(2)*0.5)*Ydata+(1-exp(-x(2)))*x(1)+(x(3)^2-X1data,...

exp(-x(2)*0.1)*Ydata+(1-exp(-x(2)))*x(1)+(x(3)^2-X2data,...

exp(-x(2)*0.15)*Ydata+(1-exp(-x(2)))*x(1)+(x(3)^2-X3data];

x0 = [1.08, 0.4, 0.1];

x = lsqnonlin(fun,x0)

[Torsten]

No, because your fitting problem needs 2 instead of 3 parameters (in other words: your formulation has one degree of freedom).

Take a simple example:

If you want to fit (Xdata,Ydata) to a linear function, you would try to determine z(1) and z(2) such that

z(1)*Ydata + z(2) - Xdata

is approximately 0.

What you do is that you try to determine x(1), x(2) and x(3) such that

x(1)*Ydata + x(2) + x(3) - Xdata

is approximately 0.

If the solver determines a solution zsol(1), zsol(2) to the first problem, the solution to the second problem would be

xsol(1) = zsol(1), 
xsol(2) + xsol(3) = zsol(2)

Thus xsol(2) and xsol(3) are not unique - the only condition is that they sum to zsol(2).

[Sofia Philippou]

Thank you for your reply, I see now what the problem is. I did not realise this before and when I asked the question I did not provide the whole equation to keep it simpler. These are my 3 Ys. How can I specify using 2 parameters instead of 3?

Y1 = exp(-x(2)*6/12)*Ydata+(1-exp(-x(2)*6/12))*x(1)+((x(3)^2)/(4*x(2)))*(1-exp(-2*x(2)*6/12))

Y2 = exp(-x(2)*12/12)*Ydata+(1-exp(-x(2)*12/12))*x(1)+((x(3)^2)/(4*x(2)))*(1- exp(-2*x(2)*12/12))

Y3 = exp(-x(2)*3/12)*Ydata+(1-exp(-x(2)*3/12))*x(1)+((x(3)^2)/(4*x(2)))*(1-exp(-2*x(2)*3/12))

Xdata as above

Thank you

[Torsten]

I think here you really need 3 parameters, but I'd determine a, b and c such that

Y1 = a^6*Ydata + (1-a^6)*(b+c*(1+a^6))
Y2 = a^12*Ydata + (1-a^12)*(b+c*(1+a^12))
Y3 = a^3*Ydata + (1-a^3)*(b+c*(1+a^3))

[Sofia Philippou]

Did you mean:

Y1 = a^6/12*Ydata + (1-a^6/12)*b+c*(1-a^2*6/12))

Y2 = a^12/12*Ydata + (1-a^12/12)*b+c*(1-a^2*12/12))

Y3 = a^3/12*Ydata + (1-a^3/12)*b+c*(1-a^2*3/12))

where a = exp(-x(2))

b = x(1)

c = (x(3)^2)/(4*x(2))

This is a silly question but a^6/12 is the same as exp(-x(2))*6/12?

Thank you again

[Torsten]

I meant

a = exp(-x(2)/12)
b = x(1)
c = (x(3)^2)/(4*x(2))

[Sofia Philippou]

Thank you.

So to confirm these are my Ys:

Y1 = a^6*Ydata + (1-a^6)*b+c*(1-a^2*6))

Y2 = a^12*Ydata + (1-a^12)*b+c*(1-a^2*12))

Y3 = a^3*Ydata + (1-a^3)*b+c*(1-a^2*3))

[Torsten]

No, I think

Y1 = a^6*Ydata + (1-a^6)*b + c*(1-a^12)
Y2 = a^12*Ydata + (1-a^12)*b + c*(1-a^24)
Y3 = a^3*Ydata + (1-a^3)*b + c*(1-a^6)

[Sofia Philippou]

Thank you! I run it in MATLAB but the results are my initial points.

Initial point is a local minimum

Any ideas?

[Torsten]

Without knowing your code and your data it's like being in the dark.

[Sofia Philippou]

xdata = ...
 [1.58, 1.58, 1.68, 1.71, 1.56, 1.59, 1.64, 1.73, 1.77, 1.76];
y1data = ...
 [1.91, 1.91, 1.97, 1.93, 1.81, 1.62, 1.66, 1.76, 1.8, 1.80];
y2data = ...
 [1.91, 1.91, 1.94, 1.9, 1.82, 1.82, 1.83, 1.89, 1.93, 1.92];
y3data = ...
 [1.58, 1.6, 1.69, 1.71, 1.57, 1.6, 1.64, 1.74, 1.79, 1.78];
a = exp(-x(2)/12);
b = x(1)
c = (x(3)^2)/(4*x(2));
fun = @(x)[a^6*xdata+(1-a^6)*b+c*(1-a^12) -y1data,...
	   a^12*xdata+(1-a^12)*b+c*(1-a^24)-y2data,...
	   a^3*xdata+(1-a^3)*b+c*(1-a^6)-y3data];
x0 = [2.08, 0.75, 0.3];
x = lsqnonlin(fun,x0)

[Torsten]

xdata =...
 [1.58, 1.58, 1.68, 1.71, 1.56, 1.59, 1.64, 1.73, 1.77, 1.76];
y1data = ...
 [1.91, 1.91, 1.97, 1.93, 1.81, 1.62, 1.66, 1.76, 1.8, 1.80];
y2data =...
 [1.91, 1.91, 1.94, 1.9, 1.82, 1.82, 1.83, 1.89, 1.93, 1.92];
y3data =...
 [1.58, 1.6, 1.69, 1.71, 1.57, 1.6, 1.64, 1.74, 1.79, 1.78];

fun = @(z)[z(1)^6*xdata+(1-z(1)^6)*z(2)+z(3)*(1-z(1)^12) -y1data,...
	   z(1)^12*xdata+(1-z(1)^12)*z(2)+z(3)*(1-z(1)^24)-y2data,...
	   z(1)^3*xdata+(1-z(1)^3)*z(2)+z(3)*(1-z(1)^6)-y3data];

x0 = [2.08, 0.75, 0.3];
z0 = [exp(-x0(2)/12), x0(1), x0(3)^2/(4*x0(2))];   

z = lsqnonlin(fun,z0)

x(1) = z(2);
x(2) = log(z(1))*(-12.0);
x(3) = sqrt(z(3)*4*x(2));

[Sofia Philippou]

Thank you that's perfect! How can I impose conditions to get positive values only? I am also expecting the results to be similar (but not the same) to the initial values.

Many thanks

[Torsten]

lb = [0, 0, 0];
ub = [1, Inf, Inf];
z = lsqnonlin(fun,z0,lb,ub)

[Sofia Philippou]

Sorry to clarify I want to impose conditions on x(1) x(2) and x(3) to be positive only values.

[Torsten]

Yes, these are the settings for z to ensure positive values for x.

[Sofia Philippou]

Ah yes thank you! My results are not what I was expecting but it could be a data problem and not a specification problem? The data in the code are the logs. When I run the code using the original values (without taking the logs) the results are what I am expecting.

[Sofia Philippou]

Apologies, taking the logs or not doesn't make a difference. I do have to take the logs because that is how my model is specified.

[Torsten]

I don't know what you mean. Taking the log of the data was not discussed before.

[Sofia Philippou]

Yes all the data are in logs. Does this change the specification?

[Torsten]

If Y1, Y2, Y3 and Ydata in your three equations from above stand for the logarithmized data, then no - it does not change anything.

[Sofia Philippou]

Yes all the data are logs. Y1, Y2, Y3 and Xdata

[Sofia Philippou]

Hi, not sure if I should post it as a separate question but will the code look like this if I want to estimate the 3 parameters not as a system of equations but as per below?

xdata =...
 [1.39, 1.22, 1.25, 1.35, 1.35];
y1data = ...
 [2.96, 2.90, 2.92, 2.96, 2.97];
y2data =...
 [1.75, 1.77, 1.78, 1.80, 1.83];
y3data =...
 [1.65, 1.69, 1.70, 1.74, 1.75];
y4data =...
 [1.71, 1.68, 1.69, 1.74, 1.75];
y5data =...
 [1.78, 1.845, 1.85, 1.87, 1.88];
lb = [0, 0, 0];
ub = [inf, inf, inf];

fun = @(x)[exp(-x(2)*3/12)*xdata+(1-exp(-x(2)*3/12))*x(1)+((x(3)^2)/(4*x(2)))*(1-exp(-2*x(2)*3/12))-y1data+exp(-x(2)*6/12)*xdata+(1-exp(-x(2)*6/12))*x(1)+((x(3)^2)/(4*x(2)))*(1-exp(-2*x(2)*6/12))-y2data+exp(-x(2)*9/12)*xdata+(1-exp(-x(2)*9/12))*x(1)+((x(3)^2)/(4*x(2)))*(1-exp(-2*x(2)*9/12))-y3data+exp(-x(2)*12/12)*xdata+(1-exp(-x(2)*12/12))*x(1)+((x(3)^2)/(4*x(2)))*(1-exp(-2*x(2)*12/12))-y4data+exp(-x(2)*18/12)*xdata+(1-exp(-x(2)*18/12))*x(1)+((x(3)^2)/(4*x(2)))*(1-exp(-2*x(2)*18/12))-y5data];
x0 = [2.08, 0.75, 0.3];
x = lsqnonlin(fun,x0,lb,ub)