classification
Title: Fix awkwardness of statistics.mode() for multimodal datasets
Type: behavior Stage: resolved
Components: Library (Lib) Versions: Python 3.8
process
Status: closed Resolution: fixed
Dependencies: Superseder:
Assigned To: steven.daprano Nosy List: Windson Yang, cheryl.sabella, francismb, rhettinger, scotchka, steven.daprano
Priority: normal Keywords: patch

Created on 2019-02-03 18:51 by rhettinger, last changed 2019-03-12 07:44 by rhettinger. This issue is now closed.

Pull Requests
URL Status Linked Edit
PR 12089 merged rhettinger, 2019-02-28 08:14
Messages (25)
msg334796 - (view) Author: Raymond Hettinger (rhettinger) * (Python committer) Date: 2019-02-03 18:51
The current code for mode() does a good deal of extra work to support its two error outcomes (empty input and multimodal input).  That latter case is informative but doesn't provide any reasonable way to find just one of those modes, where any of the most popular would suffice.  This arises in nearest neighbor algorithms for example. I suggest adding an option to the API:

   def mode(seq, *, first_tie=False):       
       if tie_goes_to_first:
           # CHOOSE FIRST x ∈ S | ∄ y ∈ S : x ≠ y ∧ count(y) > count(x)
           return return Counter(seq).most_common(1)[0][0]
       ...

Use it like this:

    >>> data = 'ABBAC'
    >>> assert mode(data, first_tie=True) == 'A'

With the current API, there is no reasonable way to get to 'A' from 'ABBAC'.

Also, the new code path is much faster than the existing code path because it extracts only the 1 most common using min() rather than the n most common which has to sort the whole items() list.  New path: O(n).  Existing path: O(n log n).

Note, the current API is somewhat awkward to use.  In general, a user can't know in advance that the data only contains a single mode.  Accordingly, every call to mode() has to be wrapped in a try-except.  And if the user just wants one of those modal values, there is no way to get to it.  See https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.mode.html for comparison.

There may be better names for the flag.  "tie_goes_to_first_encountered" seemed a bit long though ;-)
msg335147 - (view) Author: Francis MB (francismb) * Date: 2019-02-10 10:15
>> There may be better names for the flag.  "tie_goes_to_first_encountered" seemed a bit long though ;-)

Could it may be an alternative to set the mode tie case in a form like:

def mode(seq, *, case=CHOOSE_FIRST):
  [...]

(or TIE_CHOOSE_FIRST, ...) where CHOOSE_FIRST is an enum ?

Thanks!
msg335182 - (view) Author: Steven D'Aprano (steven.daprano) * (Python committer) Date: 2019-02-11 00:12
Thanks Raymond for the interesting use-case.

The original design of mode() was support only the basic form taught in 
secondary schools, namely a single unique mode for categorical data or 
discrete numerical data.

I think it is time to consider a richer interface to support more uses, 
such as estimating the mode of continuous numerical data, and the 
multi-mode case you bring up.

One reason I've been hesitant is that deciding what is the right 
behaviour is quite hard (or at least it has been for me). I think there 
are a few cases to consider:

- the existing behaviour (which may not be very useful?) which is to
  raise an exception unless the mode is unique; 

- would it be more useful to return an out-of-band value like a NAN 
  or None?

- the multi-mode case where you want some arbitrary(?) mode, possibly
  the left-most (smallest) for numeric data;

- the multi-mode case where you want all the modes.

I like Francis' suggestion to use an enum to select the behavioural, er, 
mode (pun intended). What do you think?
msg335679 - (view) Author: Raymond Hettinger (rhettinger) * (Python committer) Date: 2019-02-16 10:48
I would stick with "first_tie=False".  That caters to the common case, avoids API complications, and does something similar to what other stats packages are doing.

----- API Survey ----

Maple: """This function is only guaranteed to return one potential mode - in cases where multiple modes exist, only the first detected mode is guaranteed to be returned.""" -- https://www.maplesoft.com/support/help/maple/view.aspx?path=Statistics%2fMode

R: 'R does not have a standard in-built function to calculate mode.' -- https://www.tutorialspoint.com/r/r_mean_median_mode.htm

Matlab: "When there are multiple values occurring equally frequently, mode returns the smallest of those values. ... If A is an empty 0-by-0 matrix, mode(A) returns NaN." -- https://www.mathworks.com/help/matlab/ref/mode.html

Mathematica: "The mode of a set of data is implemented in the Wolfram Language as Commonest[data]. ... When several elements occur with equal frequency, Commonest[list,…] picks first the ones that occur first in list." -- http://mathworld.wolfram.com/Mode.html and https://reference.wolfram.com/language/ref/Commonest.html

SciPy: "If there is more than one such value, only the smallest is returned." -- https://docs.scipy.org/doc/scipy-0.19.1/reference/generated/scipy.stats.mode.html

SAS: "Most frequent value (if not unique, the smallest mode)" -- https://documentation.sas.com/?docsetId=qcug&docsetTarget=qcug_capability_sect225.htm&docsetVersion=14.2&locale=en
msg335680 - (view) Author: Francis MB (francismb) * Date: 2019-02-16 11:28
Good point Raymond!

Only a minor observation on the packages API:  


[1] SciPy: scipy.stats.mode(a, axis=0, nan_policy='propagate')
"Returns an array of the modal (most common) **value** in the passed array." --> Here it claims to return just ONE value

And use of different policies on parameters :
nan_policy : {‘propagate’, ‘raise’, ‘omit’}, optional
Defines how to handle when input contains nan. ‘propagate’ returns nan, ‘raise’ throws an error, ‘omit’ performs the calculations ignoring nan values. Default is ‘propagate’.

Equivalent one could say 'multivalue_policy'


[2] Matlab: Mode: "Most frequent **values** in array"

...returns the sample mode of A, which is the most frequently occurring *value* in A...."

IMHO it seems inconsistent *values* vs. *value* (or a doc-bug ?).


An a question:
Does it that mean that mode in that case really should potentially return an array of values, e.g. all the values with equal frequency?

In that case the user has the chance to get the first, the last or just all, ...


----
[1] https://docs.scipy.org/doc/scipy-0.19.1/reference/generated/scipy.stats.mode.html

[2] https://la.mathworks.com/help/matlab/ref/mode.html
msg335684 - (view) Author: Windson Yang (Windson Yang) * Date: 2019-02-16 13:36
IMHO, we don't need to add the option. We can return the smallest value from the **table** instead of the code below.

    if len(table) == 1:
        return table[0][0]

[1] https://github.com/python/cpython/blob/master/Lib/statistics.py#L502
msg335704 - (view) Author: Raymond Hettinger (rhettinger) * (Python committer) Date: 2019-02-16 17:44
> We can return the smallest value from the **table** instead
> of the code below.

Internally, that does too much work and then throws most of it away.

The difference between Counter(data).most_common()[1] and Counter(data).most_common(1) is that the former materializes the entire items iterator into a list and does a full sort, while the latter is memory friendly and makes a single pass with min().
msg335726 - (view) Author: Raymond Hettinger (rhettinger) * (Python committer) Date: 2019-02-16 22:52
I've been thinking about this a bit more.  ISTM that for now, it's best to get mode() working for everyday use, returning just the first mode encountered.  This keeps the signature simple (Iterable -> Scalar).

In the future, if needed, there is room to add separate functions with their own clean signatures and use cases.  For example, having a separate multimode() that always returns a list would be clearer and easier to use than setting a flag on the current scalar valued mode() function. 

Another example would be mode_estimate() that returns a scalar float
that estimates the mode given data sampled from a continuous variable. This warrants its own function because the signature, use cases, options, and implementation method are all completely different from a mode based on counting.

Categorical, binned, or ordinal data:

   mode(data: Iterable, *, first_tie=False) -> object
   multimode(data: Iterable) -> List[object]

Continuous data:

   mode(data: Iterable[Real]) -> Real
msg335746 - (view) Author: Windson Yang (Windson Yang) * Date: 2019-02-17 02:27
I only tested stats.mode() from scipy, data = 'BBAAC' should also return 'A'. But in your code **return return Counter(seq).most_common(1)[0][0]** will return B. Did I miss something?
msg335764 - (view) Author: Raymond Hettinger (rhettinger) * (Python committer) Date: 2019-02-17 09:10
> Did I miss something?

Yes.  It doesn't really matter which mode is returned as long as it is deterministically chosen.  We're proposing to return the first mode rather than the smallest mode.  

Scipy returns the smallest mode because that is convenient given that the underlying operation is np.unique() which returns unique values in sorted order [1]. 

We want to return the first mode encountered because that is convenient given that the underlying operation is max() which returns the first maximum value it encounters.

Another advantage of return-first rather than return-smallest is that our mode() would work for data values that are hashable but not orderable (i.e. frozensets).

[1] https://github.com/scipy/scipy/blob/v0.19.1/scipy/stats/stats.py#L378-L443
msg335785 - (view) Author: Raymond Hettinger (rhettinger) * (Python committer) Date: 2019-02-17 19:27
The attraction to having a first_tie=False flag is that it is backwards compatible with the current API.

However on further reflection, people would almost never want the existing behavior of raising an exception rather than returning a useful result.

So, we could make first_tie the default behavior and not have a flag at all.  We would also add a separate multimode() function with a clean signature, always returning a list.  FWIW, MS Excel evolved to this solution as well (it has MODE.SGNL and MODE.MULT).

This would give us a clean, fast API with no flags:

    mode(Iterable) -> scalar
    multimode(Iterable) -> list  

ISTM this is an all round win:

* Fixes the severe usability problems with mode().  Currently, it can't be used without a try/except. The except-clause can't distinguish between an empty data condition and no unique mode condition, nor can the except clause access the dataset counts to get to a useful answer.

* It makes mode() memory friendly and faster.  Currently, mode() has to convert the counter items into a list() and run a full sort.  However, if we only need the first mode, we can feed the items iterator to max() and make only a single pass, never having to create a list.

However, there would be an incompatibility in the API change.  It is possible that a person really does want an exception instead of a value when there are multiple modes.  In that case, this would break their code so that they have to switch to multimode() as a remediation.  OTOH, it is also possible that people have existing code that uses mode() but their code has a latent bug because they weren't using a try/except and haven't yet encountered a multi-modal dataset. 

So here are some options:

* Change mode() to return the first tie instead of raising an exception.  This is a behavior change but leaves you with the cleanest API going forward.

* Add a Deprecation warning to the current behavior of mode() when it finds multimodal data.  Then change the behavior in Python 3.9.  This is uncomfortable for one release, but does get us to the cleanest API and best performance.

* Change mode() to have a flag where first_tie defaults to False.  This is backwards compatible, but it doesn't fix latent bugs, it leaves us with on-going API complexities, and the default case remains the least usable option.

* Leave mode() as-is and just document its limitations.

For any of those options, we should still add a separate multimode() function.
msg335879 - (view) Author: Windson Yang (Windson Yang) * Date: 2019-02-19 02:27
I think right now we can 
> Change mode() to return the first tie instead of raising an exception.  This is a behavior change but leaves you with the cleanest API going forward.

as well as 
> Add a Deprecation warning to the current behavior of mode() when it finds multimodal data.  

In Python 3.9, we change the behavior.
msg336007 - (view) Author: Francis MB (francismb) * Date: 2019-02-19 19:35
>> [...] This keeps the signature simple (Iterable -> Scalar). [...]
>>
>> Categorical, binned, or ordinal data:
>>
>>   mode(data: Iterable, *, first_tie=False) -> object
>>   multimode(data: Iterable) -> List[object]

This seems reasonable to me due legacy (although I really thing that
multimode is just the real thing :-) )

>> Continuous data:
>>   mode(data: Iterable[Real]) -> Real

What should return in that case: E.g.: mode([42.0, 42.0, 42.0, 1.0, 1.0, 1.0]) ?

42.0 ? or 1.0 ? or [42.0, 1.0] ? or do I have misunderstood something ? 


Thanks!
msg336332 - (view) Author: Francis MB (francismb) * Date: 2019-02-22 17:19
Good options itemization!

>> This would give us a clean, fast API with no flags:
>>   mode(Iterable) -> scalar
>>   multimode(Iterable) -> list
[...]
>> For any of those options, we should still add a separate multimode()
>> function.
[..]
>> * Add a Deprecation warning to the current behavior of mode() when it
>> finds multimodal data.  Then change the behavior in Python 3.9.  This
>> is uncomfortable for one release, but does get us to the cleanest API
>> and best performance.

LGTM upto 3.9, shouldn't be "mode" at some point be replaced by
"multimode" ?
msg336338 - (view) Author: Raymond Hettinger (rhettinger) * (Python committer) Date: 2019-02-22 17:48
> shouldn't be "mode" at some point be replaced by "multimode" ?

No.  The signature is completely incompatible with the existing mode() function.

Like MS Excel which has two functions, MODE.SGNL and MODE.MULT, we should also have two functions, each with a clean signature and with running speed that is optimal for the desired result:

   mode(data: Iterable) -> object
       Returns a single value
       If multiple modes found, return first encountered
       Raises StatisticsError for an empty input
       Fast O(n), single pass, doesn't keep all data in memory

   multimode(data: Iterable) -> List[object]
       Always returns a list
       If multiple modes found, all are returned
       Returns an empty list for an empty input
       Slow O(n log n), loads all data in memory, full sort

For the first one, I recommend skipping deprecation and just changing the behavior that usually raises an exception for multiple modes.  In its current form, it likely to create bugs due to uncaught exceptions, and it doesn't provide any work-around if you do want only one of the modes.

By analogy, consider what would happen if we had a max() function that raised an exception when there were duplicate maximum values.  It would be almost usable.  So, what the real max() actually does is return the first encountered maximum value:

    >>> max(3, 1, 3.0)
    3
    >>> max(3.0, 1, 3)
    3.0
msg336570 - (view) Author: Steven D'Aprano (steven.daprano) * (Python committer) Date: 2019-02-25 23:12
Executive summary:

- let's change the behaviour of mode to return a single mode rather than raise an exception if there are multimodes;
- and let's do it without a depreciation period.


Further comments in no particular order:

I agree that in practice the current implementation of mode has turned out to be quite awkward.

Unfortunately I feel that just chosing one mode and returning it as "the" mode is statistically dubious (despite it being commonplace!). I've checked with some statisticians, or at least people who claim to be statisticians on the Internet, and they agree.

(The usual advice from statisticians is to always plot your data before talking about the mode or modes, as even local maxima can be meaningful.)

Nevertheless, Raymond's use-case of the nearest neighbor algorithm is persuasive. We should change the mode function to return one of the modes.

Which one? Again, people I've spoken to suggest that there's no statistical basis to picking any specific one (first found, last found, smallest, largest). I'm unhappy about documenting which mode is returned, as I expect it will be an implementation detail that could change in the future.

Raymond, are you happy to leave the specific detail of which mode is returned unspecified? Or do you have an argument for why we should specify it, and can do so without limiting future implementations?

I too strongly prefer a separate multimode function over a flag to the mode function.

As for backwards-compatibility, I've asked around trying to find anyone who claims that changing the behaviour would break their code. I haven't found anyone. I guess that means that unless any of the other core devs object, I'm happy to go with Raymond's recommendation to just change the behaviour without a depreciation period.
msg336573 - (view) Author: Steven D'Aprano (steven.daprano) * (Python committer) Date: 2019-02-25 23:19
What do people think about leaving this as an "Easy First Issue" for the sprint? If others agree that it is sufficiently easy, we can assign the task to Cheryl.

It should be fairly easy: mode calls an internal function _counts which is not public and not used anywhere else, so it can be changed completely or even discarded if needed.
msg336594 - (view) Author: Raymond Hettinger (rhettinger) * (Python committer) Date: 2019-02-26 01:58
>  If others agree that it is sufficiently easy, we can assign 
> the task to Cheryl.

It's only easy if we clearly specify what we want to occur.  Deciding what the right behavior should be is not a beginner skill.

Proposed spec:
'''
Modify the API statistics.mode to handle multimodal cases so that the first mode encountered is the one returned.  If the input is empty, raise a StatisticsError.

TestCases:
    mode([])   --> StatisticsError
    mode('aabbbcc') --> 'c'
    mode(iter('aabbbcc')) --> 'c'
    mode('eeffddddggaaaa') --> 'a'

Implementation:
    * Discard the internal _counts function.
    * Instead use Counter(data).most_common(1)[0][0]
      because that makes only a single pass over the data

Documentation:
    * Update statistics.rst and include a versionchanged directive

    * In the Whatsnew compatibility section, note this is a behavior change.
      Code that used to raise StatisticsError will now return a useful value.
      Note that the rationale for the change is that the current mode()
      behavior would unexpectedly fail when given multimodal data.

When: 
    * We want this for 3.8 so it can't wait very long
'''
msg336674 - (view) Author: Steven D'Aprano (steven.daprano) * (Python committer) Date: 2019-02-26 14:19
> Proposed spec:
> '''
> Modify the API statistics.mode to handle multimodal cases so that the 
> first mode encountered is the one returned.  If the input is empty, 
> raise a StatisticsError.

Are you happy guaranteeing that it will always be the first mode 
encountered? I'm not happy about it, but I'll follow you lead on this 
one.

> TestCases:
>     mode([])   --> StatisticsError
>     mode('aabbbcc') --> 'c'

That should be 'b'.

>     mode(iter('aabbbcc')) --> 'c'

And again 'b'.

>     mode('eeffddddggaaaa') --> 'a'

If it is first encountered, that should be 'd'.
msg336733 - (view) Author: Raymond Hettinger (rhettinger) * (Python committer) Date: 2019-02-27 07:16
> Are you happy guaranteeing that it will always be the first
> mode encountered?

Yes.  

All of the other implementations I looked at make some guarantee about which mode is returned.  Maple, Matlab, and Excel all return the first encountered.¹  That is convenient for us because it is what Counter(data).most_common(1) already does and does cheaply (single pass, no auxiliary memory).  It also matches what a number of our other tools do:

>>> max(3, 3.0)       # 3 is first encountered
3
>>> max(3.0, 3)       # 3.0 is first encountered
3.0
>>> list(dict.fromkeys('aabbaacc'))[0] # 'a' is first encountered
'a'
>>> sorted([3, 3.0])[0]  # 3 is first encountered (due to sort stability)
3
>>> sorted([3.0, 3])[0]  # 3.0 is first encountered (due to sort stability)
3.0

¹ Scipy returned the smallest value rather than first value but it algorithm was sorting based to accommodate running a parallel mode() computation on multiple columns of an array. For us, that approach would be much slow, would require more memory, and would require more bookkeeping. 

P.S. I'm no longer thinking this should be saved for Pycon sprints.  That is right at the beta 1 feature freeze.  We should aim to get this resolved well in advance of that.
msg336815 - (view) Author: Raymond Hettinger (rhettinger) * (Python committer) Date: 2019-02-28 08:17
Attached a draft PR for discussion purposes.  Let me know what you think (I'm not wedded to any part of it).
msg337626 - (view) Author: Raymond Hettinger (rhettinger) * (Python committer) Date: 2019-03-10 18:08
Steven, are you okay with applying this PR so we can put this to rest, cleanly and permanently?
msg337638 - (view) Author: Raymond Hettinger (rhettinger) * (Python committer) Date: 2019-03-10 20:28
Here's a text only link to the patch:  https://patch-diff.githubusercontent.com/raw/python/cpython/pull/12089.patch
msg337659 - (view) Author: Steven D'Aprano (steven.daprano) * (Python committer) Date: 2019-03-11 12:01
Looks good to me, I'm happy to accept it.

Thank you for your efforts Raymond, can I trouble you to do the merge yourself please, I'm still having issues using the Github website.
msg337724 - (view) Author: Raymond Hettinger (rhettinger) * (Python committer) Date: 2019-03-12 07:43
New changeset fc06a192fdc44225ef1cc879f615a81931ad0a85 by Raymond Hettinger in branch 'master':
bpo-35892: Fix mode() and add multimode() (#12089)
https://github.com/python/cpython/commit/fc06a192fdc44225ef1cc879f615a81931ad0a85
History
Date User Action Args
2019-03-12 07:44:04rhettingersetstatus: open -> closed
resolution: fixed
stage: patch review -> resolved
2019-03-12 07:43:37rhettingersetmessages: + msg337724
2019-03-11 12:01:27steven.dapranosetmessages: + msg337659
2019-03-10 20:28:46rhettingersetmessages: + msg337638
2019-03-10 20:27:25scotchkasetnosy: + scotchka
2019-03-10 18:08:49rhettingersetmessages: + msg337626
2019-02-28 08:17:02rhettingersetmessages: + msg336815
2019-02-28 08:14:55rhettingersetkeywords: + patch
stage: patch review
pull_requests: + pull_request12099
2019-02-27 07:16:35rhettingersetmessages: + msg336733
2019-02-26 14:19:22steven.dapranosetmessages: + msg336674
2019-02-26 01:58:06rhettingersetmessages: + msg336594
2019-02-25 23:19:21steven.dapranosetnosy: + cheryl.sabella
messages: + msg336573
2019-02-25 23:12:37steven.dapranosetmessages: + msg336570
2019-02-22 17:48:56rhettingersetmessages: + msg336338
2019-02-22 17:19:09francismbsetmessages: + msg336332
2019-02-19 19:35:28francismbsetmessages: + msg336007
2019-02-19 02:27:03Windson Yangsetmessages: + msg335879
2019-02-17 19:27:53rhettingersetmessages: + msg335785
2019-02-17 09:10:34rhettingersetmessages: + msg335764
2019-02-17 02:27:07Windson Yangsetmessages: + msg335746
2019-02-16 22:52:05rhettingersetmessages: + msg335726
2019-02-16 17:44:35rhettingersetmessages: + msg335704
2019-02-16 13:36:18Windson Yangsetnosy: + Windson Yang
messages: + msg335684
2019-02-16 11:28:23francismbsetmessages: + msg335680
2019-02-16 10:48:06rhettingersetmessages: + msg335679
2019-02-11 00:12:13steven.dapranosetmessages: + msg335182
2019-02-10 10:15:23francismbsetnosy: + francismb
messages: + msg335147
2019-02-03 18:51:25rhettingercreate