Nonlinear regression with 2 independent variables

I want to fit a nonlinear model to a set of experimental data. It has 1 dependent variable, i, and 2 independent variables, td and Tr. I have 7 values for Tr, 7 for td and, therefore, 49 for i. The model I want to fit is this:

i = k*Tr^m/((c+td)^n).

 So, I want to determine the values of the parameters k, m, c and n. I've been working on this for days and could not find a solution. Could anyone help me out?
Thanks in advance.
Marcus.

[dpb]

So, what have you tried and where did you run into problems? Which Toolboxes do you have (Curve Fitting, Optimization, none, ... ???)

[John D'Errico]

Or the stats toolbox?

[Marcus Souza]

I have Curve Fitting and Optimization toolboxes, but didn't use them. I tried to use fitnlm function to adjust the model, but I think my problem is probably related on how to dispose the variables. I have i as a 7x7 matrix, Tr and td as a 7x1 array (each one).

Answer 1

[Star Strider]

This is actually straightforward if not all that well-documented. You have to combine your independent variables into a matrix, then refer to the matrix columns within your regression function as the individual independent variables. This will work with fminsearch, nlinfit, and lsqcurvefit. (I don’t have the Curve Fitting Toolbox, so I can’t comment on it.) Note that here ‘i’ also has to be a column vector, so if you use nlinfit or lsqcurvefit, use it as ‘i(:)’ as I did in my code here.

Using fminsearch (core MATLAB function, so no Toolboxes required):

% Original Equation: i = k*Tr^m/((c+td)^n)
% MAPPING: b(1) = k,  b(2) = m,  b(3) = c,  b(4) = n
i_fcn = @(b,x) b(1).*x(:,2).^b(2) ./ (b(3) + x(:,1)).^b(4);         % Regression Equation
x = [td(:)  Tr(:)];                                                 % Independent Variable Matrix
b0 = rand(4,1);                                                     % Initial Parameter Estimates
SSECF = @(b) sum((i(:) - i_fcn(b,x)).^2);                           % Sum-Squared-Error Cost Function
[b_estd, SSE] = fminsearch(SSECF, b0);                              % Estimate Parameters
i_fit = i_fcn(b_estd,x);                                            % Generate Fit Line

I tested this with random data and it runs without error. It should give you the result you want.

[Marcus Souza]

Thanks!! It worked very well. Just one more thing: if I wanted to use nlinfit, how I would have to do?

[Star Strider]

My pleasure!

For nlinfit, you would use my ‘i_fcn’ function and ‘x’ matrix without modification:

[beta,R,J,CovB,MSE,ErrorModelInfo] = nlinfit(x, i(:), i_fcn, b0);

The fminsearch function is actually more robust than nlinfit since fminsearch uses a derivative-free method to estimate its parameters. It would be best to run the fminsearch code first, then use those parameter estimates as your ‘b0’ vector for nlinfit.

[Marcus Souza]

I tried to do this, but it comes with the error MODELFUN must be a function that returns a vector of fitted values the same size as Y (49-by-1). Which array I should change to solve this?

[Star Strider]

I’m not sure what the problem is, since I don’t have your data to work with. To use the parameter estimates from the fminsearch with nlinfit, the nlinfit call with the rest of the code becomes:

% % Original Equation: i = k*Tr^m/((c+td)^n)
% % MAPPING: b(1) = k,  b(2) = m,  b(3) = c,  b(4) = n
i_fcn = @(b,x) b(1).*x(:,2).^b(2) ./ (b(3) + x(:,1)).^b(4);         % Regression Equation
x = [td(:)  Tr(:)];                                                 % Independent Variable Matrix
b0 = rand(4,1);                                                     % Initial Parameter Estimates
SSECF = @(b) sum((i(:) - i_fcn(b,x)).^2);                           % Sum-Squared-Error Cost Function
[b_estd, SSE] = fminsearch(SSECF, b0);                              % Estimate Parameters
i_fit = i_fcn(b_estd,x);                                            % Generate Fit Line
[beta,R,J,CovB,MSE,ErrorModelInfo] = nlinfit(x, i(:), i_fcn, b_estd);
i_fit2 = i_fcn(beta,x);                                             % Generate Fit Line

[Marcus Souza]

Yes, it still gives the same error for the nlinfit. My data is:

i = [112.8353  145.7318  166.5839  178.0802  186.0354  192.1177  210.6576   14.9915   19.3622   22.1326   23.6600   24.7170   25.5251   27.9883   19.3381   24.9760   28.5498   30.5200   31.8834   32.9258   36.1033   22.3542   28.8714   33.0025   35.2800   36.8561   38.0611   41.7341   25.1041   32.4230   37.0623   39.6200   41.3899   42.7432   46.8680   27.5879   35.6310   40.7292   43.5400   45.4851   46.9722   51.5051   38.1440   49.2647   56.3137   60.2001   62.8893   64.9454   71.2129]';
Tr =  [2     5    10    15    20    25    50];
td = [5    10    15    20    25    30    60];

[dpb]

"MODELFUN must be a function that returns a vector of fitted values the same size as Y (49-by-1)."

I'd guess perhaps your output vector is 1x49 instead of 49x1? Check orientation, if that's it, just transpose (.'). You need to give details, not just description including such things as what the pertinent shapes actually are. This will have the side benefit of likely leading you to the problem while you're developing the information for us to try to evaluate on your own...after all, we can't see your terminal from here, only what you post.

[Star Strider]

They all have to be the same lengths. If you have measured empirical data, you should have as many dependent as independent variable observations. Since I have no idea what you’re doing, you need to resolve this.

One possible way is:

Trv =  [2     5    10    15    20    25    50];
tdv = [5    10    15    20    25    30    60];
Tr = linspace(min(Trv), max(Trv), length(i));
td = linspace(min(tdv), max(tdv), length(i));

That solves the dimension inequality problem. You need to decide if it’s what you want to do.

[Marcus Souza]

Let me explain better the problem I need to solve: I have two independent variables (Tr and td), each one of them has 7 values. I stored them in the vectors:

Tr =  [2     5    10    15    20    25    50];
td = [5    10    15    20    25    30    60];

And I have 1 dependent value (i). Each pair of value (Tr, td) leads to a value of i. Therefore, having 7 values of Tr and 7 values of td, I have 49 values of i. Initially, I stored i in matrix form:

i=[112.8353 145.7318 166.5839 178.0802 186.0354 192.1177 210.6576 ; 14.9915 19.3622 22.1326 23.6600 24.7170 25.5251 27.9883 ; 19.3381 24.9760 28.5498 30.5200 31.8834 32.9258 36.1033 ; 22.3542 28.8714 33.0025 35.2800 36.8561 38.0611 41.7341 ; 25.1041 32.4230 37.0623 39.6200 41.3899 42.7432 46.8680 ; 27.5879 35.6310 40.7292 43.5400 45.4851 46.9722 51.5051 ; 38.1440 49.2647 56.3137 60.2001 62.8893 64.9454 71.2129]

So, in the form i is organized above, each column corresponds to a different td and each line corresponds to a different Tr. For example, if Tr = 2 and td = 5, i = 112.8353 ; if Tr = 50 and td = 60, i = 71.2129. To predict this, I want to fit this distribution (preferentially using nlinfit), to the following model:

i = k*Tr^m/((c+td)^n)

So, I want to determine the parameters k, m, c and n that best fit the model. For example, k =1000, m = 0.9, c = 15 and n = 1.1, which gives me i = 1000*Tr^0.9/((15+td))^1.1. With this model, I can predict how much i will be when, for example, Tr = 200 and td = 120, simply substituting into the equation. I'm not a Matlab expert nor statistic expert, so I apologize if I was not clear before. Thanks in advance!

[Star Strider]

It would've been quite helpful to know that earlier.

I’m not certain if ‘Tr’ defines the rows and ‘td’ the columns, or if I have those reversed. Either way, it would be easy to modify this code to make them work, by simply transposing your ‘i’ matrix in the fit routines. The surface plot of your function and the fitted values suggests a decent fit:

[Trm,tdm] = meshgrid(Tr, td);
x = [tdm(:)  Trm(:)];                                               % Independent Variable Matrix
b0 = rand(4,1);                                                     % Initial Parameter Estimates
SSECF = @(b) sum((i(:) - i_fcn(b,x)).^2);                           % Sum-Squared-Error Cost Function
[b_estd, SSE] = fminsearch(SSECF, b0);                              % Estimate Parameters
i_fit = i_fcn(b_estd,x);                                            % Generate Fit Line
[beta,R,J,CovB,MSE,ErrorModelInfo] = nlinfit(x, i(:), i_fcn, b_estd);
i_fit2 = i_fcn(beta,x);                                             % Generate Fit Line

i_mtx = reshape(i, 7, 7)';
i_fit2_mtx = reshape(i_fit2, 7, 7)';

figure(1)
surf(Trm, tdm, i_mtx)
hold on
stem3(Trm, tdm, i_fit2_mtx)
hold off
grid on
xlabel('T_r')
ylabel('t_m')
zlabel('\iti\rm (\itT_r,t_d\rm)')

[Anand Shirke]

Dear Star Strider,

I have a similar problem where I have to estimate the Constants. I tried your same code and am facing some errors. I am attaching the code and error picture.

Please help.

Trv =  [2     5    10    15    20    25    50];
tdv = [5    10    15    20    25    30    60];
i=[112.8353 145.7318 166.5839 178.0802 186.0354 192.1177 210.6576 ; 14.9915 19.3622 22.1326 23.6600 24.7170 25.5251 27.9883 ; 19.3381 24.9760 28.5498 30.5200 31.8834 32.9258 36.1033 ; 22.3542 28.8714 33.0025 35.2800 36.8561 38.0611 41.7341 ; 25.1041 32.4230 37.0623 39.6200 41.3899 42.7432 46.8680 ; 27.5879 35.6310 40.7292 43.5400 45.4851 46.9722 51.5051 ; 38.1440 49.2647 56.3137 60.2001 62.8893 64.9454 71.2129];

Tr = linspace(min(Trv), max(Trv), length(i));
td = linspace(min(tdv), max(tdv), length(i));
i_fcn = @(b,x) b(1).*x(:,2).^b(2) ./ (b(3) + x(:,1)).^b(4);
x = [td(:)  Tr(:)];  
b0 = rand(4,1);
SSECF = @(b) sum((i(:) - i_fcn(b,x)).^2);
[b_estd, SSE] = fminsearch(SSECF, b0);
i_fit = i_fcn(b_estd,x);