如何计算一组循环数据的平均值?

我想计算一组循环数据的平均值。例如,我可能有几个样本从读取指南针。问题当然是如何处理包装。同样的算法可能对钟面有用。

实际的问题更加复杂——统计学在一个球体上或者在一个“环绕”的代数空间中意味着什么,例如加法群 mod n。答案可能不是唯一的,例如359度和1度的平均值可能是0度或180度,但从统计学上来说0看起来更好。

这对我来说是一个真正的编程问题,我试图让它看起来不仅仅是一个数学问题。

108900 次浏览

从角度计算单位向量,取其平均角度。

您必须更准确地定义 一般。对于两个角度的具体情况,我可以想到两种不同的情况:

  1. 「真实」平均数,即(a + b)/2% 360。
  2. 在同一个半圆内指向另外两个半圆之间的角度,例如,对于355和5,这个角度是0,而不是180。为此,您需要检查两个角度之间的差值是否大于180。如果是这样,在使用上述公式之前,将较小的角度增加360。

不过,我不明白如何将第二种选择推广到两个以上角度的情况。

我看到了这个问题-例如,如果你有一个45’角和一个315’角,“自然”平均值将是180’,但是你想要的值实际上是0’。

我觉得星蓝有线索了。只需计算每个角的(x,y)笛卡尔坐标,然后把得到的向量相加。最终矢量的角度偏移量应该是你所需要的结果。

x = y = 0
foreach angle {
x += cos(angle)
y += sin(angle)
}
average_angle = atan2(y, x)

我现在忽略了一个罗盘的方向是从北边开始,然后是顺时针方向,而“正常”的笛卡尔坐标是从 X 轴的零开始,然后是逆时针方向。不管怎样,数学都应该用同样的方法来计算。

这里有一个想法: 通过总是计算最接近的角度的平均值来迭代地构建平均值,同时保持一个权重。

另一个想法是: 找出给定角度之间的最大间隙。找到平分它的点,然后选择圆上的相反点作为参考零点,从中计算平均值。

让我们用圆周上的点来表示这些角。

我们能假设所有这些点都落在圆的同一半上吗?(否则,没有明显的方法来定义“平均角”。想想直径上的两个点,例如0度和180度-,平均温度是90度还是270度?当我们有3个或更多平均分布的点时会发生什么?)

根据这个假设,我们选择那个半圆上的任意一个点作为“原点”,并测量相对于这个原点的给定角度集合(称之为“相对角度”)。请注意,相对角度的绝对值严格小于180度。最后,取这些相对角度的平均值,得到期望的平均角度(当然是相对于我们的原点)。

这个问题在书中有详细的论述: “球体统计学”,杰弗里 · 沃森,阿肯色大学讲座 《数学科学笔记》 ,1983年约翰威立,公司,如布鲁斯 · 卡什在 http://catless.ncl.ac.uk/Risks/7.44.html#subj4上提到的。

从一组角度测量值估计平均角度的好方法 A [ i ]0 < = i

                   sum_i_from_1_to_N sin(a[i])
a = arctangent ---------------------------
sum_i_from_1_to_N cos(a[i])

Starblue 给出的方法在计算上是等价的,但是他的理由更清楚,可能在编程上更有效,而且在0情况下也工作得很好,所以值得称赞。

这个主题现在被探讨在更详细的 在维基百科上,并与其他用途,如小数部分。

Alnitak 有正确的解决方案 Nick Fortescue 的解决方案在功能上是一样的。

特殊情况下

(sum (x 分量) = 0.0 & & sum (y 分量) = 0.0)//例如2个10度和190度的角。

用0.0度作为总和

在计算上,您必须测试这种情况,因为 atan2(0. ,0.)是未定义的,并将生成一个错误。

平均角 phi _ avg 应该具有 sum _ i | phi _ avg-phi _ i | ^ 2变为最小的性质,其中差值必须是[-Pi,Pi)(因为反过来可能会更短!).通过将所有输入值标准化为[0,2Pi) ,保持运行平均 phi _ run,并选择将 | phi _ i-phi _ run | 标准化为[-Pi,Pi) ,可以很容易地实现这一点 (通过增加或减去2Pi)。上面的大多数建议做一些其他的事情,做 没有 有最小的属性,也就是说,他们平均 什么东西,但不是角度。

对于两个角度的特殊情况:

答案 ((a + b)模数360)/2。对于角度350和2,最接近的点是356,而不是176。

单位向量和三角解决方案可能太昂贵。

我从一个小小的修补中得到的是:

diff = ( ( a - b + 180 + 360 ) mod 360 ) - 180
angle = (360 + b + ( diff / 2 ) ) mod 360
  • 0,180-> 90(这个方程有两个答案: 这个方程从 a 取顺时针方向的答案)
  • 180,0-> 270(见上文)
  • 180,1-> 90.5
  • 1,180-> 90.5
  • 20,350-> 5
  • 350,20-> 5(下面的例子也适当地反过来)
  • 10,20-> 15
  • 350,2-> 356
  • 359,0-> 359.5
  • 180,180-> 180

我会使用复数的矢量方式,我的例子是 Python,它有内置的复数:

import cmath # complex math


def average_angle(list_of_angles):


# make a new list of vectors
vectors= [cmath.rect(1, angle) # length 1 for each vector
for angle in list_of_angles]


vector_sum= sum(vectors)


# no need to average, we don't care for the modulus
return cmath.phase(vector_sum)

请注意,Python 不会用 需要来构建一个临时的新向量列表,以上所有操作都可以在一个步骤中完成; 我只是选择这种方法来近似适用于其他语言的伪代码。

Ackb 是正确的,这些基于矢量的解决方案不能被认为是真正的平均角度,他们只是一个平均的单位矢量的对应物。然而,ackb 建议的解决方案在数学上似乎并不合理。

下面是一个从最小化(角[ i ]-avgAngle) ^ 2(如果需要,差异会被修正)的目标中得到的数学解,这使得它成为一个真正的角度的算术平均值。

首先,我们需要准确地观察在哪些情况下,角度之间的差异不同于它们的法数对应物之间的差异。考虑角 x 和 y,如果 y > = x-180和 y < = x + 180,那么我们可以直接使用差(x-y)。否则,如果不满足第一个条件,那么我们必须在计算中使用(y + 360)而不是 y。相应地,如果不满足第二个条件,那么我们必须使用(y-360)而不是 y。由于曲线方程我们只是在这些不等式从真到假或者反过来变化的点上最小化变化,我们可以将完整的[0,360]范围分割成一组段,由这些点分隔。然后,我们只需要找到每个片段的最小值,然后找到每个片段的最小值,也就是平均值。

下面的图片演示了在计算角度差时出现的问题。如果 x 位于灰色区域,那么就会有问题。

Angle comparisons

为了使一个变量最小化,根据曲线,我们可以对我们想要使其最小化的东西求导,然后我们找到转折点(这就是导数 = 0的地方)。

在这里,我们将应用最小平方差的思想来推导一般的算术平均公式: sum (a [ i ])/n。曲线 y = sum ((a [ i ]-x) ^ 2)可以这样最小化:

y = sum((a[i]-x)^2)
= sum(a[i]^2 - 2*a[i]*x + x^2)
= sum(a[i]^2) - 2*x*sum(a[i]) + n*x^2


dy\dx = -2*sum(a[i]) + 2*n*x


for dy/dx = 0:
-2*sum(a[i]) + 2*n*x = 0
-> n*x = sum(a[i])
-> x = sum(a[i])/n

现在把它应用到我们调整后的差异曲线上:

B = a 的子集,其中正确的(角)差 a [ i ]-x C = 其中正确(角)差(a [ i ]-360)-x 的子集 Cn = c 的大小 D = 其中正确(角)差(a [ i ] + 360)-x 的子集 Dn = d 的大小

y = sum((b[i]-x)^2) + sum(((c[i]-360)-b)^2) + sum(((d[i]+360)-c)^2)
= sum(b[i]^2 - 2*b[i]*x + x^2)
+ sum((c[i]-360)^2 - 2*(c[i]-360)*x + x^2)
+ sum((d[i]+360)^2 - 2*(d[i]+360)*x + x^2)
= sum(b[i]^2) - 2*x*sum(b[i])
+ sum((c[i]-360)^2) - 2*x*(sum(c[i]) - 360*cn)
+ sum((d[i]+360)^2) - 2*x*(sum(d[i]) + 360*dn)
+ n*x^2
= sum(b[i]^2) + sum((c[i]-360)^2) + sum((d[i]+360)^2)
- 2*x*(sum(b[i]) + sum(c[i]) + sum(d[i]))
- 2*x*(360*dn - 360*cn)
+ n*x^2
= sum(b[i]^2) + sum((c[i]-360)^2) + sum((d[i]+360)^2)
- 2*x*sum(x[i])
- 2*x*360*(dn - cn)
+ n*x^2


dy/dx = 2*n*x - 2*sum(x[i]) - 2*360*(dn - cn)


for dy/dx = 0:
2*n*x - 2*sum(x[i]) - 2*360*(dn - cn) = 0
n*x = sum(x[i]) + 360*(dn - cn)
x = (sum(x[i]) + 360*(dn - cn))/n

仅仅这样还不足以得到最小值,因为它适用于具有无界集合的正常值,所以结果肯定在集合的范围内,因此是有效的。我们需要一个范围内的最小值(由段定义)。如果最小值小于我们段的下界,那么段的最小值必须在下界(因为二次曲线只有一个转折点) ,如果最小值大于我们段的上界,那么段的最小值就在上界。在我们得到每个片段的最小值之后,我们只需找到一个最小值(sum ((b [ i ]-x) ^ 2) + sum (((c [ i ]-360)-b) ^ 2) + sum (((d [ i ] + 360)-c) ^ 2))。

下面是曲线的图像,它显示了在 x = (a [ i ] + 180)% 360的点处它是如何变化的。有问题的数据集是{65,92,230,320,250}。

Curve

下面是一个在 Java 中实现的算法,包括一些优化,其复杂度为 O (nlogn)。如果将基于比较的排序替换为基于非比较的排序(如基数排序) ,则可以将其减少为 O (n)。

static double varnc(double _mean, int _n, double _sumX, double _sumSqrX)
{
return _mean*(_n*_mean - 2*_sumX) + _sumSqrX;
}
//with lower correction
static double varlc(double _mean, int _n, double _sumX, double _sumSqrX, int _nc, double _sumC)
{
return _mean*(_n*_mean - 2*_sumX) + _sumSqrX
+ 2*360*_sumC + _nc*(-2*360*_mean + 360*360);
}
//with upper correction
static double varuc(double _mean, int _n, double _sumX, double _sumSqrX, int _nc, double _sumC)
{
return _mean*(_n*_mean - 2*_sumX) + _sumSqrX
- 2*360*_sumC + _nc*(2*360*_mean + 360*360);
}


static double[] averageAngles(double[] _angles)
{
double sumAngles;
double sumSqrAngles;


double[] lowerAngles;
double[] upperAngles;


{
List<Double> lowerAngles_ = new LinkedList<Double>();
List<Double> upperAngles_ = new LinkedList<Double>();


sumAngles = 0;
sumSqrAngles = 0;
for(double angle : _angles)
{
sumAngles += angle;
sumSqrAngles += angle*angle;
if(angle < 180)
lowerAngles_.add(angle);
else if(angle > 180)
upperAngles_.add(angle);
}




Collections.sort(lowerAngles_);
Collections.sort(upperAngles_,Collections.reverseOrder());




lowerAngles = new double[lowerAngles_.size()];
Iterator<Double> lowerAnglesIter = lowerAngles_.iterator();
for(int i = 0; i < lowerAngles_.size(); i++)
lowerAngles[i] = lowerAnglesIter.next();


upperAngles = new double[upperAngles_.size()];
Iterator<Double> upperAnglesIter = upperAngles_.iterator();
for(int i = 0; i < upperAngles_.size(); i++)
upperAngles[i] = upperAnglesIter.next();
}


List<Double> averageAngles = new LinkedList<Double>();
averageAngles.add(180d);
double variance = varnc(180,_angles.length,sumAngles,sumSqrAngles);


double lowerBound = 180;
double sumLC = 0;
for(int i = 0; i < lowerAngles.length; i++)
{
//get average for a segment based on minimum
double testAverageAngle = (sumAngles + 360*i)/_angles.length;
//minimum is outside segment range (therefore not directly relevant)
//since it is greater than lowerAngles[i], the minimum for the segment
//must lie on the boundary lowerAngles[i]
if(testAverageAngle > lowerAngles[i]+180)
testAverageAngle = lowerAngles[i];


if(testAverageAngle > lowerBound)
{
double testVariance = varlc(testAverageAngle,_angles.length,sumAngles,sumSqrAngles,i,sumLC);


if(testVariance < variance)
{
averageAngles.clear();
averageAngles.add(testAverageAngle);
variance = testVariance;
}
else if(testVariance == variance)
averageAngles.add(testAverageAngle);
}


lowerBound = lowerAngles[i];
sumLC += lowerAngles[i];
}
//Test last segment
{
//get average for a segment based on minimum
double testAverageAngle = (sumAngles + 360*lowerAngles.length)/_angles.length;
//minimum is inside segment range
//we will test average 0 (360) later
if(testAverageAngle < 360 && testAverageAngle > lowerBound)
{
double testVariance = varlc(testAverageAngle,_angles.length,sumAngles,sumSqrAngles,lowerAngles.length,sumLC);


if(testVariance < variance)
{
averageAngles.clear();
averageAngles.add(testAverageAngle);
variance = testVariance;
}
else if(testVariance == variance)
averageAngles.add(testAverageAngle);
}
}




double upperBound = 180;
double sumUC = 0;
for(int i = 0; i < upperAngles.length; i++)
{
//get average for a segment based on minimum
double testAverageAngle = (sumAngles - 360*i)/_angles.length;
//minimum is outside segment range (therefore not directly relevant)
//since it is greater than lowerAngles[i], the minimum for the segment
//must lie on the boundary lowerAngles[i]
if(testAverageAngle < upperAngles[i]-180)
testAverageAngle = upperAngles[i];


if(testAverageAngle < upperBound)
{
double testVariance = varuc(testAverageAngle,_angles.length,sumAngles,sumSqrAngles,i,sumUC);


if(testVariance < variance)
{
averageAngles.clear();
averageAngles.add(testAverageAngle);
variance = testVariance;
}
else if(testVariance == variance)
averageAngles.add(testAverageAngle);
}


upperBound = upperAngles[i];
sumUC += upperBound;
}
//Test last segment
{
//get average for a segment based on minimum
double testAverageAngle = (sumAngles - 360*upperAngles.length)/_angles.length;
//minimum is inside segment range
//we test average 0 (360) now
if(testAverageAngle < 0)
testAverageAngle = 0;


if(testAverageAngle < upperBound)
{
double testVariance = varuc(testAverageAngle,_angles.length,sumAngles,sumSqrAngles,upperAngles.length,sumUC);


if(testVariance < variance)
{
averageAngles.clear();
averageAngles.add(testAverageAngle);
variance = testVariance;
}
else if(testVariance == variance)
averageAngles.add(testAverageAngle);
}
}




double[] averageAngles_ = new double[averageAngles.size()];
Iterator<Double> averageAnglesIter = averageAngles.iterator();
for(int i = 0; i < averageAngles_.length; i++)
averageAngles_[i] = averageAnglesIter.next();




return averageAngles_;
}

一组角的算术平均值可能不符合你对平均值应该是多少的直观想法。例如,集合{179,179,0,181,181}的算术平均值是216(和144)。你马上想到的答案可能是180,然而众所周知,算术平均值受边缘值的影响很大。你还应该记住,角度不是向量,有时处理角度时,它看起来很吸引人。

这个算法当然也适用于所有服从同余关系的数量(只需要最小的调整) ,比如一天中的时间。

我还想强调的是,即使这是一个真正的角度平均值,不像矢量解,这并不一定意味着它是你应该使用的解,相应的单位向量的平均值很可能是你实际应该使用的值。

像所有的平均值一样,答案取决于度量的选择。对于给定的度量 M,[1,N ]中 k 在[-pi,pi ]中的某些角 a _ k 的平均值是 a _ M,它使 d ^ 2 _ M (a _ M,a _ k)的平方距离之和最小。对于加权平均数,只需在和中包含加权数 w _ k (使 sum _ kw _ k = 1)。就是,

A _ M = arg min _ x sum _ k w _ k d ^ 2 _ M (x,a _ k)

两种常见的度量选择是 Frobenius 度量和 Riemann 度量。对于 Frobenius 度量,存在一个直接的公式,它对应于循环统计中通常的平均方位概念。详情请参阅《旋转组中的平均值和平均值》 ,Maher Moakher,SIAM 矩阵分析与应用杂志,第24卷,第1期,2002年。
Http://link.aip.org/link/?sjmael/24/1/1

下面是 GNU Octave 3.2.4的一个函数,它执行计算:

function ma=meanangleoct(a,w,hp,ntype)
%   ma=meanangleoct(a,w,hp,ntype) returns the average of angles a
%   given weights w and half-period hp using norm type ntype
%   Ref: "Means and Averaging in the Group of Rotations",
%   Maher Moakher, SIAM Journal on Matrix Analysis and Applications,
%   Volume 24, Issue 1, 2002.


if (nargin<1) | (nargin>4), help meanangleoct, return, end
if isempty(a), error('no measurement angles'), end
la=length(a); sa=size(a);
if prod(sa)~=la, error('a must be a vector'); end
if (nargin<4) || isempty(ntype), ntype='F'; end
if ~sum(ntype==['F' 'R']), error('ntype must be F or R'), end
if (nargin<3) || isempty(hp), hp=pi; end
if (nargin<2) || isempty(w), w=1/la+0*a; end
lw=length(w); sw=size(w);
if prod(sw)~=lw, error('w must be a vector'); end
if lw~=la, error('length of w must equal length of a'), end
if sum(w)~=1, warning('resumming weights to unity'), w=w/sum(w); end


a=a(:);     % make column vector
w=w(:);     % make column vector
a=mod(a+hp,2*hp)-hp;    % reduce to central period
a=a/hp*pi;              % scale to half period pi
z=exp(i*a); % U(1) elements


% % NOTA BENE:
% % fminbnd can get hung up near the boundaries.
% % If that happens, shift the input angles a
% % forward by one half period, then shift the
% % resulting mean ma back by one half period.
% X=fminbnd(@meritfcn,-pi,pi,[],z,w,ntype);


% % seems to work better
x0=imag(log(sum(w.*z)));
X=fminbnd(@meritfcn,x0-pi,x0+pi,[],z,w,ntype);


% X=real(X);              % truncate some roundoff
X=mod(X+pi,2*pi)-pi;    % reduce to central period
ma=X*hp/pi;             % scale to half period hp


return
%%%%%%


function d2=meritfcn(x,z,w,ntype)
x=exp(i*x);
if ntype=='F'
y=x-z;
else % ntype=='R'
y=log(x'*z);
end
d2=y'*diag(w)*y;
return
%%%%%%


% %   test script
% %
% % NOTA BENE: meanangleoct(a,[],[],'R') will equal mean(a)
% % when all abs(a-b) < pi/2 for some value b
% %
% na=3, a=sort(mod(randn(1,na)+1,2)-1)*pi;
% da=diff([a a(1)+2*pi]); [mda,ndx]=min(da);
% a=circshift(a,[0 2-ndx])    % so that diff(a(2:3)) is smallest
% A=exp(i*a), B1=expm(a(1)*[0 -1; 1 0]),
% B2=expm(a(2)*[0 -1; 1 0]), B3=expm(a(3)*[0 -1; 1 0]),
% masimpl=[angle(mean(exp(i*a))) mean(a)]
% Bsum=B1+B2+B3; BmeanF=Bsum/sqrt(det(Bsum));
% % this expression for BmeanR should be correct for ordering of a above
% BmeanR=B1*(B1'*B2*(B2'*B3)^(1/2))^(2/3);
% mamtrx=real([[0 1]*logm(BmeanF)*[1 0]' [0 1]*logm(BmeanR)*[1 0]'])
% manorm=[meanangleoct(a,[],[],'F') meanangleoct(a,[],[],'R')]
% polar(a,1+0*a,'b*'), axis square, hold on
% polar(manorm(1),1,'rs'), polar(manorm(2),1,'gd'), hold off


%     Meanangleoct Version 1.0
%     Copyright (C) 2011 Alphawave Research, robjohnson@alphawaveresearch.com
%     Released under GNU GPLv3 -- see file COPYING for more info.
%
%     Meanangle is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or (at
%     your option) any later version.
%
%     Meanangle is distributed in the hope that it will be useful, but
%     WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
%     General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see `http://www.gnu.org/licenses/'.

我有一个不同于@Starblue 的方法,它可以给出上面给出的一些角度的“正确”答案,例如:

  • Angle _ avg ([350,10]) = 0
  • Angle _ avg ([-90,90,40]) = 13.333
  • Angle _ avg ([350,2]) = 356

它使用一个和来表示连续角度之间的差异。 代码(在 Matlab 中) :

function [avg] = angle_avg(angles)
last = angles(1);
sum = angles(1);
for i=2:length(angles)
diff = mod(angles(i)-angles(i-1)+ 180,360)-180
last = last + diff;
sum = sum + last;
end
avg = mod(sum/length(angles), 360);
end

没有单一的“正确答案”。我建议读这本书, 马迪亚和体育朱普,“方向性统计”,(威利,1999年) , 进行全面分析。

我想分享一个方法,我使用的微控制器没有浮点数或三角函数的能力。我仍然需要“平均”10原始轴承读数,以便平滑的变化。

  1. 检查第一个轴承的范围是270-360度还是0-90度(北两象限)
  2. 如果是的话,将该值和所有后续读数旋转180度,将所有值保持在0 < = 方位 < 360的范围内。否则就按照读数来。
  3. 一旦取得10个读数,假设没有包罗万象的情况,就可以计算数值平均值
  4. 如果180度旋转已经生效,然后旋转计算平均180度回到一个“真正的”轴承。

不是很理想,会坏的。在这种情况下,我侥幸逃脱了,因为这个装置只能非常缓慢地旋转。我会把它放在那里,以防其他人发现自己在类似的限制下工作。

英语:

  1. 制作第二个数据集,所有角度移动180。
  2. 取两个数据集的方差。
  3. 取方差最小的数据集的平均值。
  4. 如果这个平均值来自移位集,那么将答案再移180。

在巨蟒中:

一个 # numpy NX1角度数组

if np.var(A) < np.var((A-180)%360):
average = np.average(A)


else:
average = (np.average((A-180)%360)+180)%360

(只是想从参数估测或推论统计学上分享我的观点)

Nimble 的试验是得到一组角度的最小均方误差(MMSE)估计值,但它是找到一个“平均”方向的选择之一,人们也可以找到一个最小均方误差(MMAE)估计值,或者其他一些估计值作为“平均”方向,它取决于你对方向的量化误差,或者更一般地,在参数估测中,成本函数的定义。

^ MMSE/MMAE 对应于最小均方误差/绝对误差。

Ackb 说: “平均角 phi _ avg 应该具有 sum _ i | phi _ avg-phi _ i | ^ 2变为最小的性质... 它们平均一些东西,但不是角度。”

--你用均方法量化错误,这是最常见的方法之一,但不是唯一的方法。这里大多数人喜欢的答案(即单位向量之和和得到结果的角度)实际上是合理的解决方案之一。如果向量的方向被建模为冯 · 米塞斯分布,那么(可以证明)最大似然估计就是我们想要的“平均”方向。这种分布并不花哨,只是从二维高斯分布的周期性采样分布。看看伊恩。毕晓普的著作《模式识别与机器学习》(2.179)。同样,它绝不是代表“平均”方向的唯一最好的方法,然而,它是相当合理的,既有很好的理论依据又有简单的实现。

Nimble 说“ ackb 是正确的,这些基于矢量的解决方案不能被认为是角度的真平均值,它们只是单位矢量的平均值。”

——这不是真的。“单位向量对应物”揭示了向量方向的信息。角度是一个没有考虑向量长度的量,单位向量是一个带有附加信息的量,长度是1。你可以定义长度为2的“单位”向量,这并不重要。

以下是完整的解决方案: (输入为方位数组(0-360)

public static int getAvarageBearing(int[] arr)
{
double sunSin = 0;
double sunCos = 0;
int counter = 0;


for (double bearing : arr)
{
bearing *= Math.PI/180;


sunSin += Math.sin(bearing);
sunCos += Math.cos(bearing);
counter++;
}


int avBearing = INVALID_ANGLE_VALUE;
if (counter > 0)
{
double bearingInRad = Math.atan2(sunSin/counter, sunCos/counter);
avBearing = (int) (bearingInRad*180f/Math.PI);
if (avBearing<0)
avBearing += 360;
}


return avBearing;
}

下面是一个完整的 C + + 解决方案:

#include <vector>
#include <cmath>


double dAngleAvg(const vector<double>& angles) {
auto avgSin = double{ 0.0 };
auto avgCos = double{ 0.0 };
static const auto conv      = double{ 0.01745329251994 }; // PI / 180
static const auto i_conv    = double{ 57.2957795130823 }; // 180 / PI
for (const auto& theta : angles) {
avgSin += sin(theta*conv);
avgCos += cos(theta*conv);
}
avgSin /= (double)angles.size();
avgCos /= (double)angles.size();
auto ret = double{ 90.0 - atan2(avgCos, avgSin) * i_conv };
if (ret<0.0) ret += 360.0;
return fmod(ret, 360.0);
}

它以双精度矢量的形式取角,并将平均值简单地作为双精度矢量返回。角度必须是以度为单位的,当然平均值也是以度为单位的。

我在@David _ Hanak 的帮助下解决了这个问题。 正如他所说:

在同一个半圆内指向另外两个半圆之间的角度,例如,对于355和5,这个角度是0,而不是180。为此,您需要检查两个角度之间的差值是否大于180。如果是这样,在使用上述公式之前,将较小的角度增加360。

所以我计算了所有角度的平均值。然后所有小于这个的角度,增加360。然后重新计算平均值,把它们全部加起来,除以它们的长度。

        float angleY = 0f;
int count = eulerAngles.Count;


for (byte i = 0; i < count; i++)
angleY += eulerAngles[i].y;


float averageAngle = angleY / count;


angleY = 0f;
for (byte i = 0; i < count; i++)
{
float angle = eulerAngles[i].y;
if (angle < averageAngle)
angle += 360f;
angleY += angle;
}


angleY = angleY / count;

效果很好。

Python 函数:

from math import sin,cos,atan2,pi
import numpy as np
def meanangle(angles,weights=0,setting='degrees'):
'''computes the mean angle'''
if weights==0:
weights=np.ones(len(angles))
sumsin=0
sumcos=0
if setting=='degrees':
angles=np.array(angles)*pi/180
for i in range(len(angles)):
sumsin+=weights[i]/sum(weights)*sin(angles[i])
sumcos+=weights[i]/sum(weights)*cos(angles[i])
average=atan2(sumsin,sumcos)
if setting=='degrees':
average=average*180/pi
return average

你可以在 Matlab 中使用这个函数:

function retVal=DegreeAngleMean(x)


len=length(x);


sum1=0;
sum2=0;


count1=0;
count2=0;


for i=1:len
if x(i)<180
sum1=sum1+x(i);
count1=count1+1;
else
sum2=sum2+x(i);
count2=count2+1;
end
end


if (count1>0)
k1=sum1/count1;
end


if (count2>0)
k2=sum2/count2;
end


if count1>0 && count2>0
if(k2-k1 >= 180)
retVal = ((sum1+sum2)-count2*360)/len;
else
retVal = (sum1+sum2)/len;
end
elseif count1>0
retVal = k1;
else
retVal = k2;
end

您可以在下面的链接中看到任何编程语言的解决方案和一些说明: Https://rosettacode.org/wiki/averages/mean_angle

例如,C + + 解决方案:

#include<math.h>
#include<stdio.h>


double
meanAngle (double *angles, int size)
{
double y_part = 0, x_part = 0;
int i;


for (i = 0; i < size; i++)
{
x_part += cos (angles[i] * M_PI / 180);
y_part += sin (angles[i] * M_PI / 180);
}


return atan2 (y_part / size, x_part / size) * 180 / M_PI;
}


int
main ()
{
double angleSet1[] = { 350, 10 };
double angleSet2[] = { 90, 180, 270, 360};
double angleSet3[] = { 10, 20, 30};


printf ("\nMean Angle for 1st set : %lf degrees", meanAngle (angleSet1, 2));
printf ("\nMean Angle for 2nd set : %lf degrees", meanAngle (angleSet2, 4));
printf ("\nMean Angle for 3rd set : %lf degrees\n", meanAngle (angleSet3, 3));
return 0;
}

产出:

Mean Angle for 1st set : -0.000000 degrees
Mean Angle for 2nd set : -90.000000 degrees
Mean Angle for 3rd set : 20.000000 degrees

Matlab 解决方案:

function u = mean_angle(phi)
u = angle(mean(exp(i*pi*phi/180)))*180/pi;
end


mean_angle([350, 10])
ans = -2.7452e-14
mean_angle([90, 180, 270, 360])
ans = -90
mean_angle([10, 20, 30])
ans =  20.000

虽然 Starblue 的答案给出了平均单位向量的角度,但是如果你接受在0到2 * pi (或0 ° 到360 °)的范围内可能有不止一个答案,那么就有可能将算术平均值的概念扩展到角度。例如,0 ° 和180 ° 的平均值可能是90 ° 或270 ° 。

算术均值具有与输入值平方和最小的单值性质。两个单位向量之间沿单位圆的距离可以很容易地作为它们的点乘的反余弦计算出来。如果我们选择一个单位向量,通过最小化我们的向量和每个输入单位向量的点积的平方逆余弦之和,那么我们就得到了一个等价的平均值。同样,请记住,在特殊情况下可能有两个或更多的最小值。

这个概念可以扩展到任意数量的维度,因为沿单位球面的距离可以用与沿单位圆的距离完全相同的方法计算——两个单位向量的点积的逆余弦。

对于圆,我们可以用很多方法求解这个平均值,但是我提出了下面的 O (n ^ 2)算法(角度以弧度表示,我避免计算单位向量) :

var bestAverage = -1
double minimumSquareDistance
for each a1 in input
var sumA = 0;
for each a2 in input
var a = (a2 - a1) mod (2*pi) + a1
sumA += a
end for
var averageHere = sumA / input.count
var sumSqDistHere = 0
for each a2 in input
var dist = (a2 - averageHere + pi) mod (2*pi) - pi // keep within range of -pi to pi
sumSqDistHere += dist * dist
end for
if (bestAverage < 0 OR sumSqDistHere < minimumSquareDistance) // for exceptional cases, sumSqDistHere may be equal to minimumSquareDistance at least once. In these cases we will only find one of the averages
minimumSquareDistance = sumSqDistHere
bestAverage = averageHere
end if
end for
return bestAverage

如果所有的角度都在180 ° 之内,那么我们可以使用一个更简单的 O (n) + O (排序)算法(同样使用弧度,避免使用单位向量) :

sort(input)
var largestGapEnd = input[0]
var largestGapSize = (input[0] - input[input.count-1]) mod (2*pi)
for (int i = 1; i < input.count; ++i)
var gapSize = (input[i] - input[i - 1]) mod (2*pi)
if (largestGapEnd < 0 OR gapSize > largestGapSize)
largestGapSize = gapSize
largestGapEnd = input[i]
end if
end for
double sum = 0
for each angle in input
var a2 = (angle - largestGapEnd) mod (2*pi) + largestGapEnd
sum += a2
end for
return sum / input.count

要使用度,只需用180代替圆周率。如果你计划使用更多的维度,那么你很可能不得不使用一个迭代法来求解平均值。

这是一个完全的算术解决方案,使用移动平均值,并注意标准化值。如果所有的角度都在圆的一边(相互之间180 ° 以内) ,那么这个算法速度很快,并且能够给出正确的答案。

它在数学上等价于加入偏移量,将值移入范围(0,180) ,计算平均值,然后减去偏移量。

注释描述了一个特定值在任何给定时间可以取得的范围

// angles have to be in the range [0, 360) and within 180° of each other.
// n >= 1
// returns the circular average of the angles int the range [0, 360).
double meanAngle(double* angles, int n)
{
double average = angles[0];
for (int i = 1; i<n; i++)
{
// average: (0, 360)
double diff = angles[i]-average;
// diff: (-540, 540)


if (diff < -180)
diff += 360;
else if (diff >= 180)
diff -= 360;
// diff: (-180, 180)


average += diff/(i+1);
// average: (-180, 540)


if (average < 0)
average += 360;
else if (average >= 360)
average -= 360;
// average: (0, 360)
}
return average;
}

基于 Alnitak 的回答,我编写了一个 Java 方法来计算多个角度的平均值:

如果你的角度是弧度:

public static double averageAngleRadians(double... angles) {
double x = 0;
double y = 0;
for (double a : angles) {
x += Math.cos(a);
y += Math.sin(a);
}


return Math.atan2(y, x);
}

如果你的角度是角度:

public static double averageAngleDegrees(double... angles) {
double x = 0;
double y = 0;
for (double a : angles) {
x += Math.cos(Math.toRadians(a));
y += Math.sin(Math.toRadians(a));
}


return Math.toDegrees(Math.atan2(y, x));
}

问题非常简单。 1. 确保所有角度都在 -180到180度之间。 把所有非负角加起来,取它们的平均值,然后计算有多少个 加上所有的负角度,取其平均值并计算出有多少个。 3. 取正平均数减去负平均数的差 如果差值大于180,则将差值改为360减去差值。否则就改变差异的标志。请注意,差异总是非负的。 平均角等于总平均加差乘以“权重”,负数除以负数和正数之和

在巨蟒中,角度在[-180,180)之间

def add_angles(a, b):
return (a + b + 180) % 360 - 180


def average_angles(a, b):
return add_angles(a, add_angles(-a, b)/2)

详情:

对于 两个角度的平均值,有两个平均值相距180 ° ,但我们可能需要更接近的平均值。

在视觉上,蓝色(B)和绿色()的平均值产生青绿色的点:

Original

角度’环绕’(例如355 + 10 = 5) ,但是标准的算术会忽略这个分支点。 然而,如果角度 B是相对于分支点,那么(B + )/2给出最接近的平均值: 青绿点。

对于任意两个角度,我们可以旋转问题,使其中一个角度相对于分支点,执行标准平均,然后旋转回来。

rotatedreturned

这里是一些 Java 代码,平均角度,我认为它是相当健壮的。

public static double getAverageAngle(List<Double> angles)
{
// r = right (0 to 180 degrees)


// l = left (180 to 360 degrees)


double rTotal = 0;
double lTotal = 0;
double rCtr = 0;
double lCtr = 0;


for (Double angle : angles)
{
double norm = normalize(angle);
if (norm >= 180)
{
lTotal += norm;
lCtr++;
} else
{
rTotal += norm;
rCtr++;
}
}


double rAvg = rTotal / Math.max(rCtr, 1.0);
double lAvg = lTotal / Math.max(lCtr, 1.0);


if (rAvg > lAvg + 180)
{
lAvg += 360;
}
if (lAvg > rAvg + 180)
{
rAvg += 360;
}


double rPortion = rAvg * (rCtr / (rCtr + lCtr));
double lPortion = lAvg * (lCtr / (lCtr + rCtr));
return normalize(rPortion + lPortion);
}


public static double normalize(double angle)
{
double result = angle;
if (angle >= 360)
{
result = angle % 360;
}
if (angle < 0)
{
result = 360 + (angle % 360);
}
return result;
}

我已经迟到很久了,但是我觉得我应该加上我的两分钱,因为我实在找不到任何确切的答案。最后,我实现了 Mitsuta 方法的以下 Java 版本,我希望它能提供一个简单而健壮的解决方案。特别是当标准差提供一个测量色散,并且如果 sd = = 90,表明输入角度导致一个模糊的平均值时。

编辑: 实际上,我意识到我最初的实现可以进一步简化,事实上,考虑到其他答案中的所有对话和三角函数,简化得令人担忧。

/**
* The Mitsuta method
*
* @param angles Angles from 0 - 360
* @return double array containing
* 0 - mean
* 1 - sd: a measure of angular dispersion, in the range [0..360], similar to standard deviation.
* Note if sd == 90 then the mean can also be its inverse, i.e. 360 == 0, 300 == 60.
*/
public static double[] getAngleStatsMitsuta(double... angles) {
double sum = 0;
double sumsq = 0;
for (double angle : angles) {
if (angle >= 180) {
angle -= 360;
}
sum += angle;
sumsq += angle * angle;
}


double mean = sum / angles.length;
return new double[]{mean <= 0 ? 360 + mean: mean, Math.sqrt(sumsq / angles.length - (mean * mean))};
}

对于你们这些 Java 爱好者来说,你们可以用上面的方法得到一行中的平均角度。

Arrays.stream(angles).map(angle -> angle<180 ? angle: (angle-360)).sum() / angles.length;

如果有人正在寻找这方面的 JavaScript解决方案,我已经在 数学图书馆的帮助下将维基百科页面 圆形量的平均值(也在 尼克的回答中提到)中给出的示例翻译成了 JavaScript/NodeJS 代码。

如果你的角度是在 :

const maths = require('mathjs');


getAverageDegrees = (array) => {
let arrayLength = array.length;


let sinTotal = 0;
let cosTotal = 0;


for (let i = 0; i < arrayLength; i++) {
sinTotal += maths.sin(array[i] * (maths.pi / 180));
cosTotal += maths.cos(array[i] * (maths.pi / 180));
}


let averageDirection = maths.atan(sinTotal / cosTotal) * (180 / maths.pi);


if (cosTotal < 0) {
averageDirection += 180;
} else if (sinTotal < 0) {
averageDirection += 360;
}


return averageDirection;
}

为了从一组 指南针方向中找到平均方向,这个解决方案对我来说非常有效。我已经在大范围的方向性数据(0-360度)上测试了它,它看起来非常健壮。

或者,如果你的角度是在 弧度:

const maths = require('mathjs');
getAverageRadians = (array) => {
let arrayLength = array.length;


let sinTotal = 0;
let cosTotal = 0;


for (let i = 0; i < arrayLength; i++) {
sinTotal += maths.sin(array[i]);
cosTotal += maths.cos(array[i]);
}


let averageDirection = maths.atan(sinTotal / cosTotal);


if (cosTotal < 0) {
averageDirection += 180;
} else if (sinTotal < 0) {
averageDirection += 360;
}


return averageDirection;
}

希望这些解决方案对我遇到类似编程挑战的人有所帮助。