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
--
http://mail.python.org/mailman/listinfo/python-list
