Nonlinear regression; minimze absolute relative error

Greetings, I would appreciate if anyone can guide me through this. I am trying to fit the following data:

• aw = 0, 0.1145, 0.2260, 0.3232, 0.4355, 0.5251, 0.5779, 0.7083, 0.7554, 0.8432, 0.9044
• X = 0, 0.0149, 0.0353, 0.0565, 0.0590, 0.0754, 0.0906, 0.2688, 0.3326, 0.5300, 0.7894

to the following equation (Eq. 1) using the "fit" function

• Xm*C*K*aw/((1-K*aw)*(1-K*aw+C*K*aw)); parameters K, C, Xm (Eq. 1)

When the absolute relative error (Eq. 2) is minimized with the "Solver" add-in of Microsoft Excel, the fit yields K = 0.9146, C = 8.1415, Xm = 0.0244

Abs((Obs-Pred)/Obs) (Eq. 2)

Nevertheless, I cannot reproduce this K, C, and Xm values using the MATLAB "fit" function. I've already tried to modify other "fit" function arguments, and I ran out of ideas:

• Changing the 'Algorithm'
• Specifying 'Robsut' and 'LAR'
• Setting the Excel solution as the 'StartPoint'
• Constraining the 'Upper' and 'Lower' bounds of every parameter, but paramter Xm always takes the upper bound value
• Setting higher tolerance for model and coefficient values

Thank you,

Vinicio

[Star Strider]

I don’t have the Curve Fitting Toolbox so I can’t help you with that, but I do have nlinfit and lsqcurvefit. With starting parameter estimates of ( [0.1 0.1 0.1] ) and with ‘MaxIter’, 10000, I get:

Xm =  12.9417e+000 
C  =   6.1165e-003 
K  = 809.4168e-003

Mean squared residual = 323.5635e-006

and parameter confidence intervals:

ci =

      12.9417e+000    12.9417e+000
       4.6496e-003     7.5834e-003
     776.2376e-003   842.5960e-003

If you get parameter estimates close to these with fit, I would trust them. I suggest you plot the function with the estimated parameter values against your original data.

[Vinicio Serment]

Well I appreciate your comment, but that's a similar unwanted solution you get when minimizing the least squares using the "fit" function.

Parameter Xm physically represents the "monolayer moisture" which in this case should be close to zero as the independent variable (aw) is cloze to zero, just as in the solution obtained with the Microsoft Excel Solver

Let me try different settings for the two functions you recommend and please let me know if you know anything about minimizing the absolute relative error,

Thanks again, have a nice day

[Matt Tearle]

As Star Strider says, look at the quality of the fit you get from the two solutions:

f = @(C,K,Xm,aw) Xm*C*K*aw./((1-K*aw).*(1-K*aw+C*K*aw));

errorfun = @(C,K,Xm) nanmean(abs((X - f(C,K,Xm,aw))./X));
awfine = linspace(0,1);

K = 0.9146;
C = 8.1415;
Xm = 0.0244;

plot(aw,X,'o',awfine,f(C,K,Xm,awfine))
errorfun(C,K,Xm)

Xm =  12.9417e+000;
C  =   6.1165e-003;
K  = 809.4168e-003;

figure
plot(aw,X,'o',awfine,f(C,K,Xm,awfine))
errorfun(C,K,Xm)

The solution fit is finding has a lower error than the one you're looking for from Excel.

Given your comments about the physical interpretation, it looks like you need to include some kind of constraint on Xm. Are you doing that somehow in Excel? If so, what is the constraint you're imposing?

[Star Strider]

As Matt Tearle notes, if you’re constraining the parameters, we need to know the constraints you’re imposing. We can’t guess what you’re thinking.

I’m obviously missing something. This looks like a reasonably good fit to me:

[Vinicio Serment]

Hello guys, thanks again for trying to help, and sorry that I messed it up greatly by not specifying the possible constrains and posted other values for K, C, and Xm, I'm terribly sorry.

I totally agree with both of you that the regression you both got has a way better fit, but it lacks of physical interpretability. If X = is the moisture of a food system, the theoretical interpretation of Xm corresponds to the layer of water molecules most strongly attached to the food surface, so it cannot be higher than our highest X value.

The values of K, C and Xm in Excel were obtained with good initial estimates after a visual inspection, therefore there was no need not apply any constrains for the Solver solution to minimize the absolute relative error. If constrains were to be applied, Xm = 0.2 could be defined as the upper bound based on the plot of X (yaxis) vs aw (xaxis).

Still, I previously did not know how to anidate a function as Matt Tearle did, so I finally could use the 'fminsearch' function as in the code below and get the results we needed. The lack of fit of the model to this particular data will be part of a discussion.

Once more, thanks for trying to help and for your patience. Have a good weekend.

Vinicio

    % code
aw = [0, 0.1145, 0.2260, 0.3232, 0.4355, 0.5251, 0.5779, 0.7083, 0.7554, 0.8432, 0.9044];
X = [0, 0.0149, 0.0353, 0.0565, 0.0590, 0.0754, 0.0906, 0.2688, 0.3326, 0.5300, 0.7894];
p0 = [1 5 0.03];
% p(1) = K, p(2) = C, p(3) = Xm
f = @(p, aw)(p(1).*p(2).*p(3).*aw./((1-p(1).*aw).*(1-p(1).*aw+p(2).*p(1).*aw)));
errorfun = @(p) nanmean(abs((X-f(p, aw))./X));
opts.FinDiffType = 'central';
opts.MaxIter = 10000;
pars = fminsearch(errorfun, p0, opts)
plot(aw, X, 'db', 'MarkerFaceColor', 'b');
hold all
plot(aw, f(pars, aw), '-^r', 'MarkerFaceColor', 'r')
hold all
legend('Data', 'Fit', 'Location', 'NorthWest')
axis([0 1 0 1])
xlabel('Water activity (a_W)')
ylabel('Moisture (X)')

end

[Star Strider]

When I use your unusual cost function ( errorfun ), I get a value of 5.6621e-003 for the mean squared error (mean of the squared residuals) and parameter estimates:

Xm =  1.0434e+000     
C  =  3.0440e+000    
K  = 44.8857e-003

If I use the sum-of-squares cost function instead:

errorfun = @(p) sum((X-f(p, aw)).^2);

I get a mean squared error of 327.9325e-006 and parameter estimates:

Xm = 847.9914e-003   
C  = 110.4277e-003   
K  = 693.9700e-003

I can’t account for the differences in the parameter estimates between fminsearch and nlinfit with similar mean-squared-errors, other than perhaps nlinfit and the other curve-fitting algorithms are more robust and managed to find the global minimum, and fminsearch found a local minimum. I’ll leave this for those who know more about the details of the respective algorithms.