【问题标题】:Finding median of a set of circular data查找一组循环数据的中位数
【发布时间】:2018-03-09 03:38:06
【问题描述】:

我想编写一个 C++ 函数来查找循环数据数组的中值。 例如,考虑罗盘的读数,假设读数在 [0,360) 范围内。虽然 1 和 359 看起来很远,但由于读数的循环性质,它们非常接近。

普通数据中N个元素的求中位数如下。 1.对N个元素的数据进行排序(升序或降序) 2. 如果 N 为奇数,则中位数为已排序数组中的第 (N+1)/2 个元素。 3. 如果 N 是偶数,则中位数是排序数组中第 N/2 和 N/2+1 个元素的平均值。

但是,循环数据中的环绕问题将问题带到了不同的维度,并且解决方案不平凡。

此处解释了从循环数据中查找均值的类似问题How do you calculate the average of a set of circular data? 上面链接中的建议是找到每个角度对应的单位向量,然后求平均值。但是,中位数需要对数据进行排序,而向量的排序在这种情况下没有任何意义。因此,我认为我们不能使用建议的方案来找到中位数!

【问题讨论】:

  • 其实我认为中位数的概念并没有自然地延伸到这个案例中。我认为您需要一个额外的条件,例如:中位数,使其周围的分布最小
  • 考虑例如[ 0 180],其中90270 一样好...或[0 60 120 180 270] 有5 种可能的解决方案,甚至“最小传播”条件也无济于事唯一地选择其中之一
  • 您能否给出循环值中位数的定义并解释其含义?
  • 也许你会在math.stackexchange.com获得更多运气
  • 但是360和0没有区别,我觉得@tobi303是对的。

标签: c++ algorithm median


【解决方案1】:

使用您的角度数据点向量(即从 0 到 259 的数字向量),创建两个新向量,我将它们称为 xy。这两个新向量分别是角度数据点的正弦和余弦。

也就是说,x[n] = cos(data[n])y[n] = sin(data[n]) 其中data 是您的角度数据向量,n 是有很多数据点。

接下来,将x 向量中的所有值相加得到一个值,称之为sum_x,并将y 向量中的所有值相加得到另一个值,称之为@ 987654330@.

现在你可以做正切逆(例如atan(sum_y/sum_x))来获得一个新值。而这个值是非常有意义的。该值基本上告诉您数据“指向”的方向,即大部分数据存在的位置。注意:您必须小心除以 0(当 sum_x=0 时)和出现不确定形式时(当 sum_x=0 和 sum_y=0 时)。不确定的形式只是意味着你的数据是均匀分布的,在这种情况下中位数是没有意义的,当sum_x=0 但sum_y!=0 时,它实际上是atan(inf)atan(-inf),两者都是已知。

编辑:

在此之后,我之前的答案需要进行一些调整。

从这里开始,很容易。取您在上一步中获得的值 (atan(sum_y/sum_x)) 并将该值加上 180 度。这是您的数据开始和结束的参考点。从这里,您可以将此参考点作为起点和终点对角度数据进行排序,并找到该数据的中位数。

【讨论】:

  • 巧妙但遗憾的是不正确。您的 arctan(sum_y / sum_x) 是平均点,您断言中值是最接近该点的值没有统计依据。
  • 我明白你的意思,但我不完全确定这是不是真的。在循环数据中,中位数和平均值的关系更为密切。
  • 对我来说,以与平均值相反的点为原点,旋转变换你的数据,使原点的值为零,取该集合的中值,然后再变换回来是一个很好的定义。但我没有引用这种方法。
  • 但我认为现在这里有足够的材料可以投赞成票。来了。
  • 这是个好主意。我现在修改我的答案。
【解决方案2】:

中位数的两个属性允许发明两种不同的算法来寻找中位数。

1) 中位数最小化到所有其他元素的绝对距离之和 -- O(n^2) 算法:

for (i = 0; i < N; i++)
{
     sum = 0;
     for (j = 0; j < N; j++)
        sum += abs(item[i] - item[j]) % 360;
     if (sum < best_so_far) { best_so_far = sum; index = i; }
}

2) 中位数满足一半的项目少,一半的项目大

  • 对项目进行排序
  • 找到第一组项目 (i=0...I),满足 I i + 180
  • 如果不满足中位数条件,则前进 i 或 I。
  • 排序需要 O(N*log N),下一次扫描需要 O(N)

当然,在周期性数据中,所有项目(以及数据点之间的所有项目)都可以成为中位数的合适候选者。

【讨论】:

    【解决方案3】:

    不可能将中位数的概念规范地扩展到循环数据。为简单起见,让我们考虑[0 10) 中的数字,并以(已排序的)集合{ 1 3 5 7 8 } 为例。根据您旋转数组的方式,您会获得不同的中位数值:

    1 3 5 7 8    -> 5
    3 5 7 8 1    -> 7
    5 7 8 1 3    -> 8
    ...etc...
    

    任何一个都和另一个一样好。

    声称不可能在循环数据上定义中位数。我只是声称“正常”中位数不能以有意义的方式扩展到这种情况,而不添加额外的约束或做出任意选择。

    【讨论】:

      【解决方案4】:

      有关圆形中位数的定义和讨论,请参阅

      N.I.费舍尔的“循环数据的统计分析”,剑桥大学。 1993年出版

      以及围绕方程 2.32 和 2.33 的讨论。对于多模态或各向同性数据,可能不存在唯一的中位数。

      找到一个将数据分成 2 个相等组的轴,并选择该轴在角度较小值处的末端。如果样本量是奇数,则中位数将是一个数据点,否则它将是 2 个数据点的中点。

      有其他语言的包(例如 R、MatLab)可以帮助为您编写的任何函数提供测试值。

      例如 https://www.rdocumentation.org/packages/circular/versions/0.4-93

      具体见median.circularmedianHL.circular

      贝伦斯,菲利普。 “CircStat:循环统计的 MATLAB 工具箱”。统计软件杂志 31,没有。 1(2009 年 9 月 23 日):1-21。 https://doi.org/10.18637/jss.v031.i10.

      看看circ_median

      【讨论】:

        【解决方案5】:

        实际上,我对这个话题的思考多于健康,所以我将在这里分享我的想法和发现。也许有人会遇到类似的问题并发现这很有用。

        我已经很多年没用过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。

        使用这些权重,计算加权循环平均值(与循环平均值相同,但在相加之前将每个正弦和余弦乘以该点的权重)并重复该过程。重复直到经过足够多的迭代或直到结果不再发生太大变化。

        这个算法的问题是,如果当前解正好落在一个数据点上,它就会被零除。即使距离不完全为零,如果您足够接近该点,解决方案也会停止移动,因为与所有其他重量相比,重量会变得巨大。这可以通过在除以距离之前添加一个小的固定偏移量来解决。这将使解决方案变得次优,但至少不会停在错误的点上。

        除非偏移量相对较大,否则仍然需要一些迭代才能将自己从错误点中挖掘出来,并且最终的解决方案越差,偏移量越大。所以最好的方法可能是从一个相当大的偏移量开始,然后在每次下一次迭代中逐渐减小它。

        【讨论】:

        • 在 OP 提到的平均主题中,有一种算法可能与我提出的弧中位数方法很好地配对,用于多个候选人的最终平均:stackoverflow.com/a/3651941/4247453 与循环平均值不同方法,该算法通过弧长而不是直线距离来选择平均值。但我还没有真正研究过这个算法,我不确定它是否能保证从弧中位数给出的候选中产生一个唯一的点。
        猜你喜欢
        • 2022-12-05
        • 1970-01-01
        • 2014-01-12
        • 1970-01-01
        • 1970-01-01
        • 1970-01-01
        • 2016-08-04
        • 1970-01-01
        • 1970-01-01
        相关资源
        最近更新 更多