ENH: stats.circmedian: function for calculating the circular median - #24166

I concur in the choice of arc-distance as the preferred default. It's the most convenient way to link to "medians" in other spaces (minimizing the L1 distance of deviations, or its equivalent in the space).

convention works for me.

I re-read two articles, ¹ and ², to look for results on the structure of angles that minimize the total arc length.

It is shown in [1] (Lemma 5) that there exists at least one sample point that is a minimizer. In addition, if the sample size is odd, then a minimizer is necessarily a sample point ([1], Theorem 4).

If the sample size is even, the results seem less clean. If the minimizers form a unique arc of length less than $pi$, it seems customary to take the mid-point of the arc as the circular median. Otherwise, [2] (Definition 4.1) proposes to pick an arbitrary sample point minimizer, say one with the smallest value.

This PR returns the circular mean of sample point minimizers. I couldn't find a result stating that the circular mean of sample point minimizers remains a minimizer if it is well-defined. However, this is probably true given several sources appear to do this. Update: This is false. See counterexample in the second comment below.

I'm open to argument names other than 'convention' (e.g. kind, variant, definition), but method is more commonly used when there are multiple ways of computing/approximating a quantity with a single definition. convention is what we used to distinguish between Pearson and Fisher kurtosis in the new RV infrastructure.

I think "variant" might be more accurate here, since "arc-distance" and "bisecting" correspond to quite different definitions of median, and there are yet more variants (such as the "separable" and the "geometric" medians). "Convention" might be more appropriate for things that differ only in (constant) offset and scale, such as kg vs pound, 0- vs 1-based index, Fisher vs Pearson kurtosis, left-hand vs right-hand, etc. Here the various definitions of median all try to convey the "preferred direction" but are otherwise loosely related apparently.

Separately, the handling of non-unique "bisecting" median seems tricky. For example, by eye-balling the illustrations of Figure 2 in ¹, it seems doubtful that the circular mean of bisecting sample points remain a bisecting point. And I didn't find any source that elaborate on the choice of a bisecting median when the candidates are non-unique. Since the "bisecting" median seems not widely supported in other statistical software, maybe would it be easier to only implement the "arc-distance" variant for now and keep the "variant" argument for future extension?

A third point of thought is about how to document the return value of the function. I think it could make life easier for the user if they are told what to expect for the return value instead of how the return value is computed. If the user asks "what is returned if the minimizer is not unique?", the current documentation would tell them "the circular mean of sample minimizers or NAN if such mean is not well-defined", which is a faithful description of the implementation but might arguably fall short of answering the question. I wonder if it's correct to assert that the return value (of the current implementation) is the mid point of all minimizers if they form a (single) arc of length less than $\pi$, and NAN otherwise. Update: It is not. See counterexample in the next comment.

Using a semi-exhaustive approach, I found a curious example where the current implementation (as well as that of R's circular package) produces a circular median that neither minimises the total arc distance nor bisects the datapoints.

The sample is 1, 2, 5, 9, 9 with a period of 12. (Think of those as hours on a 12-hour clock.)

Using a semi-exhaustive search, the arg min of total arc distance is found to be the set $\{1, 9\}$. The minimum total arc distance is 13. According to the current implementation, the circular mean of the sample subset $\{1, 9, 9\}$ is 10; the corresponding total arc distance is 14, larger than the minimum, and the diameter that connects 10 o'clock and 4 o'clock does not bisect the sample points (there are 3 points below and 2 points above).

This example shows that taking the circular mean of sample point minimisers may yield an answer that neither minimises the total arc distance nor bisects the data points, and hence is not a median. If we want the returned value to have the median property, a possible workaround is to return an arbitrary sample point minimiser when the sample size is odd. Alternatively, we could keep the current implementation and put a note in the documentation that the return value might not be a minimiser or bisector.

import numpy as np

def arc_distance_median(x: np.array, period: float) -> list[tuple[float, float]]:
    """Given an array `x` of angular measurements with period `period`, return
    a list of disjoint, closed intervals such that every direction within
    these intervals minimize the total arc distance to the measurements."""
    x = np.asarray(x)
    if len(x) == 0:
        raise [(0, period)]
    
    y = np.unique(x % period)
    if len(y) == 1:
        return [(y[0], y[0])]
    
    y = np.unique(np.concat([y, y + period / 2]) % period)  # include antipodals
    y = np.unique(
        np.concat([
            y,
            (y[1:] + y[:-1]) / 2,
            [(y[-1] + y[0] + period) / 2 % period],
        ])
    )  # include mid points
    
    d = np.abs(y[:, np.newaxis] - x) % period
    d = np.where(d > period / 2, period - d, d)
    d = np.sum(d, axis=-1)

    dmin = min(d)
    if np.all(d == dmin):
        return [(0, period)]
    
    intervals = []
    i0 = np.argmax(d > dmin)
    interval = None
    m = len(y)
    for j in range(m):
        i = (i0 + j + 1) % len(d)
        if d[i] == dmin:
            if interval is None:
                interval = (float(y[i]), float(y[i]))
            else:
                interval = (interval[0], float(y[i]))
        elif interval is not None:
            intervals.append(interval)
            interval = None
    return intervals

def simulate(period: int, n: int):
    x = np.random.randint(0, period, size=n)
    med = arc_distance_median(x, period)
    assert len(med) >= 1
    if med == [(0, period)]:  # every point is a median
        print(f"{x=}, {med=}")
    elif n % 2 == 1 and not (len(med) == 1 and med[0][0] == med[0][1]):
        # odd sample; non-unique minimizer, or interval minimizer (should not happen)
        print(f"{x=}, {med=}")
    elif n % 2 == 0 and len(med) > 1:
        # even sample; non-unique interval
        print(f"{x=}, {med=}")
    else:
        pass  # nothing special to print

def main():
    np.random.seed(1234)
    for _ in range(100):
        simulate(period=12, n=5)
        # simulate(period=13, n=6)

if __name__ == "__main__":
    main()

x=array([0, 5, 1, 5, 9]), med=[(1.0, 1.0), (5.0, 5.0)]
x=array([10,  4,  2,  7,  9]), med=[(7.0, 7.0), (9.0, 9.0)]
x=array([ 7, 11,  7,  2,  9]), med=[(7.0, 7.0), (9.0, 9.0)]
x=array([0, 5, 4, 9, 2]), med=[(2.0, 2.0), (4.0, 4.0)]
x=array([9, 1, 3, 5, 9]), med=[(1.0, 1.0), (5.0, 5.0), (9.0, 9.0)]
x=array([6, 8, 6, 1, 9]), med=[(6.0, 6.0), (8.0, 8.0)]
x=array([ 8,  1,  9, 11,  4]), med=[(9.0, 9.0), (11.0, 11.0)]
x=array([6, 8, 9, 1, 5]), med=[(6.0, 6.0), (8.0, 8.0)]
x=array([ 5,  7, 10,  0,  0]), med=[(10.0, 10.0), (0.0, 0.0)]
x=array([ 6,  1,  3, 10,  5]), med=[(3.0, 3.0), (5.0, 5.0)]
x=array([ 2, 10,  6, 11,  6]), med=[(6.0, 6.0), (10.0, 10.0)]
x=array([5, 9, 1, 9, 2]), med=[(1.0, 1.0), (9.0, 9.0)]

library(circular)
median(circular(c(1, 2, 5, 9, 9) * 30, units = "degrees"))

Circular Data: 
Type = angles 
Units = degrees 
Template = none 
Modulo = asis 
Zero = 0 
Rotation = counter 
[1] -60
attr(,"medians")
[1]  30 270

Let me recap the discussion above and then propose a solution. The goal is to compute a sample circular median, an angle that minimizes the mean circular deviation. If the sample is large or drawn from a continuous distribution, it is almost certain that the sample circular median is unique. The question is what to do for a small sample drawn from a discrete distribution, which could be relevant e.g. if the sample data are binned.

When the minimizer for mean circular deviation is not unique, most (all?) statistical software and two academic papers choose to take the circular average of minimizers in the sample as the circular median. However, there are (non-contrived) examples where such an average do not minimize the mean circular deviation. Other academic literature note the possibility of multiple median, and one proposes to choose an arbitrary minimizer as the circular median.

Based on the above, I'd propose the following for the circmedian function:

• Specification: Return an angle that minimizes the mean circular deviation. If such angles are not unique, an arbitrary minimizer is returned.

• Implementation: If the sample size is odd and there is a unique minimizer, return that. If the sample size is even and there are exactly two minimizers in the sample and they are not antipodal, return their mid point. Otherwise, return the first minimizer in the sample.

• Post condition (accuracy guarantee): The mean circular deviation from the returned angle is equal to the (true) minimum circular deviation from any angle up to numerical noise (of a few epsilon).

Basically, we offer weak but sure guarantee in the specification, and make a best effort in the implementation. The weak guarantee may be strengthened later if there turns out to be a need to (which I doubt). This also eliminates the need to deal with epsilon in the implementation. (After all, the return value is extremely sensitive to numerical error when the circular mean is near zero, and it is therefore near impossible to guarantee the accuracy of the returned angle anyway. On the other hand, it's relatively easy to guarantee the optimality of the returned mean circular deviation.)

Are you happy with this PR @mdhaber ? Then I would take a detailed look.

Yup, I'm happy with it. I didn't see anything in the other comments that I wanted to follow up on. Thanks!

I've planned to take biteopt for a spin at some point. Can't review/merge, but I can add a +1.

To simplify review (and merge conflict resolution), I moved the existing circular/directional stats functions back to _morestats.py and kept circmedian in the new _circstats.py. The existing functions can be moved over in a separate PR.

GPU test failure is in main, so I pushed a lint-only commit.

Thanks @dschmitz89! At some point, this new input validation should probably be shared with the other circular statistics functions - maybe when they are moved into this _circstats.py file.

Thanks @dschmitz89!
There are a number of follow-ups to consider. Would you be interested in reviewing any?

• migrating the other circular statistics functions to _circstats.py (MAINT: stats: migrate circular statistics functions to _circstats.py #26006)
• sharing the input validation from this PR ( https://github.com/mdhaber/scipy/pull/new/circstats_iv)
• adding MArray support to circmedian
• adding the geometric median (top post) as a method (https://github.com/mdhaber/scipy/pull/new/circmedian_geometric)

Thought I might also try to interest you in Maritz-Jarrett standard errors (gh-24757).
In any case, thanks for taking a look at this!