bugGNU Octave - Bugs: bug #36718, moment, skewness, and kurtosis...

 
 

bug #36718: moment, skewness, and kurtosis differences with Matlab

Submitter:  Alexander Barth <abarth>
Submitted:  Sun 24 Jun 2012 08:46:06 PM UTC
   
 
Category:  Interpreter Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Matlab Compatibility
Status:  Fixed Assigned to:  None
Originator Name:  Alexander Barth Open/Closed:  * Closed
Release:  * 3.6.0 Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Thu 07 Nov 2013 05:41:05 PM UTC, comment #30: 

I checked in a final changeset to have the moment function default to central moments instead of raw moments (http://hg.savannah.gnu.org/hgweb/octave/rev/58b39152b0f6).  All three changes to moment, skewness, and kurtosis will be part of version 3.8.0.  Closing report.

Rik <rik5>
Group administrator
Sat 19 Oct 2013 02:39:09 PM UTC, comment #29: 

I checked in your changeset here (http://hg.savannah.gnu.org/hgweb/octave/rev/f0e777cf348f).  I also announced in the NEWS file that the definition of the kurtosis function has changed from "excess kurtosis" to just kurtosis.  The NEWS file is a ChangeLog for user-visible changes that are important or that might change behavior from previous versions.

Rik <rik5>
Group administrator
Sat 19 Oct 2013 10:13:20 AM UTC, comment #28: 

Here is a new version of the patch :

  • the tests have been modified as discussed,
  • the documentation has been updated.




(file #29418)

Julien Bect <jbect>
Sat 19 Oct 2013 06:27:12 AM UTC, comment #27: 

For compatibility with MAtlab's kurtosis function, I have changed the definition kurtosis (x).

The new value of kurtosis (x) is equal to the former one ("excess kurtosis") + 3.

Should I mention that in the documentation ? or somewhere else ?

Julien Bect <jbect>
Sat 19 Oct 2013 06:14:01 AM UTC, comment #26: 

Ok, I'll raise te tolerance.

Julien Bect <jbect>
Fri 18 Oct 2013 08:23:55 PM UTC, comment #25: 

That is strange, but I'm inclined to accept the results of the arbitrary precision arithmetic.  Raising the tolerance is fine.

Rik <rik5>
Group administrator
Fri 18 Oct 2013 08:00:59 PM UTC, comment #24: 

This is weird, because I get

5.437731792528892

on my computer with both R 3.0.0, Matlab R2012a and my latest Octave build. On the other hand, using mpmath with 30 digits of precision I get the following interval:

mpi('5.43773179252889008331093791989766', '5.43773179252889008331093791992922')

which means that the value that you report is the correct one. I will change the current reference values for more precise ones:

2.9786509002956195 (sample kurtosis)
5.4377317925288901 (with bias corrections)

but then I need to raise the tolerance to 8 * eps on my computer to get the tests to pass. Is that OK ?

Here is the Python code that I use:


from mpmath import iv
iv.dps = 30

x = iv.matrix ([1, 2, 3, 4, 5, 10])
n = iv.mpf('6')
x = x - sum (x) / n

x2 = iv.matrix ([0, 0, 0, 0, 0, 0])
for i in range(0, 6):
    x2[i] = x[i] * x[i];

x4 = iv.matrix ([0, 0, 0, 0, 0, 0])
for i in range(0, 6):
    x4[i] = x2[i] * x2[i];

mu4 = sum (x4) / n
mu2 = sum (x2) / n

k1 = mu4 / (mu2 * mu2)

C = (n - 1) / ((n - 2) * (n - 3))
k0 = 3 + C * ((n + 1) * k1 - 3 * (n - 1))



Julien Bect <jbect>
Fri 18 Oct 2013 05:41:02 PM UTC, comment #23: 

I ran 'test kurtosis' and I get a failing test


 ***** assert (kurtosis ([1:5 10; 1:5 10],  0, 2), 5.437731792528892 * [1; 1], eps)
!!!!! test failed
ASSERT errors for:  assert (kurtosis ([1:5, 10; 1:5, 10], 0, 2),5.437731792528892 * [1; 1],eps)

  Location  |  Observed  |  Expected  |  Reason
    (1)         5.4377       5.4377      Abs err 1.7764e-15 exceeds tol 2.2204e-16
    (2)         5.4377       5.4377      Abs err 1.7764e-15 exceeds tol 2.2204e-16


How did you determine the expected value for this test?  Did you use an arbitrary precision calculator?  The results for Octave show a difference of '2' in the last significant digit.


format long
kurtosis ([1:5 10; 1:5 10],  0, 2)
ans =
   5.437731792528890
   5.437731792528890



Rik <rik5>
Group administrator
Fri 18 Oct 2013 05:20:27 PM UTC, comment #22: 

Here is a new version of the patch with the two modifications that you suggest. I will fix the documentation during the week-end.

(file #29407)

Julien Bect <jbect>
Fri 18 Oct 2013 03:20:14 PM UTC, comment #21: 

The patch looks good.  I might change


y(idx) = 3 + C * ((n + 1) * y(idx) - 3 * (n - 1));


to


y = 3 + C * ((n + 1) * y - 3 * (n - 1));


Since it "looks" more like a straightforward equation.  You don't have to worry about the values not in idx because they have been assigned the value of NaN, and any further math operation with a NaN value produces a NaN.

Also, this was my mistake in one of the tests for division-by-zero.  I save and restore the warning state but don't actually force the "divide-by-zero" state to on.  I've corrected that in skewness so that the test reads


## Verify no "divide-by-zero" warnings
%!test
%! wstate = warning ("query", "Octave:divide-by-zero");
%! warning ("on", "Octave:divide-by-zero");
%! unwind_protect
%!   lastwarn ("");  # clear last warning
%!   skewness (1);
%!   assert (lastwarn (), "");
%! unwind_protect_cleanup
%!   warning (wstate, "Octave:divide-by-zero");
%! end_unwind_protect


Rik <rik5>
Group administrator
Fri 18 Oct 2013 11:29:43 AM UTC, comment #20: 

Here is a first version of my patch for kurtosis.

It's not ready to be pushed yet, the documentation needs to be rewritten, but the code and the tests are ready for review.

@++


(file #29401)

Julien Bect <jbect>
Mon 14 Oct 2013 07:02:25 AM UTC, comment #19: 

Hello Rik,

I will do my best to provide a similar patch for kurtosis before Friday.

@++
Julien


Julien Bect <jbect>
Mon 14 Oct 2013 01:23:53 AM UTC, comment #18: 

Julien,

We are getting ready for the 3.8 release.  If you have any further patches for kurtosis or moment can I get them within the next week?

Thanks,
Rik

Rik <rik5>
Group administrator
Mon 07 Oct 2013 08:14:02 PM UTC, comment #17: 

I added the fix to the comment here (http://hg.savannah.gnu.org/hgweb/octave/rev/a53a97faf61a).

Rik <rik5>
Group administrator
Mon 07 Oct 2013 07:12:19 PM UTC, comment #16: 

Oooops, here is the patch, sorry about that.

(file #29309)

Julien Bect <jbect>
Mon 07 Oct 2013 07:11:45 PM UTC, comment #15: 

Another tiny patch on skewness.m (Fix comment).

Julien Bect <jbect>
Mon 07 Oct 2013 02:52:31 PM UTC, comment #14: 

Julien,

I added your changes here (http://hg.savannah.gnu.org/hgweb/octave/rev/8334144458a4).

I thought there was a potential for a "division-by-zero" error which is why I had coded skewness the way it was a long time ago.  I added a %!test to verify that if the code is ever changed again we don't re-introduce the same error.

Rik <rik5>
Group administrator
Mon 07 Oct 2013 01:28:00 PM UTC, comment #13: 

Rik,

Please find attached three patches for skewness.m :

a) Your coding style fixes make some test fail on my computer (that's because you replaced some explicite tolerance values by eps). I simply replace eps by 2 * eps.

b) Because of some changes you made in input argument checking, weird (non scalar) values of flag were able to go though undetected. I fixed that in a manner, I hope, you will find in agreement with Octave core coding conventions.

c) Get rid of a "division by zero" warning that I hadn't noticed before.

@++
Julien


(file #29306, file #29307, file #29308)

Julien Bect <jbect>
Mon 07 Oct 2013 02:58:03 AM UTC, comment #12: 

Julien,

I reviewed your patch and committed it here (http://hg.savannah.gnu.org/hgweb/octave/rev/ab75e72c5b36).  I followed it up with a changeset which modified the coding style to reflect core Octave (http://hg.savannah.gnu.org/hgweb/octave/rev/7004c733412f).  I'm attaching the diff file to this bug report which I annotated with the reason for each change.  I think that will make it useful for learning.  It should be straightforward to do kurtosis now, although I would use the correction factor as shown by the Mathworks in their documentation (http://www.mathworks.com/help/stats/kurtosis.html).

Cheers,
Rik

(file #29303)

Rik <rik5>
Group administrator
Sun 06 Oct 2013 10:35:45 AM UTC, comment #11: 

Hello Rik,

I can make a similar patch for kurtosis, if this one is OK.

Please tell me whenever you have time to look at it.

@++
Julien

Julien Bect <jbect>
Wed 02 Oct 2013 09:19:03 PM UTC, comment #10: 

Here is another version of the same patch, with a better documentation (the TeX output looks ok now).


(file #29288)

Julien Bect <jbect>
Wed 02 Oct 2013 06:26:26 PM UTC, comment #9: 

Ok, I'm trying now to improve the TeX documentation.

Julien Bect <jbect>
Wed 02 Oct 2013 05:20:14 PM UTC, comment #8: 

Julien,

The TeX part shows up in:
octave/doc/interpreter/octave.pdf

This is the Octave Manual in PDF form.

Michael Godfrey <godfrey>
Group Member
Wed 02 Oct 2013 11:11:04 AM UTC, comment #7: 

Ok, here is a first version of my patch for skewness.m, for review.

It solves several problems:

1) Add the missing 'flag' argument for Matlab compatibility

2) Return NaN when the standard deviation is zero (also for Matlab compatibility)

3) Deal when the type of the input argument (single/double) in a more transparent manner, in order to make it easir to solve bug #32168 (related to dataframes).

I have updated the documentation, but I probably need some help to make it right, since this is my first attempt at writing texinfo doc. In particular, I have tried to update the tex-specific part of the doc, but I don't known where I should check that the result is correct...

@++
Julien




(file #29282)

Julien Bect <jbect>
Tue 01 Oct 2013 10:53:15 PM UTC, comment #6: 

Execellent, I'll review it whenever it's ready.

Rik <rik5>
Group administrator
Tue 01 Oct 2013 09:25:21 PM UTC, comment #5: 

Hello Rik,

I have a patch almost ready for skewness.m that fixes the problem discussed here, fixes another compatibility issue (skewness should return NaN when the standard deviation is zero) and should make it easy to close bug #32168 at the same time.

I just need to update the documentation, I'll upload the patch in a few days at most.

@++
Julien

Julien Bect <jbect>
Thu 26 Sep 2013 10:34:06 PM UTC, comment #4: 

I had simply forgotten about this.  I'll see if I can fix this in within the next week or so.

Rik <rik5>
Group administrator
Thu 26 Sep 2013 07:54:33 PM UTC, comment #3: 

I checked now with version 3.7.7+ and the problem is still present. The status is "Need Info". Is there more information that I can provide?

Alexander Barth <abarth>
Sun 01 Jul 2012 07:30:00 PM UTC, comment #2: 

Here are the results with bias-correction:

>> x = [1:5].^2;
>> skewness (x, 0)


ans =

    0.6996

>> kurtosis (x, 0)


ans =

    2.3577

Thanks and best regards,
Alexander

Alexander Barth <abarth>
Fri 29 Jun 2012 06:13:30 PM UTC, comment #1: 

I swear I updated this a few years ago.  Maybe Matlab has been slowly changing their defaults.

On reading the documentation I see that skewness and kurtosis have a flag which decides whether bias correction is applied.  Rather than just changing Octave to use the current default, it would be better to add support for the bias flag and then set the flag to the default value used by Matlab.

Could you run the following tests to see how Matlab is applying bias correction?


x = [1:5].^2;
skewness (x, 0)
kurtosis (x, 0)


Rik <rik5>
Group administrator
Sun 24 Jun 2012 08:46:06 PM UTC, original submission:  

In matlab (R2010a), I get:

>> x = [1:5].^2;
>> moment(x,3)

ans =
  303.6000

>> skewness(x)

ans =
    0.4693

% this is (std(x,1)^(-3)) * mean ((x - mean(x)).^3)

>> kurtosis(x)

ans =
    1.8394

while in octave 3.6.0:

moment(x,3) 
ans =  4103

In matlab, moments are central while octave uses raw moment per default (central moments can be activated with moment(x,3,'c'))

skewness(x)
ans =  0.33580

This is (std(x,0)^(-3)) * mean ((x - mean(x)).^3). The normalization of std is thus different.

kurtosis(x)
ans = -1.8228

Again a different normalization of std, and a minus 3 in octave (excess kurtosis).

Concerning the normalization of std (1/n versus 1/(n-1)), I have seen at least a third definition on the web.
(http://www.itl.nist.gov/div898/handbook/eda/section3/eda35b.htm). What a mess!

I would propose that we follow the matlab convention here.

Attached are some patches, with the necessary changes in the calculations to the functions (but help text and test are not updated yet, ... bed time now).

Alexander Barth <abarth>

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #29418:  jbPatchKurtosis3.patch added by jbect (7KiB - text/x-patch)
file #29407:  jbPatchKurtosis2.patch added by jbect (4KiB - text/x-patch)
file #29401:  jbPatchKurtosis.patch added by jbect (4KiB - text/x-patch)
file #29309:  jbPatchSkewness6.patch added by jbect (969B - text/x-patch)
file #29306:  jbPatchSkewness5a.patch added by jbect (1KiB - text/x-patch)
file #29307:  jbPatchSkewness5b.patch added by jbect (1KiB - text/x-patch)
file #29308:  jbPatchSkewness5c.patch added by jbect (855B - text/x-patch)
file #29303:  diff.annotate added by rik5 (8KiB - application/octet-stream)
file #29288:  jbPatchSkewness4.patch added by jbect (6KiB - text/x-patch)
file #29282:  jbPatchSkewness3.patch added by jbect (6KiB - text/x-patch)
file #26100:  kurtosis.m.patch added by abarth (708B - text/x-patch)
file #26101:  skewness.m.patch added by abarth (668B - text/x-patch)
file #26102:  moment.m.patch added by abarth (751B - text/x-patch)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by godfrey (Posted a comment)
  • -email is unavailable- added by jbect (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by abarth (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 17 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2013-11-07 rik5 StatusNeed Info Fixed
        Open/ClosedOpen Closed
    2013-10-19 jbect Attached File- Added jbPatchKurtosis3.patch, #29418
    2013-10-18 jbect Attached File- Added jbPatchKurtosis2.patch, #29407
    2013-10-18 jbect Attached File- Added jbPatchKurtosis.patch, #29401
    2013-10-07 jbect Attached File- Added jbPatchSkewness6.patch, #29309
    2013-10-07 jbect Attached File- Added jbPatchSkewness5a.patch, #29306
        Attached File- Added jbPatchSkewness5b.patch, #29307
        Attached File- Added jbPatchSkewness5c.patch, #29308
    2013-10-07 rik5 Attached File- Added diff.annotate, #29303
    2013-10-02 jbect Attached File- Added jbPatchSkewness4.patch, #29288
    2013-10-02 jbect Attached File- Added jbPatchSkewness3.patch, #29282
    2012-06-29 rik5 StatusNone Need Info
        Summarymoment, skewness and kurtosis moment, skewness, and kurtosis differences with Matlab
    2012-06-24 abarth Attached File- Added kurtosis.m.patch, #26100
        Attached File- Added skewness.m.patch, #26101
        Attached File- Added moment.m.patch, #26102

    Back to the top

    Powered by Savane 3.13-4448.
    Corresponding source code