bugGNU Octave - Bugs: bug #47369, [octave forge] (statistics) gevfit...

 
 

bug #47369: [octave forge] (statistics) gevfit malfunctioning (because of fminunc failure)

Submitter:  Muhali <muhali>
Submitted:  Wed 09 Mar 2016 10:20:23 AM UTC
   
 
Category:  Octave Package Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * dev Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Wed 28 Sep 2022 12:59:44 AM UTC, comment #12: 

This has been fixed with commit 3ac6fcb pushed at https://github.com/gnu-octave/statistics/commit/3ac6fcb69c7b32589dcb184407aae14f60a05058
Replaced fminunc with fminsearch for minimization. Added BISTs, clarified help, made it Matlab compatible. This can be closed.

Andreas Bertsatos <pr0m1th3as>
Fri 12 Aug 2022 01:45:12 PM UTC, comment #11: 

Would is help to introduce `nelder_mead' as a private function in statistics in order to solve this issue?

Andreas Bertsatos <pr0m1th3as>
Mon 08 Aug 2022 06:57:33 PM UTC, comment #10: 


comment #9:

> your first example doesn't have a function f to feed into nelder_mead_min.


On line 63 there is

f = @(p) gevlike(p, data);

and nelder_mead seems to be working (= outputs the same as ML), but fminunc is not (at least in this example).

Muhali <muhali>
Mon 08 Aug 2022 05:29:30 PM UTC, comment #9: 

your first example doesn't have a function f to feed into nelder_mead_min.

I'm mainly asking if the issue is with gevfit in the statistics package (which seems to only handle initial param guessing then call other functions), gevlike (cost function passed to fminunc by gevfit, also in the statistics package), or something specific to fminunc (not in the statistics package). 

if the latter, then the most useful course wolud be to
(1) retitle this report removing reference to the statistics package.
(2) come up with a minimal fminunc only test case that highlights the specific problem
(3) come up with a more descriptive title for the bug report.

Nicholas Jankowski <nrjank>
Group Member
Mon 08 Aug 2022 04:08:03 PM UTC, comment #8: 

this is definitely a fminunc issue. If you replace (in the 1st example) the fminunc line with

paramhat = nelder_mead_min(f, paramguess);

you get the ML result.

Muhali <muhali>
Sat 06 Aug 2022 02:26:12 PM UTC, comment #7: 

looking at the status of this report:

with octave 7.2, statistics 1.4.3,
that the comment #0 supplied tests in both matlab and octave produce:


x = [549,664,382,540,856,1190,432,542,383,251,967,876,1190,302,740,981,620,423,753,801,1100,457,763,572,526,430,459,597,1310,612,896]';

[PARAMHAT, PARAMCI] = gevfit (x)

Octave:
PARAMHAT =
  -6.9804e-03   2.1645e+02   5.5777e+02

PARAMCI =
    -0.3402   151.3145   467.3575
     0.3262   281.5900   648.1833

Matlab 2022a:

PARAMHAT =
   -0.0092  219.5077  556.1071
PARAMCI =
   -0.3542  160.9012  465.2270
    0.3359  299.4610  646.9872


testing the Octave function example in the gevfit help:


data = 1:50;
[pfit, pci] = gevfit (data)

Octave:
pfit =
   -0.4407   15.1923   21.5309

pci =
   -0.7532   11.0777   16.5688
   -0.1282   19.3070   26.4930


Matlab 2022a:
pfit =
   -0.4407   15.1924   21.5306
pci =
   -0.7532   11.5878   16.5686
   -0.1282   19.9183   26.4926



there is discrepancy between octave and matlab in both examples, although in the second example it is much less.

is this a gevfit issue or a fminunc issue? algorithm or convergence detection? perhaps someone more familiar with the underlying deficiency can improve the title/summary?

Nicholas Jankowski <nrjank>
Group Member
Wed 24 Nov 2021 07:49:55 PM UTC, comment #6: 

Your comparison with Matlab is misleading, because Matlab's gevfit has no equivalent of the paramguess input (so the [0 1 0] is not treated as the initial guess). I think it is useful for now to keep paramguess in our version to allow users more control. But with an equivalent input configuration as Matlab (i.e. no paramguess), our version also converges.

In this particular case and for an initial guess of [0 1 0], nelder_mead_min would do a better job of finding the optimum than fminunc. However, it isn't clear whether this applies generally.

I now added a warning based on the status flag returned by fminunc. However, for the case gevfit(x, [0 1 0]), the flag indicates convergence, which may be a bug in fminunc, so there's no warning if gevfit is called with 1 output. If gevfit is called with 2 outputs, there is a warning, because then the Fisher information matrix is computed (and found to be not positive definite).

I also changed the returned vector orientation to be the same as Matlab.

https://sourceforge.net/p/octave/statistics/ci/b11f43821a99b693e9e92218c83f304c4d46a469/

Nir Krakauer <nir_krakauer>
Wed 24 Nov 2021 06:10:24 PM UTC, comment #5: 

this is still present in the Octave 6.4.0, statistics 1.4.2.

the comment 0 test:

matlab 2021a:

 p = gevfit(x, [0 1 0])

 p =
   -0.0092  219.5077  556.1071


octave:

p = gevfit(x, [0 1 0])
p =

     7.3956    35.1377   255.7373


using the suggested optim package function:

p =

  -9.1576e-03   2.1951e+02   5.5611e+02



No warning is produced as suggested by Nir in comment #3, so not sure which repository change he was referring to and if it made it into the released version.

checking the 'automatic param estimation', I notice the default is a column not a row vector, and so get:

Octave:

>> p = gevfit(x)
p =

  -6.9804e-03
   2.1645e+02
   5.5777e+02


Matlab:

>> p = gevfit(x)
p =
   -0.0092  219.5077  556.1071


so in addition to the output discrepancy this is also a format compatibility concern.  that latter part might be worth a quick correction separate from solving the former problem.

Nicholas Jankowski <nrjank>
Group Member
Wed 09 Mar 2016 08:24:57 PM UTC, comment #4: 

Thanks to both, and thanks for the fix. I agree with you that one should stick to the octave core as much as possible. But if an optimization fails so gravely (see attached cdf.png) without issuing any warning, then the 'normal' user will be turned away immediately, especially since EV analysis is quite relevant.

And as mentioned, it seems (the core function) fminunc has trouble with gevlike, without throwing errors. I suspect the gradient calculation fails near the domain boundary.


Muhali <muhali>
Wed 09 Mar 2016 08:02:57 PM UTC, comment #3: 

Hi Muhali-- I agree with Arno that it's preferable to stick to functions in core Octave as much as possible so that Statistics doesn't require other Forge packages. It is to be expected that any optimization algorithm used in gevfit could fail to converge in some cases, because the objective function is quite nonlinear. Whether the optimization converges to a likelihood maximum or not for particular data typically depends on the initial value, which can be user-supplied or determined automatically.

In response to your problem, I improved the automatically estimated initial value so that gevfit (in the repository) now converges properly for your example with automatic choice of initial value (p = gevfit(x);). p = gevfit(x, [0 1 0]); will still not converge, but now there is a warning that the optimization didn't work.

Best, Nir

Nir Krakauer <nir_krakauer>
Wed 09 Mar 2016 06:34:50 PM UTC, comment #2: 

Thanks for the bug report. However, I am reluctant to replace fminunc with nelder_mead_min (and add the optim package dependency). Which of the two optimization algorithms performs better strongly depends on the data. In this case the simplex method performed better, in other cases the gradient based method will perform better.

Arno Onken <asnelt>
Wed 09 Mar 2016 05:54:25 PM UTC, comment #1: 

Adding maintainer to CC.

Carnë Draug <carandraug>
Group Member
Wed 09 Mar 2016 10:20:23 AM UTC, original submission:  

With the present version of the statistics package (1.2.4), the gevfit function returns unusable results.

The reason is that fminunc fails to find an appropriate minimum, with nelder_mead_min, for example, obviously working. I actually checked the fminunc solution, and indeed the objective function decreased in the immediate eps neighborhood. So it seems that this is in fact a fminunc bug.

The failure is demonstrated by the following case:


pkg load statistics optim
load foo.ot
source gevfit.m
p = gevfit(x, [0 1 0]);


The corresponding files are attached (with 3 lines added to gevfit that show the sub-optimal solution as compared to nelder_mead).

Muhali <muhali>

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #36587:  cdf.png added by muhali (18KiB - image/png)
file #36575:  x.ot added by muhali (208B - application/octet-stream)
file #36576:  gevfit.m added by muhali (3KiB - text/x-objcsrc)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by pr0m1th3as (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 nir_krakauer (Posted a comment)
  • -email is unavailable- added by asnelt (Posted a comment)
  • -email is unavailable- added by carandraug (Arno Onken - maintainer of statistics package)
  • -email is unavailable- added by muhali (Submitted the item)
  • -email is unavailable- added by muhali
  •  

    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 9 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2022-09-28 nrjank StatusNeed Info Fixed
        Open/ClosedOpen Closed
    2022-08-06 nrjank StatusNone Need Info
    2017-08-13 jwe Summarystatistics: gevfit malfunctioning (because of fminunc failure) [octave forge] (statistics) gevfit malfunctioning (because of fminunc failure)
    2016-03-09 muhali Attached File- Added cdf.png, #36587
    2016-03-09 carandraug Carbon-Copy- Added -email is unavailable-
    2016-03-09 muhali Attached File- Added x.ot, #36575
        Attached File- Added gevfit.m, #36576
    2016-03-09 muhali Carbon-Copy- Added asnelt

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code