medians for degree measurements
Robert Kern
robert.kern at gmail.com
Sat Jan 23 01:45:46 CET 2010
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
More information about the Python-list
mailing list