Nonlinear regression with nested models and two or more independent variables

Hello everyone, I would appreciate if anyone could help me.

I'm trying to find parameters a0, a1, n0, and n1 that best describe my response variable (z) as a function of two independent variables (x, y).

My equation set is:

• z(x, y) = -a(y) * x ^ [n(y)]
• a(y) = ln{1 + exp[a0 * (y - a1)]}
• n(y) = n0*exp(-n1 * y)

My dataset is:

• x = [0 1 2 3 4 12 24 0 1 2 3 4 8 16 24 0 1 2 3 4 0 0.5 1 2 3 4];
• y = [100 100 100 100 100 100 100 200 200 200 200 200 200 200 200 300 300 300 300 300 400 400 400 400 400 400];
• z = [0 0.59 0.76 0.33 -0.2 -0.28 -0.3 0 -1.27 -0.11 -0.3 -1.03 -1.21 -1.55 -2.6 0 -0.58 -0.73 -2.48 -2.56 0 -0.5 -3.4 -3.69 -4.78 -3.96];

I already tried the following code to define my z(x, y) function:

a = @(ap, y)(log(1 + exp(-ap(1)*(y - ap(2)))))
n = @(np, y)(np(1)*exp(-np(2)*y))
zp = {a, n}
z = @(zp, x)(-zp(1)*x^(zp(2)))

Apparently, the function z(x, y) is correctly defined. I also tried the following code to define parameter vectos "a" and "n" to evaluate the z(x, y), but it did not work and I ignore how to do this.

ap0 = [150 2];
np0 = [2 0.1];
S(pars, t)

After evaluating the function, I would appreciate if someone taught how to use the "fminsearch" function when using two independent variables and two different parameter arryas (a, n)

Thanks you for your time,

Vinicio

Answer 1

[Star Strider]

The regression functions characteristically accept only two input arguments: the parameter vector and the independent variable vector (or array). All the parameters have to be expressed as elements of the parameter vector. When you have more than one independent variable, combine them both in one matrix and then reference them by column (my convention) within the function. Here, I combined both in ‘xy’. The first comment line maps your parameters and independent variables to the ones I used in my function. The second comment line is your original function with the substitutions I made. I was able to express the entire equation as an anonymous function. I used fminsearch because you requested it. It is appropriate for small parameter vector sizes, and since it uses the derivative-free approach for its oprimisation, is robust.

This works:

% b(1) = a0,  b(2) = a1,  b(3) = n0,  b(4) = n1,  x = xy(:,1),  y = xy(:,2)
% z = @(b, xy)  -log(1 + exp[a0 .* (y - a1)]) .* x.^(n0*exp(-n1 * y));
zf = @(b, xy)  -log(1 + exp(b(1) .* (xy(:,2) - b(2)))) .* xy(:,1).^(b(3)*exp(-b(4) * xy(:,2)));

x = [0 1 2 3 4 12 24 0 1 2 3 4 8 16 24 0 1 2 3 4 0 0.5 1 2 3 4];
y = [100 100 100 100 100 100 100 200 200 200 200 200 200 200 200 300 300 300 300 300 400 400 400 400 400 400];

xy = [x' y'];
z = [0 0.59 0.76 0.33 -0.2 -0.28 -0.3 0 -1.27 -0.11 -0.3 -1.03 -1.21 -1.55 -2.6 0 -0.58 -0.73 -2.48 -2.56 0 -0.5 -3.4 -3.69 -4.78 -3.96]';

OLSCF = @(b) sum((z - zf(b,xy)).^2);     % Ordinary Least Squares Cost Function

[B, fv] = fminsearch(OLSCF, [1; 1; 1; 1]*0.01);

figure(1)
stem3(x, y, z, 'filled', 'b')
hold on
stem3(x, y, zf(B,xy), 'filled', 'r')
hold off
grid on
legend('Data', 'Fit')

It gives parameter estimates:

a0 = 0.0062
a1 = 0.023
b0 = 0.081
b1 = 0.012

I plotted your data and the fit using stem3 because it’s easier to see how the variables compare on that plot.

I asked fminsearch to output both the estimated parameter vector ‘B’ and the final function value, ‘fv’, that here corresponds to the sum squared error.

[Vinicio Serment]

Thank you Star Strider, that really worked well!

You mentioned "I used fminsearch because you requested it". It looks that there are better options in your opinion. Which of the other minimization functions do you recommend?

By the way, can this same approach be used for an *.m file function instead of an anonymous defined function?

Thanks again Star Strider,

Vinicio

[Star Strider]

My pleasure!

Other options specifically for curve fitting are the Statistics Toolbox function nlinfit (also fitnlm) and the Optimization Toolbox function lsqcurvefit. For small problems (a maximum of about 7 parameters), fminsearch works well. It just requires an extra line to define the cost function (in my code ‘OLSCF’).

An .m-file function file would work as well. The cost function changes to:

OLSCF = @(b) sum((z - zf(b,xy)).^2);

with the function file defined as:

function y = zf( b, xy )
    y = -log(1 + exp(b(1) .* (xy(:,2) - b(2)))) .* xy(:,1).^(b(3)*exp(-b(4) * xy(:,2)));
end

For the nlinfit and lsqcurvefit, the call would be:

B0 = ones(4,1)*0.01;
B = lsqcurvefit(@zf, B0, xy, z);

and for nlinfit:

B = nlinfit(xy, z, @zf, B0);

(I tested all of these to be certain they work.)