Comparing two curve fits (using AIC?)

Hello all,

Let's say for a moment I have some a priori ideas about a family of functions that might best describe a particular data set. After I fit the data with each candidate, I can simply look at the outputs (e.g., RMSE, r²) and chose the one with the best values. However, is there a way to gain some degree of confidence about that selection? For example, let's say I have reason to believe that in a set of data like the below, the true distribution is described by a low-order integer exponent in the function y=b×xⁿ+c. So, I might capture some data and then try to fit it as follows:

% Generate some data
b = 1; % True coefficient
n = 3; % True exponent
c = 0; % True y-intercept
rx = randi(100, [100, 1])/10;
rn = 3 + rand([100, 1])/4 - 0.125;
rc = (rand([100, 1]) * 500) - 250 + c;
ry = b * rx.^(rn) + rc;

% Generate and fit the models
m = @(b, c, n, x) b * x.^n + c;
m1 = @(p, x) m(p(1), p(2), 1, x); % n=1
m2 = @(p, x) m(p(1), p(2), 2, x); % n=2
m3 = @(p, x) m(p(1), p(2), 3, x); % n=3

c1 = lsqcurvefit(m1, [1, 0], rx, ry)
c2 = lsqcurvefit(m2, [1, 0], rx, ry)
c3 = lsqcurvefit(m3, [1, 0], rx, ry)

% Plot the results
x = linspace(min(rx), max(rx));
plot(rx, ry, 'ok', x, m1(c1, x), '-g', x, m2(c2, x), '-r', x, m3(c3, x), '-b')

% Determine which model is most likely to be correct.  Use AIC?

In the above example, m3 should generally fit the data best, unless the random number generation is very unlucky, because that's what I set it to.

For those who are familiar with Prism, there I might perform this test by starting an "analysis" and chosing compare, and then selecting "for each data set, which of two equations (models) fits best" and then chosing the "Akaike's Information Criterion" test. After running that on one set of data, I got

This is very helpful as it lets me know how certain I am about the fit choice (i.e., 86% sure n=3 is better than n=2). What's the best way to do something similar in MATLAB?

Thanks in advance.

Answer 1

[James Akula]

I belive the following functions perform the AIC test, as shown in https://www.mathworks.com/matlabcentral/answers/1893055-comparing-two-curve-fits-using-aic#comment_2561700. p1 is the probablility that the first model is the correct one.

% y   = y data
% o   = observed
% a   = "answer" to curve fit 
% m   = modeled
% 1/2 = model IDs
df     = @(oy, a) numel(oy) - numel(a);
SS     = @(my, oy) sum(((my - oy).^2));
DifAIC = @(m1y, m2y, oy, a1, a2) numel(oy) .* ...
    log(SS(m1y, oy)./SS(m2y, oy)) + 2 .* (df(oy, a2) - df(oy, a1));
p1 = @(m1y, m2y, oy, a1, a2) ...
    1 - exp(DifAIC(m1y, m2y, oy, a1, a2)/2)./...
    (1 + exp(DifAIC(m1y, m2y, oy, a1, a2)/2));

[James Akula]

I should have cited my source! Please see https://www.graphpad.com/guides/prism/latest/curve-fitting/reg_how_the_aicc_computations_work.htm for more details.

Answer 2

[Bora Eryilmaz]

The lsqcurvefit returns the residual norm of the fit. The model that has the least residual norm would be the "best" one by this criterion:

% Generate some data
b = 1; % True coefficient
n = 3; % True exponent
c = 0; % True y-intercept
rx = randi(100, [100, 1])/10;
rn = 3 + rand([100, 1])/4 - 0.125;
rc = (rand([100, 1]) * 500) - 250 + c;
ry = b * rx.^(rn) + rc;

% Generate and fit the models
m = @(b, c, n, x) b * x.^n + c;
m1 = @(p, x) m(p(1), p(2), 1, x); % n=1
m2 = @(p, x) m(p(1), p(2), 2, x); % n=2
m3 = @(p, x) m(p(1), p(2), 3, x); % n=3

[c1, resnorm1] = lsqcurvefit(m1, [1, 0], rx, ry);

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

[c2, resnorm2] = lsqcurvefit(m2, [1, 0], rx, ry);

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

[c3, resnorm3] = lsqcurvefit(m3, [1, 0], rx, ry);

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

% Plot the results
x = linspace(min(rx), max(rx));
plot(rx, ry, 'ok', x, m1(c1, x), '-g', x, m2(c2, x), '-r', x, m3(c3, x), '-b')

% Determine which model is most likely to be correct.  Use AIC?
[~, idx] = min([resnorm1 resnorm2 resnorm3]) % Index of best model

idx = 3

[Bjorn Gustavsson]

Sure but the Akaike Information Criterion (and its Bayesian IC and all their siblings) add a penalty cost to avoid(reduce the risk of/handle) the problem of overfitting, I think that's the feature the OP wonders about. (I'm by no means an expert on the model-selection problem, but comment since you're Mathworks stadd just to say: matlab has a couple of aic-functions but in some specialised add-on toolboxes and not in the statistics toolbox where it should have its "most basic home" - as it seems to me...)

[James Akula]

Thanks for the responses, both.

As noted above, it is trivial to find which curve fits the best by looking at straightforward statistics like the residual error. My real question, perhaps not well stated, is about confidence that one function is better than the other.

In the above instance, for example, let's say someone told me, "These data should be described the function y=b×x²+c." Now, I go and produce the red line in the plot above and, sure enough, it looks pretty good. However, after playing with the data, I also create the blue line and see that it fits just a bit better.

Now, I want to go back to that someone and say, "Actually, I think y=b×x³+c is better." I know this person is invested in n=2, and so I would like to be able to say how certain I am that we should change formulae to n=3.

[Bjorn Gustavsson]

(I'm getting out on dodgy looking thin ice over scary deep waters here) In your example case it should be as "simple" as to look at the parameter-covariance-matrix - since both models have the same number of parameters if you extend the model to also fit for the exponent. If the best value for the exponent is 3 with a standard deviation of "much less than 1" then this ought to show that it is better than an exponent of 2. When it comes to chosing between models with different numbers of parameters my "statistics vision" is not good enough to bring any enlightenment - I just go with the best AICc-model.

[Bora Eryilmaz]

I am not an expert on AIC criterion but as your models have the same number of parameters, it may not be very applicable in this scenario. In term of its model ranking power, AIC (with same number of parameters) would give similar results to looking at the residual norm. Assuming a normal distribution of modeling errors, log likelihood needed for AIC is basically the same as the residual norm.

In your scenario, you can rank the models looking at the residual norm of a "test dataset", and not the "training data" as was done above.

[James Akula]

Whether it is applicable or not, I cannot say for certain, but my understanding is that, indeed, it is. If I am correct, for functions with different numbers of degrees of freedom, you can use an F test, but AIC is "right" here (when the degrees of freedom are the same). I belive I have worked it out. Here is an example of code that seems to calculate the probabilities as expected.

%% Generate some data
rng(10)
b = 1; % True coefficient
n = 3; % True exponent
c = 0; % True y-intercept
rx = randi(100, [100, 1])/10;
rn = 3 + rand([100, 1])/4 - 0.125;
rc = (rand([100, 1]) * 500) - 250 + c;
rb = rand([100, 1]) + 0.5;
ry = rb .* rx.^(rn) + rc;

%% Generate and test the models
m  = @(b, c, n, x) b * x.^(n) + c;
m2 = @(p, x) m(p(1), p(2), 2, x);
m3 = @(p, x) m(p(1), p(2), 3, x);

f2 = @(p) rmse(m2(p, rx), ry);
f3 = @(p) rmse(m3(p, rx), ry);

a2 = fminsearch(f2, [1, 0])';
a3 = fminsearch(f3, [1, 0])';

%% Plot the results
x = linspace(min(rx), max(rx));
plot(rx, ry, 'ok', x, m2(a2, x), '-r', x, m3(a3, x), '-b')

%% AIC test
% o = observed data
% a = "answer" to curve fit 
% m = modeled data
df     = @(oy, a) numel(oy) - numel(a);
SS     = @(my, oy) sum(((my - oy).^2));
DifAIC = @(m1y, m2y, oy, a1, a2) numel(oy) .* ...
    log(SS(m1y, oy)./SS(m2y, oy)) + 2 .* (df(oy, a2) - df(oy, a1));

p1 = @(m1y, m2y, oy, a1, a2) ...
    1 - exp(DifAIC(m1y, m2y, oy, a1, a2)/2)./...
    (1 + exp(DifAIC(m1y, m2y, oy, a1, a2)/2));

pm2 = p1(m2(a2, rx), m3(a3, rx), ry, a2, a3);
pm3 = p1(m3(a3, rx), m2(a2, rx), ry, a3, a2);

disp(array2table(round([[a2, a3];[pm2, pm3] * 100], 4, 'significant'), ...
    'VariableNames', {'n = 2', 'n = 3'}, ...
    'RowNames', {'b', 'c', '% is Correct'}));

                    n = 2     n = 3 
                    ______    ______

    b                8.192    0.8739
    c               -34.61     23.43
    % is Correct     20.98     79.02

[Bjorn Gustavsson]

It really makes me queasy to say that one model is true with some probability here, when it is so simple to use the exponent as another free parameter and then calculate the parameter covariance-matrix to get an estimate of the uncertainty of the exponent from its standard deviation. I know this is a simplified toy-model example, but it really makes me uneasy...

[James Akula]

I hear ya. But I do feel like if you have a priori reasons to think something is either one thing or another, and you are hopitn the data will guide you to the truth, then this approach perhaps makes sense.

[Bjorn Gustavsson]

@James Akula, sure, and my suggestion might work well enough for this (with reasonable certainty) simplified toy-example, but can run into severe problems in more complicated real-world cases. I mainly thought it was worth mentioning as a reminder for the case where someone were to have this specific choise to make. (and I don't need anyone to argue this quandry with - I come down on all three of the two sides in this discussion and have frequent and angry arguments with myself from last week...)