medians for degree measurements

[Robert Kern]

On 2010-01-22 18:09 PM, Steve Howell wrote:
> I just saw the thread for medians, and it reminded me of a problem
> that I need to solve.  We are writing some Python software for
> sailing, and we need to detect when we've departed from the median
> heading on the leg.  Calculating arithmetic medians is
> straightforward, but compass bearings add a twist.
>
> The numerical median of 1, 2, 3, 4, 5, 6, 359 is 4.  But for
> navigational purposes you would actually order the numbers 359, 1, 2,
> 3, 4, 5, 6, so the desired median heading of the boat is actually 3.
>
> Of course, you could express 359 better as -1 degrees to north, then
> the sequence would be -1, 1, 2, 3, 4, 5, and 6.  And you'd be good.
>
> But that trick does not generalize if you go south instead, as you
> have similar issues with -179, 174, 175, 176, 177, 178, and 179.
>
> Has anybody solved this in Python, either for compass bearings or a
> different domain?  I can think of kind of a brute force solution where
> you keep rotating the sequence until the endpoints are closest
> together mod 360, but I wonder if there is something more elegant.

Brute force doesn't suck too much when using numpy to do the heavy lifting.

In [158]: import numpy as np

# Worst case. A von Mises distribution "centered" around the "South pole"
# at pi/-pi.
In [159]: theta = np.random.vonmises(np.pi, 1.0, size=1000)

# Complex division makes circular distances easier to compute.
In [160]: z = np.exp(1j * theta)

# The newaxis bit plus numpy's broadcasting yields a 1000x1000 array of
# the quotient of each pair of points in the dataset.
In [161]: zdiv = z / z[:,np.newaxis]

# Convert back to angular distances. Take the absolute value. Sum along one
# dimension. Find the index of the element with the smallest mean absolute
# circular deviation when used as the putative median.
In [162]: i = abs(np.arctan2(zdiv.imag, zdiv.real)).sum(axis=1).argmin()

# As expected, the result is close to the mode of the distribution.
In [163]: theta[i]
Out[163]: 3.0886753423606521

-- 
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma
  that is made terrible by our own mad attempt to interpret it as though it had
  an underlying truth."
   -- Umberto Eco