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

 
 

bug #54009: [octave forge] (statistics) incorrect result for gppdf

Submitter:  Muhali <muhali>
Submitted:  Wed 30 May 2018 11:20:57 AM UTC
   
 
Category:  Octave Package Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  lostbard
Originator Name:  Open/Closed:  * Closed
Release:  * 4.4.0 Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Mon 30 Jul 2018 02:18:03 AM UTC, comment #8: 

Closing as fixed - changes will appear in next release of statistics

John Donoghue <lostbard>
Group Member
Wed 18 Jul 2018 03:45:36 PM UTC, comment #7: 
John Donoghue <lostbard>
Group Member
Tue 17 Jul 2018 05:11:06 PM UTC, comment #6: 

All the other gpXXXX functions are also wrong.

In addition, at a quick look, bbsXXXX, cauchy_XXXXX functions (although not available as is in matlab) are the same format and would make sense to also be swapped.

John Donoghue <lostbard>
Group Member
Tue 17 Jul 2018 07:27:02 AM UTC, comment #5: 

works for me now, thanks. But recall that (LOCATION, SCALE, SHAPE) is still in wrong order (compared to gev, or Matlab).

Muhali <muhali>
Thu 12 Jul 2018 07:15:53 PM UTC, comment #4: 
John Donoghue <lostbard>
Group Member
Thu 12 Jul 2018 07:10:50 PM UTC, comment #3: 

The last 2 changed lines in your diff should be ./ scale(j)

John Donoghue <lostbard>
Group Member
Thu 12 Jul 2018 06:58:46 PM UTC, comment #2: 

Applying that change in the code causes 14 of the tests to fail.

Not sure if that's because the tests are not valid or something else.



For example:

x = [-Inf, -1, 0, 1/2, 1, Inf]

gppdf(x, zeros(1,6), ones(1,6), ones(1,6))



John Donoghue <lostbard>
Group Member
Wed 27 Jun 2018 08:37:59 AM UTC, comment #1: 

here is more evidence, if needed:

octave> Q = integral(@normpdf, -Inf, Inf)
Q =  1.000
octave> Q = integral(@(x) gampdf(x, 1, 1), -Inf, Inf)
Q =  1.00
octave> Q = integral(@(x) gevpdf(x, 1, 2, 0), -Inf, Inf)
Q =  1.000
octave> Q = integral(@(x) gppdf(x, 0, 2, 1), -Inf, Inf)
Q =  2.00
octave> Q = integral(@(x) gppdf(x, 0, 3, 1), -Inf, Inf)
Q =  3.00
octave> Q = integral(@(x) gppdf(x, 0, 3, 2), -Inf, Inf)
Q =  3.00


Please note also that the order of arguments (LOCATION, SCALE, SHAPE) is reversed relative to Matlab.

I attach a patch that fixes the scale issue.


--- /usr/local/src/octave-statistics/inst/gppdf.m
+++ /tmp/gppdf.m
@@ -69,13 +69,13 @@

     j = k & (shape > 0) & (z >= 0);
     if (any (j))
-      pdf(j) = (shape * z(j) + 1).^(-(shape + 1) / shape);
+      pdf(j) = (shape * z(j) + 1).^(-(shape + 1) / shape) ./ scale;
     endif

     if (shape < 0)
       j = k & (shape < 0) & (0 <= z) & (z <= -1. / shape);
       if (any (j))
-        pdf(j) = (shape * z(j) + 1).^(-(shape + 1) / shape);
+        pdf(j) = (shape * z(j) + 1).^(-(shape + 1) / shape) ./ scale;
       endif
     endif
   else
@@ -88,13 +88,13 @@

     j = k & (shape > 0) & (z >= 0);
     if (any (j))
-      pdf(j) = (shape(j) .* z(j) + 1).^(-(shape(j) + 1) ./ shape(j));
+      pdf(j) = (shape(j) .* z(j) + 1).^(-(shape(j) + 1) ./ shape(j)) ./ scale;
     endif

     if (any (shape < 0))
       j = k & (shape < 0) & (0 <= z) & (z <= -1 ./ shape);
       if (any (j))
-        pdf(j) = (shape(j) .* z(j) + 1).^(-(shape(j) + 1) ./ shape(j));
+        pdf(j) = (shape(j) .* z(j) + 1).^(-(shape(j) + 1) ./ shape(j)) ./ scale;
       endif
     endif
   endif


Muhali <muhali>
Wed 30 May 2018 11:20:57 AM UTC, original submission:  

The Generalized Pareto function gppdf produces incorrect results. It misses a factor 1/sigma, which is evident if one compares the relevant line of code

pdf(j) = (shape(j) .* z(j) + 1).^(-(shape(j) + 1) ./ shape(j));

with (e.g.) the corresponding Wikipedia entry.

Muhali <muhali>

 

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

Attach Files:
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by lostbard (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
    2018-07-30 lostbard StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2018-07-12 lostbard StatusNone Ready For Test
        Assigned toNone lostbard
    2018-06-04 mtmiller Carbon-CopyRemoved 80942 -
    2018-06-04 mtmiller Summary(statistics) incorrect result for gppdf [octave forge] (statistics) incorrect result for gppdf

    Back to the top

    Powered by Savane 3.13-758e.
    Corresponding source code