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

 
 

bug #38063: [octave forge] (optim) expfit name collision with statistics pkg, possible accuracy issues

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

Add a New Comment Rich Markup
   

Jump to the original submission

Sat 01 Jan 2022 06:33:52 PM UTC, comment #15: 

Renamed optims expfit to pronyfit, closing this bug.

Olaf Till <i7tiol>
Group Member
Tue 16 Nov 2021 01:19:22 AM UTC, comment #14: 

retitling this bug to indicate that the current state of the problem includes the function needing to be renamed (for matlab compatibility and naming issues with the function now pushed to the statistics package), as well as general performance concerns.

I believe the current state needs:
(1) short term - a function rename for the optim version of the function, and/or (2) a complete rework to implement the matlab 'fit' function that maybe could borrow some of the code from optim's expfit.

Nicholas Jankowski <nrjank>
Group Member
Sat 07 Aug 2021 08:23:50 AM UTC, comment #13: 

note patch #10092 provides a matlab compatible expfit for the statistics package

Nicholas Jankowski <nrjank>
Group Member
Thu 05 Aug 2021 03:12:22 PM UTC, comment #12: 

For Matlab compatibility, it would be a good idea to rename our expfit, to for example fit_expn

Nir Krakauer <nir_krakauer>
Thu 05 Aug 2021 03:04:32 PM UTC, comment #11: 

that's a good point about the conflicting "unimplemented" listing.

Nicholas Jankowski <nrjank>
Group Member
Thu 05 Aug 2021 08:02:02 AM UTC, comment #10: 

because expfit in Matlab is about the exponential distribution, and thus totally different, I suggest to at least flag this as a Matlab incompatibility bug. But I generally favor renaming.

Note also that 'help expfit' tells us

error: help: the 'expfit' function belongs to the statistics package from Octave
Forge but has not yet been implemented.

which indicates that it should be implemented as a distribution fit as in ML, not as an optim-function.

Muhali <muhali>
Wed 04 Aug 2021 08:02:59 PM UTC, comment #9: 

this bug is still sitting open and the conditions haven't changed judging by the lack of anything significant in the changelog for expfit.  some compatibility namespace issues with Matlab were identified in comment #7 and 8. 

octave has no prony implementation, which would fit into octave-signal as it is a filter design function based on Prony fitting. (and is listed as a missing function: https://wiki.octave.org/Signal_package)
matlab doc: https://www.mathworks.com/help/signal/ref/prony.html

matlab's expfit function is part of the Statistics and Machine Learning toolbox, and is currently missing from Octave's statistics package, that is supposed to calculate the mean value of a set of data assumed to be taken from an exponential distribution and optionally provide the confidence intervals. (https://wiki.octave.org/Statistics_package)
matlab doc: https://www.mathworks.com/help/stats/expfit.html

expfit in matlab seems to be a fairly simple function to implement, (albeit sitting in a long line of unimplemented statistics functions.) but if implemented would produce a namespace collision.

the function expfit seems to best match the matlab function fit, which has 1st and 2nd order exponentials as available models exp1 and exp2, expfit being something that could be an 'expn' model for it.  however, octave currently doesn't have fit implemented either.

I'd recommend closing this bug, which seems to mainly be about the performance of expfit but hasn't been changed in 8 years, and/or retitling or making a new bug about the incompatible namespace issue.

Nicholas Jankowski <nrjank>
Group Member
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:

http://www.mathworks.com/help/signal/ref/prony.html

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>
Group Member
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>
Group 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 =

  -0.093230
  -0.119135
   2.146504
   2.146838

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 =

   2.0000
   1.3000

c =

  -0.50002
   1.99994

rms =  0.0018253
octave:21>


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>
Group 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:

https://mailman.cae.wisc.edu/pipermail/octave-maintainers/2013-April/032944.html

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>
Group 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:
   
   
Comment:
   

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 nir_krakauer (Posted a comment)
  • -email is unavailable- added by nrjank (Posted a comment)
  • -email is unavailable- added by jwe (Updated the item)
  • -email is unavailable- added by i7tiol (Posted a comment)
  • -email is unavailable- added by muhali (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 6 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2022-01-03 mtmiller Carbon-CopyRemoved 80942 -
    2022-01-01 i7tiol Open/ClosedOpen Closed
    2021-11-16 nrjank Summary[octave forge] (optim) expfit inaccurate (incorrect) [octave forge] (optim) expfit name collision with statistics pkg, possible accuracy issues
    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.13-f8d8.
    Corresponding source code