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

 
 

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

Submitted by:  Alexander Barth <abarth>
Submitted on:  Sun 24 Jun 2012 08:46:06 PM UTC  
 
Category: InterpreterSeverity: 3 - Normal
Priority: 5 - NormalItem Group: Matlab Compatibility
Status: FixedAssigned to: None
Originator Name: Alexander BarthOpen/Closed: Closed
Release: 3.6.0Operating System: GNU/Linux

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

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

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

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.

Rik <rik5>
Project 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

to

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

Rik <rik5>
Project 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>
Project 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>
Project 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>
Project 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>
Project 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>
Project 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>
Project 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>
Project 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?

Rik <rik5>
Project 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 File(s):
   
   
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
  • -unavailable- added by godfrey (Posted a comment)
  • -unavailable- added by jbect (Posted a comment)
  • -unavailable- added by rik5 (Posted a comment)
  • -unavailable- added by abarth (Submitted the item)
  •  

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

    Only project members can vote.

     

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follow 17 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Thu 07 Nov 2013 05:41:05 PM UTCrik5StatusNeed Info=>Fixed
      Open/ClosedOpen=>Closed
    Sat 19 Oct 2013 10:13:19 AM UTCjbectAttached File-=>Added jbPatchKurtosis3.patch, #29418
    Fri 18 Oct 2013 05:20:27 PM UTCjbectAttached File-=>Added jbPatchKurtosis2.patch, #29407
    Fri 18 Oct 2013 11:29:43 AM UTCjbectAttached File-=>Added jbPatchKurtosis.patch, #29401
    Mon 07 Oct 2013 07:12:19 PM UTCjbectAttached File-=>Added jbPatchSkewness6.patch, #29309
    Mon 07 Oct 2013 01:28:00 PM UTCjbectAttached File-=>Added jbPatchSkewness5a.patch, #29306
      Attached File-=>Added jbPatchSkewness5b.patch, #29307
      Attached File-=>Added jbPatchSkewness5c.patch, #29308
    Mon 07 Oct 2013 02:58:03 AM UTCrik5Attached File-=>Added diff.annotate, #29303
    Wed 02 Oct 2013 09:19:03 PM UTCjbectAttached File-=>Added jbPatchSkewness4.patch, #29288
    Wed 02 Oct 2013 11:11:04 AM UTCjbectAttached File-=>Added jbPatchSkewness3.patch, #29282
    Fri 29 Jun 2012 06:13:30 PM UTCrik5StatusNone=>Need Info
      Summarymoment, skewness and kurtosis=>moment, skewness, and kurtosis differences with Matlab
    Sun 24 Jun 2012 08:46:06 PM UTCabarthAttached 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.1-cleanup1