bugGNU Octave - Bugs: bug #38063, [octave forge] (optim) expfit...


bug #38063: [octave forge] (optim) expfit inaccurate (incorrect)

Submitted by:  Muhali <muhali>
Submitted on:  Fri 11 Jan 2013 02:15:26 PM UTC  
Category:  Octave Forge Package Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Inaccurate Result
Status:  Need Info Assigned to:  None
Originator Name:  Open/Closed:  Open
Release:  dev Operating System:  Any

Add a New Comment (Rich Markup)

You are not logged in

Please log in, so followups can be emailed to you.


( Jump to the original submission)

Wed 10 Apr 2013 06:27:19 AM UTC, comment #8: 

>> Can we get someone to run some expfit test cases in Matlab and compare results?

expfit in ML is about the 1-parameter exponential distribution, a*exp(a*x), unlike expfit in octave.

prony in ML is about impulse response functions, and I have the impression that is what the current expfit actually does.

Muhali <muhali>
Wed 10 Apr 2013 01:05:45 AM UTC, comment #7: 

The name prony is taken, there is an existing Matlab function with that name:


Maybe bits of this function can be used to implement that, but that's a separate issue.

Can we get someone to run some expfit test cases in Matlab and compare results?

Mike Miller <mtmiller>
Project Administrator
Tue 09 Apr 2013 08:27:19 PM UTC, comment #6: 

This leaves the most frequently encountered, hence most important case of a simple exponential fit, with expfit utterly failing.

So if I were a natural or life scientist, reading in the expfit help that 'This method is known to be more stable than 'brute-force' non-linear least squares fitting', I would probably be somewhat clueless about what octave is offering.

To avoid this, what is wrong with renaming 'expfit' to 'prony' (plus correcting the help text) and have 'expfit' reserved for the case above and use a simple nonlinear regression, or just not have any 'expfit' at all?

Muhali <muhali>
Tue 09 Apr 2013 05:20:52 PM UTC, comment #5: 

Ok, I confirmed that nonlin_curvefit finds the parameters with suitable initial values.

But now the situation, as I see it, is:

- little noise:

-- expfit finds the correct parameters

-- nonlin_residmin finds them, too, but only with suitable initial values (since you don't know the correct parameters beforehand, this may be a problem)

- much noise:

-- expfit does not find the correct parameters and not even a good fit

-- nonlin_residmin also does not find the correct parameters, but achieves a good fit

So with little noise having expfit may be an advantage. With much (realistic?) noise possibly nothing helps; except if you only want to fit, not to find the parameters -- then you have to take a general optimizer.

In other words, if you need the correct parameters, not only the fit, you are either at the loss due to much (realistic?) noise, or you should better use expfit.

I still don't see how renaming expfit helps to clarify this situation.

It is a matter of ones emphasis how bad it is that expfit doesn't achieve a good fit if it can't find the correct parameters anyway; maybe this is even a feature of the algorithm (which I havn't checked). In life- (and probably other) science, one is often concerned with finding the correct parameters of a hopefully suitable model; if this isn't successful and a good fit is still needed, one can still switch to a general optimizer. If you feel that expfit should still give a good fit in this case, the right way is probably to check whether the algorithm is principally able to provide this and make a respective correction.

Olaf Till <i7tiol>
Project Member
Tue 09 Apr 2013 04:22:19 PM UTC, comment #4: 

You are giving [-1; -1; 1; 1] as an initial guess for nonlin_curvefit, that is, two negative exponentials, so no wonder the iteration gets stuck. That the fit still looks ok is because -0.5exp(2*x) dominates for larger x. Using MaxIter = 1000 and [0; 0; 1; 1] as initial guess makes it work.

I would still prefer renaming expfit, because who wants to have a function of that name that gives largely incorrect estimates in the most simple case of y = a*exp(b*x) with realistic noise?

Muhali <muhali>
Tue 09 Apr 2013 01:26:58 PM UTC, comment #3: 

In my tests, MaxIter was not limiting (it was set to 100 and there were 34 and 32 iterations) and yet general optimization, though giving a good fit, did not come up with reasonably correct parameters, while expfit did the latter for the small noise case. (I'm pasting only the small-noise case.)

octave:1> x0 = 0; step = 0.05; xend = 5; x = x0:step:xend;
octave:2> y = 2*exp(1.3*x)-0.5*exp(2*x);
octave:4> serror = (rand(1,length(y))-0.5)*3e-4;
octave:16> [p, fy, cvg, outp] = nonlin_curvefit (@ (p, x) p(1) * exp (p(3) * x) + p(2) * exp (p(4) * x), [-1; -1; 1; 1], x, y + serror, optimset ("maxiter", 100));
octave:17> p
p =


octave:18> cvg
cvg =  3
octave:19> outp
outp =

  scalar structure containing the fields:

    niter =  34

octave:20> [alpha,c,rms] = expfit(2,x0,step,y+serror)
alpha =


c =


rms =  0.0018253

I still suspect a misunderstanding. Your emphasis seems to be on curve fit, while possibly in the case of the models dealt with by expfit, parameter determination could be more important. My argument was that in the case of the example with the large noise, although you say this is more like real world, there is possibly no optimizer at all which can figure out the correct model parameters ('correct' in the sense of 'underlying' (assuming the model reflects the real conditions producing the data), not in the sense of producing a good fit).

Olaf Till <i7tiol>
Project Member
Tue 09 Apr 2013 10:10:19 AM UTC, comment #2: 

>> - The example in the report fits a degree-1 model to data produced with a degree-2 mode

This doesn't matter. You get similar (bad) expfit results with degree-1 data, e.g. y=-0.5*exp(2*x). Moreover, with deg=2 expfit produces non-real parameters.

>> - Although general nonlinear fitting gave a good fit, the determined model parameters (for a degree-2 model) were vastly incorrect

That is simply because the default MaxIter parameter in nonlin_curvefit is set to 20; increasing that to 1000 gives a perfect fit (in reasonable time).

Again, if there is no noise there is not much to fit, and real-world cases look more like the figure attached to the original post, for which expfit gets it totally wrong.

I suggest to at least add a warning that expfit should not be used for classical fitting problems. How about renaming it to e.g. prony_fit.m?

Muhali <muhali>
Tue 09 Apr 2013 08:28:16 AM UTC, comment #1: 

Not knowing this algorithm, I asked at the maintainers list:


Summary of the thread and some tests:

- The example in the report fits a degree-1 model to data produced with a degree-2 model. But fitting a degree-2 model still gives bad results with this data for expfit.

- Although general nonlinear fitting gave a good fit, the determined model parameters (for a degree-2 model) were vastly incorrect, even for the not so noisy data in the help-text-example where expfit is able to determine the correct parameters.

- Looking at the noise in the data in the report and the used model I'd feel that probably no optimizer is able to figure out reasonably correct parameters.

So, since expfit, but not easily the general optimizer, was able to determine correct parameters in the example with little noise, and this is probably just not possible with the example with much noise, I'd see no need to deprecate or change expfit at the moment. This assumes that the emphasis with expfit is on parameter determination, not just curve fitting.

Olaf Till <i7tiol>
Project Member
Fri 11 Jan 2013 02:15:26 PM UTC, original submission:  

I get very bad results from the expfit estimator for exp fitting. The function has a demo that uses extremely low levels of noise. With more realistic levels the fit becomes useless. In this example

pkg load optim

## expfit
x0 = 0; step = 0.05; xend = 5; x = x0:step:xend;
y = 2*exp(1.3*x)-0.5*exp(2*x);
error = (rand(1,length(y))-0.5)*3e3;
[alpha,c,rms] = expfit(1,x0,step,y+error);

## nonlin_curvefit
fun = @(p,x) p(2)*exp(p(1)*x) ; p0 = [1;-1] ;
[p, fy, cvg, outp] = nonlin_curvefit(fun, p0, x, y+error)

hold on
plot(x, y=y+error, "k.")
plot(x, c*exp(alpha*x), "r-;expfit;")
plot(x, fy, "b-;nonlin_curvefit;")
legend("location", "southwest")
print foo.png

expfit is compared with a more conventional iterative technique, showing the discrepancies (see attached figure; the demo noise level had to be expanded by 7(!) orders of magnitude).

Apparently expfit is based on Prony's method, a DSP technique similar to Fourier analysis. Either that method cannot be used for fitting noisy data, or its implementation in expfit is incorrect. In either case I consider the expfit function buggy, and recommend not using it.

Muhali <muhali>


(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:

Attached Files
file #27243:  foo.png added by muhali (46KiB - image/png)


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -email is unavailable- added by jwe (Updated the item)
  • -email is unavailable- added by mtmiller (Posted a comment)
  • -email is unavailable- added by i7tiol (Posted a comment)
  • -email is unavailable- added by muhali (Submitted the item)

    Do you think this task is very important?
    If so, you can add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.




    Follow 3 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2017-08-13 jwe Summaryoptim: expfit inaccurate (incorrect) => [octave forge] (optim) expfit inaccurate (incorrect)
    2013-04-10 mtmiller StatusNone => Need Info
    2013-01-11 muhali Attached File- => Added foo.png, #27243

    Back to the top

    Powered by Savane 3.6