Issue 27181: Add geometric mean to statistics module (original) (raw)
Created on 2016-06-02 12:24 by cool-RR, last changed 2022-04-11 14:58 by admin. This issue is now closed.
Messages (45)
Author: Raymond Hettinger (rhettinger) * 
Date: 2016-06-02 21:04
Steven, this seems like a reasonable suggestion (though I would expect someone else will immediately suggest a harmonic mean as well). Is this within the scope of what you were trying to do with the statistics module?
Author: Steven D'Aprano (steven.daprano) *
Date: 2016-06-02 22:04
On Thu, Jun 02, 2016 at 09:04:54PM +0000, Raymond Hettinger wrote:
Steven, this seems like a reasonable suggestion (though I would expect someone else will immediately suggest a harmonic mean as well). Is this within the scope of what you were trying to do with the statistics module?
Yes, I think it is reasonable too. I'll aim to get this in to 3.6.
Author: Ram Rachum (cool-RR) *
Date: 2016-06-03 06:00
To complicate things further...
I implemented a geometric mean on my own, and then I figured out what I really want is a weighted geometric mean, so I implemented that for myself. If you'd want to include that, that'll be cool. Actually I'm not sure if the goal of the statistics module is to be comprehensive or minimal. I'm hoping it's meant to be comprehensive. But then I'd guess there would be a lot of things you'd want to add except my little feature.
Author: Ram Rachum (cool-RR) *
Date: 2016-06-04 14:14
And of course, if the goal of the statistics module is to be comprehensive, one should ask himself what should be the difference between this new module and a mature statistics module like scipy.stats, and whether we should try to copy the features of off scipy.stats.
Author: Mark Dickinson (mark.dickinson) *
Date: 2016-06-09 08:23
Choice of algorithm is a bit tricky here. There are a couple of obvious algorithms that work mathematically but result in significant accuracy loss in an IEEE 754 floating-point implementation: one is exp(mean(map(log, my_numbers))), where the log calls can introduce significant loss of information, and the other is prod(x**(1./len(my_numbers)) for x in my_numbers), where the **(1./n) operation similarly discards information. A better algorithm numerically is prod(my_numbers)**(1./len(my_numbers)), but that's likely to overflow quickly for large datasets (and/or datasets containing large values).
I'd suggest something along the lines of prod(my_numbers)**(1./len(my_numbers)), but keeping track of the exponent of the product separately and renormalizing where necessary to avoid overflow.
There are also algorithms for improved accuracy in a product, along the same lines as the algorithm used in fsum. See e.g., the paper "Accurate Floating-Point Product and Exponentiation" by Stef Graillat. [1] (I didn't know about this paper: I just saw a reference to it in a StackOverflow comment [2], which reminded me of this issue.)
[1] http://www-pequan.lip6.fr/~graillat/papers/IEEE-TC-prod.pdf [2] http://stackoverflow.com/questions/37715250/safe-computation-of-geometric-mean
Author: Mark Dickinson (mark.dickinson) *
Date: 2016-06-09 09:24
On the other hand, apparently exp(mean(log(...))) is good enough for SciPy: its current implementation looks like this:
def gmean(a, axis=0): a, axis = _chk_asarray(a, axis) log_a = ma.log(a) return ma.exp(log_a.mean(axis=axis))
Author: Steven D'Aprano (steven.daprano) *
Date: 2016-06-09 11:59
On Thu, Jun 09, 2016 at 09:24:04AM +0000, Mark Dickinson wrote:
On the other hand, apparently
exp(mean(log(...)))is good enough for SciPy:
Hmm, well, I don't have SciPy installed, but I've found that despite their (well-deserved) reputation, numpy (and presumably scipy) often have rather naive algorithms that can lose accuracy rather spectacularly.
py> statistics.mean([1e50, 2e-50, -1e50, 2e-50]) 1e-50 py> np.mean(np.array([1e50, 2e-50, -1e50, 2e-50])) 5e-51
py> statistics.mean([1e50, 2e-50, -1e50, 2e-50]*1000) 1e-50 py> np.mean(np.array([1e50, 2e-50, -1e50, 2e-50]*1000)) 5.0000000000000002e-54
On the other hand, np is probably a hundred times (or more) faster, so I suppose accuracy/speed makes a good trade off.
Author: Mark Dickinson (mark.dickinson) *
Date: 2016-06-09 12:55
Hmm, well, I don't have SciPy installed, but I've found that despite their (well-deserved) reputation, numpy (and presumably scipy) often have rather naive algorithms that can lose accuracy rather spectacularly.
Agreed. And as Ram Rachum hinted, there seems little point aiming to duplicate things that already exist in the de facto standard scientific libraries. So I think there's a place for a non-naive carefully computed geometric mean in the std. lib. statistics module, but I wouldn't see the point of simply adding an exp-mean-log implementation (not that anyone is advocating that).
Author: Steven D'Aprano (steven.daprano) *
Date: 2016-07-09 05:49
Does anyone have any strong feeling about the name for these functions?
gmean and hmean;
geometric_mean and harmonic_mean
And "subcontrary_mean" is not an option :-)
Author: Raymond Hettinger (rhettinger) *
Date: 2016-07-09 05:57
I would like to see them spelled-out: geometric_mean and harmonic_mean
Author: Mark Dickinson (mark.dickinson) *
Date: 2016-07-09 09:15
I would like to see them spelled-out: geometric_mean and harmonic_mean
+1
Author: Roundup Robot (python-dev) 
Date: 2016-08-09 04:18
New changeset 9eb5edfcf604 by Steven D'Aprano in branch 'default': Issue27181 add geometric mean. https://hg.python.org/cpython/rev/9eb5edfcf604
Author: Ram Rachum (cool-RR) *
Date: 2016-08-09 06:40
Thanks for the patch Steven! I won't comment about the code because I don't know enough about these algorithms, but I'm thinking, since you also did a refactoring of the statistics module, maybe these should be two separate patches/commits so it'll be easy to see which part is the new feature and which part is moving existing code around?
Author: Ram Rachum (cool-RR) *
Date: 2016-08-09 06:44
Also... I like the detailed docstrings with the real-life examples! That stuff helps when coding and using an unfamiliar function (since I see the docs in a panel of my IDE), so I wish I'd see more detailed docstrings like these ones in the standard library. For geometric_mean, maybe I'd add one sentence that describes how the geometric mean is calculated.
Author: Steven D'Aprano (steven.daprano) *
Date: 2016-08-09 08:00
On Tue, Aug 09, 2016 at 06:44:22AM +0000, Ram Rachum wrote:
For
geometric_mean, maybe I'd add one sentence that describes how the geometric mean is calculated.
What do you mean? As in, the mathematical definition of geometric mean?
Or do you mean a one sentence description of the algorithm?
Author: Ram Rachum (cool-RR) *
Date: 2016-08-09 08:44
I meant the mathematical definition.
Author: Martin Panter (martin.panter) *
Date: 2016-08-12 01:04
Tests fail on a Power PC buildbot:
http://buildbot.python.org/all/builders/PPC64LE%20Fedora%203.x/builds/1476/steps/test/logs/stdio
FAIL: testExactPowers (test.test_statistics.Test_Nth_Root) (i=29, n=11)
Traceback (most recent call last): File "/home/shager/cpython-buildarea/3.x.edelsohn-fedora-ppc64le/build/Lib/test/test_statistics.py", line 1216, in testExactPowers self.assertEqual(self.nroot(x, n), i) AssertionError: 29.000000000000004 != 29
====================================================================== FAIL: testExactPowersNegatives (test.test_statistics.Test_Nth_Root) (i=-29, n=11)
Traceback (most recent call last): File "/home/shager/cpython-buildarea/3.x.edelsohn-fedora-ppc64le/build/Lib/test/test_statistics.py", line 1228, in testExactPowersNegatives self.assertEqual(self.nroot(x, n), i) AssertionError: -29.000000000000004 != -29
Author: Mark Dickinson (mark.dickinson) *
Date: 2016-08-12 07:44
What no patch for pre-commit review?!
For computing nth roots, it may be worth special-casing the case n=2: for floats, math.sqrt is likely to be faster and more precise than an ad-hoc algorithm. (Indeed, I'd expect it to be perfectly correctly rounded on the vast majority of current machines.)
Author: Steven D'Aprano (steven.daprano) *
Date: 2016-08-12 11:07
I thought about special-casing n=2 to math.sqrt, but as that's an implementation detail I can make that change at any time. According to my testing, math.pow(x, 0.5) is no worse than sqrt, so I'm not sure if there's any advantage to having yet another branch.
I'd be interested in special-casing n=3 to math.cbrt (if and when it exists) now that its a standard C99 function.
Author: Mark Dickinson (mark.dickinson) *
Date: 2016-08-12 11:45
According to my testing, math.pow(x, 0.5) is no worse than sqrt.
It certainly is worse than sqrt, both in terms of speed and accuracy. Whether the difference is enough to make it worth special-casing is another question, of course, and as you say, that can happen later.
Author: Steven D'Aprano (steven.daprano) *
Date: 2016-08-14 02:31
I've created a new issue to track the loss of accuracy on PowerPC: http://bugs.python.org/issue27761
Author: Mark Dickinson (mark.dickinson) *
Date: 2016-08-14 09:16
A failing case:
statistics.geometric_mean([0.7 for _ in range(5000)]) Traceback (most recent call last): File "/Users/mdickinson/Python/cpython-git/Lib/statistics.py", line 362, in float_nroot isinfinity = math.isinf(x) OverflowError: int too large to convert to float
During handling of the above exception, another exception occurred:
Traceback (most recent call last): File "", line 1, in File "/Users/mdickinson/Python/cpython-git/Lib/statistics.py", line 595, in geometric_mean s = 2p * _nth_root(2q, n) File "/Users/mdickinson/Python/cpython-git/Lib/statistics.py", line 346, in nth_root return _nroot_NS.float_nroot(x, n) File "/Users/mdickinson/Python/cpython-git/Lib/statistics.py", line 364, in float_nroot return _nroot_NS.bignum_nroot(x, n) File "/Users/mdickinson/Python/cpython-git/Lib/statistics.py", line 489, in bignum_nroot b = 2q * _nroot_NS.nroot(2r, n) File "/Users/mdickinson/Python/cpython-git/Lib/statistics.py", line 384, in nroot r1 = math.pow(x, 1.0/n) OverflowError: int too large to convert to float
Author: Ned Deily (ned.deily) *
Date: 2016-08-15 23:05
FTR, multiple platforms are failing in various ways, not just PPC64, so Issue27761 was expanded to cover them and has been marked as a "release blocker".
Author: STINNER Victor (vstinner) *
Date: 2016-08-16 20:22
Failure on s390x Debian 3.x:
http://buildbot.python.org/all/builders/s390x%20Debian%203.x/builds/1455/steps/test/logs/stdio
====================================================================== FAIL: testExactPowers (test.test_statistics.Test_Nth_Root) (i=29, n=11)
Traceback (most recent call last): File "/home/dje/cpython-buildarea/3.x.edelsohn-debian-z/build/Lib/test/test_statistics.py", line 1216, in testExactPowers self.assertEqual(self.nroot(x, n), i) AssertionError: 29.000000000000004 != 29
====================================================================== FAIL: testExactPowersNegatives (test.test_statistics.Test_Nth_Root) (i=-29, n=11)
Traceback (most recent call last): File "/home/dje/cpython-buildarea/3.x.edelsohn-debian-z/build/Lib/test/test_statistics.py", line 1228, in testExactPowersNegatives self.assertEqual(self.nroot(x, n), i) AssertionError: -29.000000000000004 != -29
Author: Roundup Robot (python-dev)
Date: 2016-08-16 20:22
New changeset 54288b160243 by Victor Stinner in branch 'default': Issue #27181: Skip tests known to fail until a fix is found https://hg.python.org/cpython/rev/54288b160243
Author: STINNER Victor (vstinner) *
Date: 2016-08-16 20:24
I would like to use buildbots to check for regressions, but I see a lot of red buildbots, so buildbots became useless :-/
I skipped failing test_statistics tests, since failures are known.
I put the priority to "release blocker".
I suggest to either revert the change or find a fix before 3.6b1.
Author: Kubilay Kocak (koobs)
Date: 2016-08-17 08:41
For posterity, the following failure was observed on all (9/10/11(current) FreeBSD buildbots:
====================================================================== FAIL: testFraction (test.test_statistics.Test_Nth_Root)
Traceback (most recent call last): File "/usr/home/buildbot/python/3.x.koobs-freebsd9/build/Lib/test/test_statistics.py", line 1247, in testFraction self.assertEqual(self.nroot(x**12, 12), float(x)) AssertionError: 1.1866666666666665 != 1.1866666666666668
Author: Mark Dickinson (mark.dickinson) *
Date: 2016-08-17 17:18
self.assertEqual(self.nroot(x**12, 12), float(x))AssertionError: 1.1866666666666665 != 1.1866666666666668
That looks like a case where the test should simply be weakened to an assertAlmostEqual with a suitable tolerance; there's no strong reason to expect that nroot will give a faithfully rounded result in this case or any other.
Author: Steven D'Aprano (steven.daprano) *
Date: 2016-09-12 02:50
As discussed with Ned by email, I'm currently unable to build 3.6 and won't have time to work on this before b1. As discussed on #27761 my tests here are too strict and should be loosened, e.g. from assertEqual to assertAlmostEqual. Ned wrote:
"If you are only planning to make changes to the tests themselves, I think that can wait for b2."
I have no plans to change the publicly visible interface of geometric_mean.
Author: Mark Dickinson (mark.dickinson) *
Date: 2016-09-12 15:35
Steven: any thoughts about the
statistics.geometric_mean(0.7 for _ in range(5000))
failure? Should I open a separate bug report for that, or would you rather address it as part of this issue?
Author: STINNER Victor (vstinner) *
Date: 2016-09-12 15:39
statistics.geometric_mean([0.7 for _ in range(5000)]) Traceback (most recent call last): File "/Users/mdickinson/Python/cpython-git/Lib/statistics.py", line 362, in float_nroot isinfinity = math.isinf(x) OverflowError: int too large to convert to float
=> see also issue #27975
Author: Steven D'Aprano (steven.daprano) *
Date: 2016-09-13 02:28
On Mon, Sep 12, 2016 at 03:35:14PM +0000, Mark Dickinson wrote:
statistics.geometric_mean(0.7 for _ in range(5000))
I've raised a new ticket #28111
Author: Steven D'Aprano (steven.daprano) *
Date: 2016-10-04 16:14
I'm sorry to say that due to technical difficulties, geometric mean is not going to be in a fit state for beta 2 of 3.6, and so is going to be removed and delayed until 3.7.
Author: Roundup Robot (python-dev)
Date: 2016-10-04 16:25
New changeset 9dce0e41bedd by Steven D'Aprano in branch 'default': Issue #27181 remove geometric_mean and defer for 3.7. https://hg.python.org/cpython/rev/9dce0e41bedd
Author: Roundup Robot (python-dev)
Date: 2016-10-04 18:52
New changeset de0fa478c22e by Steven D'Aprano in branch '3.6': Issue #27181 remove geometric_mean and defer for 3.7. https://hg.python.org/cpython/rev/de0fa478c22e
Author: Ned Deily (ned.deily) *
Date: 2016-10-04 18:56
Thanks, Steven. Actually, we needed to remove geometric_mean from the 3.6 branch, not the default branch (which will become 3.7). I backported your removal patch to 3.6. Feel free to reapply geometric_mean to the default branch at your leisure.
Author: Cheryl Sabella (cheryl.sabella) *
Date: 2017-08-27 14:05
I was wondering if this has been taken up again for 3.7? Thanks!
Author: Cheryl Sabella (cheryl.sabella) *
Date: 2019-02-16 00:25
Updating the version in case this wanted to be considered for 3.8.
Author: Raymond Hettinger (rhettinger) *
Date: 2019-02-16 20:55
Updating the version in case this wanted to be considered for 3.8.
Yes. It would be nice to get this wrapped-up.
Author: Raymond Hettinger (rhettinger) *
Date: 2019-03-24 18:24
Almost three years have passed.
In the spirit of "perfect is the enemy of good", would it be reasonable to start with a simple, fast implementation using exp-mean-log? Then if someone wants to make it more accurate later, they can do so.
In some quick tests, I don't see much of an accuracy loss. It looks to be plenty good enough to use as a starting point:
--- Accuracy experiments ---
from decimal import Decimal from functools import reduce from operator import mul from random import expovariate, triangular from statistics import fmean
https://www.wolframalpha.com/input/?i=geometric+mean+12,+17,+13,+5,+120,+7
data = [12, 17, 13, 5, 120, 7] print(reduce(mul, map(Decimal, data)) ** (Decimal(1) / len(data))) 14.94412420173971227234687688 exp(fmean(map(log, map(fabs, data)))) 14.944124201739715
data = [expovariate(50.0) for i in range(1_000)] print(reduce(mul, map(Decimal, data)) ** (Decimal(1) / len(data))) 0.01140902688569587677205587938 exp(fmean(map(log, map(fabs, data)))) 0.011409026885695879
data = [triangular(2000.0, 3000.0, 2200.0) for i in range(10_000)] print(reduce(mul, map(Decimal, data)) ** (Decimal(1) / len(data))) 2388.381301718524160840023868 exp(fmean(map(log, map(fabs, data)))) 2388.3813017185225
data = [lognormvariate(20.0, 3.0) for i in range(100_000)] min(data), max(data) (2421.506538652375, 137887726484094.5) print(reduce(mul, map(Decimal, data)) ** (Decimal(1) / len(data))) 484709306.8805352290183838500 exp(fmean(map(log, map(fabs, data)))) 484709306.8805349
Author: Steven D'Aprano (steven.daprano) *
Date: 2019-03-28 23:37
In the spirit of "perfect is the enemy of good", would it be reasonable to start with a simple, fast implementation using exp-mean-log? Then if someone wants to make it more accurate later, they can do so.
I think that is a reasonable idea. On the basis that something is better than nothing, go ahead. We can discuss accuracy and speed issues later.
Getting some tricky cases down for reference:
older (removed) implementation
py> geometric_mean([7]*2) 7.0 py> geometric_mean([7]*15) 7.0
Raymond's newer (faster) implementation
py> exp(fmean(map(log, [7]*2))) 6.999999999999999 py> exp(fmean(map(log, [7]*15))) 6.999999999999999
py> geometric_mean([3,27]) 9.0 py> geometric_mean([3,27]*5) 9.0
py> exp(fmean(map(log, [3,27]))) 9.000000000000002 py> exp(fmean(map(log, [3,27]*5))) 8.999999999999998
py> x = 2.5e15 py> geometric_mean([x]*100) 2500000000000000.0 py> exp(fmean(map(log, [x]*100))) 2499999999999999.5
On the other hand, sometimes rounding errors work in our favour:
py> geometric_mean([1e50, 1e-50]) # people might expect 1.0 0.9999999999999998 py> 1e-50 == 1/(1e50) # even though they aren't quite inverses False
py> exp(fmean(map(log, [1e50, 1e-50]))) 1.0
Author: Raymond Hettinger (rhettinger) *
Date: 2019-03-29 03:55
On the basis that something is better than nothing, go ahead. We can discuss accuracy and speed issues later.
Thanks. I'll put together a PR for your consideration.
Author: Raymond Hettinger (rhettinger) *
Date: 2019-04-02 07:13
Steven, how does this look?
https://patch-diff.githubusercontent.com/raw/python/cpython/pull/12638.diff
Author: Raymond Hettinger (rhettinger) *
Date: 2019-04-07 16:20
New changeset 6463ba3061bd311413d2951dc83c565907e10459 by Raymond Hettinger in branch 'master': bpo-27181: Add statistics.geometric_mean() (GH-12638) https://github.com/python/cpython/commit/6463ba3061bd311413d2951dc83c565907e10459
Author: Raymond Hettinger (rhettinger) *
Date: 2019-04-07 16:21
Feel free to reopen this if something further needed to be changed or discussed.
History
Date
User
Action
Args
2022-04-11 14:58:31
admin
set
github: 71368
2019-04-07 16:21:07
rhettinger
set
status: open -> closed
resolution: fixed
messages: +
stage: patch review -> resolved
2019-04-07 16:20:06
rhettinger
set
messages: +
2019-04-02 07:13:37
rhettinger
set
messages: +
2019-03-31 08:49:05
rhettinger
set
keywords: + patch
pull_requests: + <pull%5Frequest12570>
2019-03-29 03:55:43
rhettinger
set
messages: +
2019-03-29 00:06:14
vstinner
set
nosy: - vstinner
2019-03-28 23:37:04
steven.daprano
set
messages: +
2019-03-24 18:24:21
rhettinger
set
messages: +
2019-02-16 21:16:08
ned.deily
set
nosy: - ned.deily
2019-02-16 20:55:21
rhettinger
set
messages: +
2019-02-16 00:25:41
cheryl.sabella
set
messages: +
versions: + Python 3.8, - Python 3.7
2017-08-27 14:05:18
cheryl.sabella
set
nosy: + cheryl.sabella
messages: +
2016-10-04 18:56:20
ned.deily
set
messages: +
stage: patch review
2016-10-04 18:52:26
python-dev
set
messages: +
2016-10-04 16:25:29
python-dev
set
messages: +
2016-10-04 16:14:43
steven.daprano
set
priority: release blocker ->
messages: +
versions: + Python 3.7, - Python 3.6
2016-09-13 02:28:03
steven.daprano
set
messages: +
2016-09-12 15:39:40
vstinner
set
messages: +
2016-09-12 15:35:13
mark.dickinson
set
messages: +
2016-09-12 02:50:47
steven.daprano
set
messages: +
2016-08-17 17🔞15
mark.dickinson
set
messages: +
2016-08-17 08:41:31
koobs
set
nosy: + koobs
messages: +
2016-08-16 20:24:35
vstinner
set
priority: normal -> release blocker
messages: +
2016-08-16 20:22:57
python-dev
set
messages: +
2016-08-16 20:22:07
vstinner
set
nosy: + vstinner
messages: +
2016-08-15 23:05:40
ned.deily
set
nosy: + ned.deily
messages: +
2016-08-14 09:16:11
mark.dickinson
set
messages: +
2016-08-14 02:31:20
steven.daprano
set
messages: +
2016-08-12 11:45:25
mark.dickinson
set
messages: +
2016-08-12 11:07:42
steven.daprano
set
messages: +
2016-08-12 07:44:40
mark.dickinson
set
messages: +
2016-08-12 01:04:21
martin.panter
set
nosy: + martin.panter
messages: +
2016-08-09 08:44:08
cool-RR
set
messages: +
2016-08-09 08:00:07
steven.daprano
set
messages: +
2016-08-09 06:44:22
cool-RR
set
messages: +
2016-08-09 06:40:53
cool-RR
set
messages: +
2016-08-09 04🔞51
python-dev
set
nosy: + python-dev
messages: +
2016-07-09 09:15:14
mark.dickinson
set
messages: +
2016-07-09 05:57:46
rhettinger
set
messages: +
2016-07-09 05:49:17
steven.daprano
set
messages: +
2016-06-09 12:55:29
mark.dickinson
set
messages: +
2016-06-09 11:59:05
steven.daprano
set
messages: +
2016-06-09 09:24:04
mark.dickinson
set
messages: +
2016-06-09 08:23:34
mark.dickinson
set
nosy: + mark.dickinson
messages: +
2016-06-04 14:14:26
cool-RR
set
messages: +
2016-06-03 06:00:58
cool-RR
set
messages: +
2016-06-02 22:04:12
steven.daprano
set
messages: +
2016-06-02 21:04:54
rhettinger
set
assignee: steven.daprano
messages: +
nosy: + rhettinger
2016-06-02 12:38:33
xiang.zhang
set
nosy: + steven.daprano
2016-06-02 12:24:07
cool-RR
create