实际上,我对这个话题的思考多于健康,所以我将在这里分享我的想法和发现。也许有人会遇到类似的问题并发现这很有用。
我已经很多年没用过C++了,如果我用C#写了所有的代码,请原谅我。我相信一个流利的 C++ 演讲者可以很容易地翻译算法。
循环平均值
首先,让我们定义circular mean。它是通过将您的点转换为弧度来计算的,其中您的周期(256、360 或其他任何值 - 被解释为与零相同的值)被缩放为 2*pi。然后计算这些弧度值的正弦和余弦。这些是单位圆上值的 y 和 x 坐标。然后将所有正弦和余弦相加并计算 atan2。这为您提供了平均角度,可以通过除以比例因子轻松转换回您的数据点。
var scalingFactor = 2 * Math.PI / period;
var sines = 0.0;
var cosines = 0.0;
foreach (var value in inputs)
{
var radians = value * scalingFactor;
sines += Math.Sin(radians);
cosines += Math.Cos(radians);
}
var circularMean = Math.Atan2(sines, cosines) / scalingFactor;
if (circularMean >= 0)
return circularMean;
else
return circularMean + period;
边际圆形中位数
最简单的圆形中位数方法只是处理圆形均值的一种修改方法。
可以以类似的方式计算圆形中位数,只需找到正弦和余弦的中位数而不是总和,然后计算其中的 atan2。这样,您就可以找到圆点的marginal median 并获取其角度。
var scalingFactor = 2 * Math.PI / period;
var sines = new List<double>();
var cosines = new List<double>();
foreach (var value in inputs)
{
var radians = value * scalingFactor;
sines.Add(Math.Sin(radians));
cosines.Add(Math.Cos(radians));
}
var circularMedian = Math.Atan2(Median(sines), Median(cosines)) / scalingFactor;
if (circularMedian >= 0)
return circularMedian;
else
return circularMedian + period;
这种方法是 O(n),对异常值具有鲁棒性并且实现起来非常简单。它可能非常适合您的目的,但它有一个问题:旋转输入点会给您带来不同的结果。根据输入数据的分布情况,这可能是也可能不是问题。
圆弧中线
要理解这种其他方法,您需要停止考虑“这是如何计算的”的均值和中位数,而是考虑结果值实际代表的含义。
对于非循环数据,您可以通过将所有值相加并除以元素数来获得平均值。然而,这个数字代表的是与数据元素的所有平方距离之和最小的值。 (我听说统计学家将此值称为位置的 L2 估计值,但统计学家可能应该确认或否认这一点。)
对于中位数也是如此。如果所有数据都已排序(理想情况下,使用 O(n) selection algorithm,如 C++ 中的 nth_element),您可以通过查找最终位于中间的数据元素来获得它。但是,这个数字是一个值,它具有到数据元素的所有绝对(非平方!)距离的最小总和。 (据说,这个值称为位置的 L1 估计值。)
对循环数据进行排序并不能帮助您找到中间点,因此考虑中位数的通常方法不起作用,但您仍然可以找到使与所有数据点的绝对距离总和最小化的点。这是我想出的算法,它在 O(n) 时间内运行,假设输入数据被标准化为 >= 0 和
它通过遍历所有数据点并跟踪距离总和来工作。当您向右数据点移动距离 D 时,到所有左侧点的距离总和增加 D*LeftCount,到所有右侧点的距离总和减少 D*RightCount。然后,如果某些左侧点现在实际上是右侧点,因为它们的左侧距离大于period/2,则减去它们之前的距离并添加新的正确距离。
为了将当前总和与最佳总和进行比较,我添加了一些容差以防止不精确的浮点运算。
可能有多个或无限多个满足最小距离条件的点。对于偶数个值的非圆形中位数,中位数可以是两个中心值之间的任何值。它通常被认为是这两个中心值的平均值,所以我对这个中值算法采取了类似的方法。我找到所有最小化距离的数据点,然后计算这些点的圆形平均值。
// Requires a sorted list with values normalized to [0,period).
// Doing an initialization pass:
// * candidate is the lowest number
// * finding the index where the circle with this candidate starts
// * calculating the score for this candidate - the sum of absolute distances
// * counting the number of values to the left of the candidate
int i;
var candidate = list[0];
var distanceSum = 0.0;
for (i = 1; i < list.Count; ++i)
{
if (list[i] >= candidate + period / 2)
break;
distanceSum += list[i] - candidate;
}
var leftCount = list.Count - i;
var circleStart = i;
if (circleStart == list.Count)
circleStart = 0;
else
for (; i < list.Count; ++i)
distanceSum += candidate + period - list[i];
var previousCandidate = candidate;
var bestCandidates = new List<double> { candidate };
var bestDistanceSum = distanceSum;
var equalityTolerance = period * 1e-10;
for (i = 1; i < list.Count; ++i)
{
candidate = list[i];
// A formula for correcting the distance given the movement to the right.
// It doesn't take into account that some values may have wrapped to the other side of the circle.
++leftCount;
distanceSum += (2 * leftCount - list.Count) * (candidate - previousCandidate);
// Counting all the values that wrapped to the other side of the circle
// and correcting the sum of distances from the candidate.
if (i <= circleStart)
while (list[circleStart] < candidate + period / 2)
{
--leftCount;
distanceSum += 2 * (list[circleStart] - candidate) - period;
++circleStart;
if (circleStart == list.Count)
{
circleStart = 0;
break; // Letting the next loop continue.
}
}
if (i > circleStart)
while (list[circleStart] < candidate - period / 2)
{
--leftCount;
distanceSum += 2 * (list[circleStart] - candidate) + period;
++circleStart;
}
// Comparing current sum to the best one, using the given tolerance.
if (distanceSum <= bestDistanceSum + equalityTolerance)
{
if (distanceSum >= bestDistanceSum - equalityTolerance)
{
// The numbers are close, so using their average as the next best.
bestDistanceSum = (bestCandidates.Count * bestDistanceSum + distanceSum) / (bestCandidates.Count + 1);
}
else
{
// The new number is significantly better, clearing.
bestDistanceSum = distanceSum;
bestCandidates.Clear();
}
bestCandidates.Add(candidate);
}
previousCandidate = candidate;
}
if (bestCandidates.Count == 1)
return bestCandidates[0];
else
return CircularMean(bestCandidates, period);
几何圆形中线
在之前的算法中存在一个不一致之处,即中位数相对于循环平均值的定义方式。圆形平均值最小化圆上点之间的平方欧几里得距离之和。换句话说,它查看连接圆上点的直线,穿过圆。
圆弧中位数,正如我在上面计算的那样,着眼于圆弧的距离:通过在圆的周长上移动,而不是在它们之间画一条直线,这些点之间的距离。
我已经考虑过如何解决这个问题,如果它困扰你,但我还没有真正做过任何实验,所以我不能声称以下方法有效。简而言之,我相信你可以使用Iteratively reweighted least squares algorithm (IRLS)的修改,这是通常用来计算geometric medians的。
这个想法是选择一个起始值(例如,上面显示的圆形平均值或弧形中值),并计算到每个点的欧几里得距离:Di = sqrt(dxi^2 + dyi^2)。圆形平均值将使这些距离的平方最小化,因此每个点的权重应该抵消平方并重置为 D:Wi = Di / Di^2,即 Wi = 1 / Di。
使用这些权重,计算加权循环平均值(与循环平均值相同,但在相加之前将每个正弦和余弦乘以该点的权重)并重复该过程。重复直到经过足够多的迭代或直到结果不再发生太大变化。
这个算法的问题是,如果当前解正好落在一个数据点上,它就会被零除。即使距离不完全为零,如果您足够接近该点,解决方案也会停止移动,因为与所有其他重量相比,重量会变得巨大。这可以通过在除以距离之前添加一个小的固定偏移量来解决。这将使解决方案变得次优,但至少不会停在错误的点上。
除非偏移量相对较大,否则仍然需要一些迭代才能将自己从错误点中挖掘出来,并且最终的解决方案越差,偏移量越大。所以最好的方法可能是从一个相当大的偏移量开始,然后在每次下一次迭代中逐渐减小它。