Nonlinear regression with nested models and two or more independent variables
10 views (last 30 days)
Show older comments
Vinicio Serment
on 5 Dec 2014
Commented: Star Strider
on 8 Dec 2014
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
0 Comments
Accepted Answer
Star Strider
on 5 Dec 2014
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.
2 Comments
Star Strider
on 8 Dec 2014
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.)
More Answers (0)
See Also
Categories
Find more on Descriptive Statistics in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!