Multiple variables in non linear regression

Hi,

Just downloaded Matlab with the ambition of trying to fit an equation i have to my data through adding curve fitting parameters.

The problem:

I have 3 independent vairables x1, x2, x3 and my output variable i am trying to fit them to is y

thus the equation looks like y=x1*x2*x3 and i want to add the following parameters A,m,n so the equation will looks roughly like

y=A*x1*(x2^m)*(x3^n) in order to fit my data

i know i cant use the curve fitting tool and am in need of help

Thanks in advance!

i am also using the 1 month trial version

[Sebastian Daneli]

Try Rstudio, it's a free software for statistics, very popular, works like Matlab. You also have to download R, the programming language.

https://rstudio.com/

https://www.r-project.org/

[Andrew Doyle]

is it straightforward to use? as i am already confused enough

[Sebastian Daneli]

No, not really. Excel is on the other hand much "easier", at least somewhat more "straightforward"

[Andrew Doyle]

i've tried looking for excel help but the equations generally look like y= ax+bx

and arent to the power or for multiple independent variables

Answer 1

[Star Strider]

Try this:

x1 = rand(1,10);                                                % Independent Variable Vector
x2 = rand(1,10);                                                % Independent Variable Vector
x3 = rand(1,10);                                                % Independent Variable Vector
y  = rand(1,10);                                                % Dependent Variable Vector

xm = [x1(:) x2(:), x3(:)];                                      % X Matrix
yfcn = @(b,x) b(1).*x(:,1).*(x(:,2).^b(2)).*(x(:,3).^b(3));     % Objective Function

b0 = rand(3,1);
B = lsqcurvefit(yfcn,b0,xm,y(:))

That runs without error.

Make appropriate changes in the initial parameter estimates vector ‘b0’ to be compatible with your data.

EDIT —

I did not see that you had posted your data while I was writing my Answer.

Try this (with ‘y’ shortened by one element to be compatible with ‘xm’):

y = y(1:end-1);

xm = [x1(:) x2(:), x3(:)];                                      % X Matrix
yfcn = @(b,x) b(1).*x(:,1).*(x(:,2).^b(2)).*(x(:,3).^b(3));     % Objective Function

b0 = rand(3,1);
B = lsqcurvefit(yfcn,b0,xm,y(:))

producing:

B =
          0.21033007874519
         -1.94721481237717
         0.123415136896506

Plotting it is not possible (would require four dimensions) so you need to determine if these are reasonable. If not, experiment with different values for the elements in the ‘b0’ vector to get the appropriate parameter values.

That consideration aside, and since ‘x2’ is constant, one way to plot it would be to leave out ‘x2’ from the plot and then select the view parameters to plot it as a 2D plot:

figure
plot3(x1, x3, y, 'p')
hold on
plot3(x1, x3, yfcn(B,xm), '-r')
hold off
grid on
xlabel('x_1')
ylabel('x_3')
view(0,90)

producing:

This appears to be a reasonably good fit!

[Andrew Doyle]

i made a mistake and inputted some of the data wrong, i have corrected it and used your script but this appears

[Andrew Doyle]

x1 = [1.36475E-05

1.78359E-05

2.00991E-05

3.15843E-05

3.32467E-05

1.33488E-05

2.30703E-05]; % Independent Variable Vector

x2 = [130000

130000

130000

130000

130000

130000

130000]; % Independent Variable Vector

x3 = [26.1

26.1

24.8

24.8

24.8

45.6

45.6]; % Independent Variable Vector

y = [10.8

11.9

11

14

13.3

15.9

13.8

13.4]; % Dependent Variable Vector

xm = [x1(:) x2(:), x3(:)]; % X Matrix

yfcn = @(b,x) b(1).*x(:,1).*(x(:,2).^b(2)).*(x(:,3).^b(3)); % Objective Function

b0 = rand(3,1);

B = lsqcurvefit(yfcn,b0,xm,y(:))

[Andrew Doyle]

>> untitled3

Error using lsqcurvefit (line 271)

Function value and YDATA sizes are not equal.

Error in untitled3 (line 33)

B = lsqcurvefit(yfcn,b0,xm,y(:))

[Andrew Doyle]

for x1, i just forgot to do 1/the values but the correct data is there now

[Star Strider]

Please see my edited Answer using your data (likely posted while you were posting your Comment). It is necessary to shorten ‘y’ by one element so that it has the same row size as ‘xm’.

Then my code works, and appears tror give a decent result!

[Andrew Doyle]

ill try this now thanks!

[Star Strider]

My pleasure!

[Andrew Doyle]

get the same error: basically this is my problem i am trying to solve

Basically my shear values for my 1st set of data is

tmax = 10.8 , 11.9 , 11 ,14 , 13.3 ,15.9 ,13.8

and i need an ewquation to predict tmax, an equation already exists which is t(f) = P / ( pi * d * L)

where d (diamater) and L (length)

but to predict the max shear the equation needs to be in terms of just the material properties.

So i am using the same equation but substituting P = E(frp) * fc

(Elastic Modulus of frp) and (compressive strength of concrete) as they have the biggest impact on the pull out out, (pulling out an frp rod out of a cylinder block)

thus t(f) = ( Efrp * fc ) / ( pi * d * L )

when i do this i get the results

t(f) = 500.10 , 720.15 , 712.7, 1425.5, 1425.5 , 1258.19, 1887.29

so if i then do t(f) = ( Efrp * fc ) / ( pi * d * L * 100)

you get t(f) = 5.00 , 7.20 , 7.13 , 14.26 , 14.26 ,12.58 , 18.87

which as you can see are fairly close to the actual values hence why i am trying to add the curve fitting parameters Efrp^m and fc^n and using A instead of 100 to get a value that is more suited

the other values are as follow:

E = 130000 ,130000 ,130000 ,130000 , 130000 ,130000 , 130000

This is only a constant for this set of data, it varies more in the other FRP rods

fc = 26.1 , 26.1 , 24.8 , 24.8 , 24.8 ,45.6 , 45.6

d = 12 ,10 ,12 , 12 ,12 ,10 , 10

l = 180 , 150 ,120 ,60 ,60 , 150 , 100

pi = 3.141

apologoies about the long message trying to be as detailed to help you guys

[Star Strider]

‘get the same error’

I am using the equation you specified as the objective function, and the data you provided.

I have no idea what you are doing otherwise. The data must have the same dimensions so that when combined in my ‘xm’ matrix (as I wrote it), it and the ‘y(:)’ vector have the same row length. I have no control over that without having the data you are actually using.

My code works with the data you posted, after editing the variables to be certain the row lengths were the same. That is the best I can do.

[Andrew Doyle]

ill just keep trying but thankyou so much! appreciate your time

[Star Strider]

My pleasure.

If you want to attach your data, I will do my best to fit it. I do not need to know the other details, my background being medicine and biomedical engineering, not mechanical engineeering or material science.

[Andrew Doyle]

Hi again, i have the script running as follow and have re written the function so that it goes

y = a*( ( (x1^m)*(x2^n) ) / ( x3 ) )

and have inputted all the data again correctly and for some reason the data doesnt match the y data the script i have wrote

[Andrew Doyle]

x1 = [130000

130000

130000

130000

130000

130000

130000

130000

120000

144000

155000

155000

155000

155000

155000

155000

155000

155000

130000

130000

130000

130000

130000

124000

124000

124000

124000

124000]

x2=[26.1

26.1

24.8

24.8

24.8

45.6

45.6

45.6

20.7

20.7

42.7

42.7

42.7

42.7

42.7

42.7

42.7

42.7

60

60

60

60

60

60

60

60

60

60]

x3=[6784.56

4711.5

4523.04

2261.52

2261.52

4711.5

3141

1570.5

4276.03176

5704.3701

4276.03176

4276.03176

4276.03176

1435.31136

2123.06472

2840.7204

3558.37608

8522.1612

353.3625

706.725

1060.0875

1413.45

1766.8125

282.69

565.38

848.07

1130.76

1413.45]

y=[10.8

11.9

11

14

13.3

15.9

13.8

13.4

8.4

7.5

18.8

22.3

18.4

29.9

26.9

22.3

20.1

15.1

36

32

28

24

25

33

30

27

23

21]

y = 2.*((x1.^3).*(x2.^2)./(x3.')); % y vectors with A=2, m=3, n=2

xdata = [x1 x2 x3];

ydata = [y]

model = @(A, m, n, x1, x2, x3) A.*((x1.^m).*(x2.^n)./(x3.'));

model_lsq = @(param, xdata) model(param(1), param(2), param(3), ...

xdata(:,1), xdata(:,2), xdata(:,3));

sol = lsqcurvefit(model_lsq, rand(1,3), xdata, y)

A = sol(1);

m = sol(2);

n = sol(3);

[Andrew Doyle]

and this is the answers i get

sol =

0.5678 0.0759 0.0540

but if you do

0.5678 * ( (x1^0.0759 * x2^0.054) / x3 )

the answers do not match up with y

[Star Strider]

It will likely not be a perfect match. I do not know if it is necessary to constrain the parameters (for example to be positive). I used the ga function with an unconstrained problem and with 100 different runs, and the best parameter set was:

7.437845749605835   1.070000000000000  -1.195894462198460

with a residual norm of:

54.741895805282986

with the function and residual norm calculated as:

xm = [x1(:) x2(:), x3(:)];                                          % X Matrix
yfcn = @(b,x) b(1) .* (x(:,1).^b(2) .* x(:,2).^b(3)) ./ x(:,3);     % Objective Function

resnorm = @(theta) norm(y-yfcn(theta,xm));

Constrained to force all the parameters to be , the parameters and residual norm are:

   0.006956686973572   1.252512910246849   0.000000071048737
  61.819773481310285

.

Answer 2

[Ameer Hamza]

If you have the optimization toolbox, then you can use lsqcurvefit. The following shows an example using random data

x1 = rand(100,1);
x2 = rand(100,1);
x3 = rand(100,1);
y = 2.*x1.*(x2.^3).*(x3.^2); % random y vectors with A=2, m=3, n=2
xdata = [x1 x2 x3];

model = @(A, m, n, x1, x2, x3) A.*x1.*(x2.^m).*(x3.^n);
model_lsq = @(param, xdata) model(param(1), param(2), param(3), ...
    xdata(:,1), xdata(:,2), xdata(:,3));

sol = lsqcurvefit(model_lsq, rand(1,3), xdata, y)
A = sol(1);
m = sol(2);
n = sol(3);

Result

>> A
A =
    2.0000
>> m
m =
    3.0000
>> n
n =
    2.0000

[Andrew Doyle]

How would i input the data as i wont be usng random data.

would it simply be

x1 = [1, 2, 3 ]

etc?

Thankyou!

[Ameer Hamza]

Yes, you can just type or paste the values. However, these must be column vector, so you need to put transpose operator (.') at the end like this

x1 = [1, 2, 3, ....].';
x2 = [9, 8, 7, ....].';

[Andrew Doyle]

the script has worked however the equation does not link to the y values i need?

would i not need to input the y values so the script knows to make A, m, n ewual to the correct constants?

[Ameer Hamza]

How did you put y-values? Can you show your values of x1, x2, x3, and y.

[Andrew Doyle]

x1 = [73273.248
56066.85
49753.44
31661.28
30078.216
74912.85
43345.8];
x2 = [130000
130000
130000
130000
130000
130000
130000];
x3 = [26.1
26.1
24.8
24.8
24.8
45.6
45.6];
y=[10.8
11.9
11
14
13.3
15.9
13.8
13.4];
 
y = 2.*x1.*(x2.^3).*(x3.^2); % y vectors with A=2, m=3, n=2
xdata = [x1 x2 x3 y];
model = @(A, m, n, x1, x2, x3) A.*x1.*(x2.^m).*(x3.^n);
model_lsq = @(param, xdata) model(param(1), param(2), param(3), ...
    xdata(:,1), xdata(:,2), xdata(:,3));
sol = lsqcurvefit(model_lsq, rand(1,3), xdata, y)
A = sol(1);
m = sol(2);
n = sol(3);

[Andrew Doyle]

basically my equation is

t = ( E * f ) / ( pi * d * l )

but only E and f need to be to the power and i wanted to add a constant to help thus

t = A* ( ( E^m * f^n) / pi * d * l ) )

so i thought clumping variables together may be best hecne y=a*x1*x2^m*x3^n

where x1= 1 / pi * d * l )

do you think it would be easier to just put the whole eqaution in?

[Andrew Doyle]

as i think the way i have tried to clump the variables together may be the issue

[Andrew Doyle]

could be possible to write it as

y = A ( x1^m * x2^n ) / (x3)

where x1 - E, x2- f and x3 - pi*d*l

[Andrew Doyle]

tried the above with the script i sent you and this is the error

Error using lsqcurvefit (line 271)

Function value and YDATA sizes are not equal.

[Ameer Hamza]

Can you exactly paste the code you are trying to run?

[Andrew Doyle]

x1 = [73273.248
56066.85
49753.44
31661.28
30078.216
74912.85
43345.8];
x2 = [130000
130000
130000
130000
130000
130000
130000];
x3 = [26.1
26.1
24.8
24.8
24.8
45.6
45.6];
y=[10.8
11.9
11
14
13.3
15.9
13.8
13.4];
 
y = 2.*x1.*(x2.^3).*(x3.^2); % y vectors with A=2, m=3, n=2
xdata = [x1 x2 x3 y];
model = @(A, m, n, x1, x2, x3) A.*x1.*(x2.^m).*(x3.^n);
model_lsq = @(param, xdata) model(param(1), param(2), param(3), ...
    xdata(:,1), xdata(:,2), xdata(:,3));
sol = lsqcurvefit(model_lsq, rand(1,3), xdata, y)
A = sol(1);
m = sol(2);
n = sol(3);

[Andrew Doyle]

>> untitled2
Initial point is a local minimum.
Optimization completed because the size of the gradient at the initial point 
is less than the value of the optimality tolerance.

<stopping criteria details>
sol =
    0.2769    0.0462    0.0971

[Andrew Doyle]

but when i use them values for A, m, n i get a y value of 9.65x10(-5)

which is way off as you can see in my y data the first y should be 10.8

[Ameer Hamza]

Yes. The message indicates that you got a local minimum. That is why you get wrong 'y' values from the model. This is happening due to the very different scaling of x1, x2, and x3. Can you share the actual values of independent variables in this equation: t = A* (( E^m * f^n) / pi * d * l )). That will help in suggesting a suitable solution.

[the cyclist]

x2 is a constant vector, which means that x2.^m is a constant vector, so that term is just multiplying by a constant factor. So, don't need that term.

Also, your y vector has 8 elements, and your x vector has 7 elements. Not sure what to do about that.

[Andrew Doyle]

x2 is a constant for this short set of data, but i have alot more data where x2 varies i am just trying to make a script that works before going crazy with lots of data.

Basically my shear values for my 1st set of data is

tmax = 10.8 , 11.9 , 11 ,14 , 13.3 ,15.9 ,13.8

and i need an ewquation to predict tmax, an equation already exists which is t(f) = P / ( pi * d * L)

where d (diamater) and L (length)

but to predict the max shear the equation needs to be in terms of just the material properties.

So i am using the same equation but substituting P = E(frp) * fc

(Elastic Modulus of frp) and (compressive strength of concrete) as they have the biggest impact on the pull out out, (pulling out an frp rod out of a cylinder block)

thus t(f) = ( Efrp * fc ) / ( pi * d * L )

when i do this i get the results

t(f) = 500.10 , 720.15 , 712.7, 1425.5, 1425.5 , 1258.19, 1887.29

so if i then do t(f) = ( Efrp * fc ) / ( pi * d * L * 100)

you get t(f) = 5.00 , 7.20 , 7.13 , 14.26 , 14.26 ,12.58 , 18.87

which as you can see are fairly close to the actual values hence why i am trying to add the curve fitting parameters Efrp^m and fc^n and using A instead of 100 to get a value that is more suited

the other values are as follow:

E = 130000 ,130000 ,130000 ,130000 , 130000 ,130000 , 130000

This is only a constant for this set of data, it varies more in the other FRP rods

fc = 26.1 , 26.1 , 24.8 , 24.8 , 24.8 ,45.6 , 45.6

d = 12 ,10 ,12 , 12 ,12 ,10 , 10

l = 180 , 150 ,120 ,60 ,60 , 150 , 100

pi = 3.141

apologoies about the long message trying to be as detailed to help you guys