我想计算一组循环数据的平均值。例如,我可能有几个指南针读数的样本。问题当然是如何处理环绕。相同的算法可能对表盘有用。
实际问题更复杂 - 统计在球体或“环绕”的代数空间中意味着什么,例如加法组 mod n。答案可能不是唯一的,例如 359 度和 1 度的平均值可能是 0 度或 180 度,但统计上 0 看起来更好。
这对我来说是一个真正的编程问题,我试图让它看起来不仅仅是一个数学问题。
从角度计算单位向量并取其平均值的角度。
这个问题在书中进行了详细研究:“球体统计”,Geoffrey S. Watson,阿肯色大学数学科学讲义,1983 年 John Wiley & Sons, Inc.,如http://catless.ncl 所述。 ac.uk/Risks/7.44.html#subj4,作者 Bruce Karsh。
从一组角度测量值 a[i] 0<=i 估计平均角度 A 的好方法
sum_i_from_1_to_N sin(a[i])
a = arctangent ---------------------------
sum_i_from_1_to_N cos(a[i])
starblue 给出的方法在计算上是等价的,但他的理由更清楚,可能在编程上更有效,并且在零情况下也能很好地工作,所以对他表示敬意。
该主题现在在 Wikipedia 上进行了更详细的探索,并具有其他用途,例如分数部分。
我看到了问题 - 例如,如果您有 45' 角和 315' 角,“自然”平均值将为 180',但您想要的值实际上是 0'。
我认为 Starblue 是在做某事。只需计算每个角度的 (x, y) 笛卡尔坐标,然后将这些结果向量相加。最终向量的角度偏移应该是您需要的结果。
x = y = 0
foreach angle {
x += cos(angle)
y += sin(angle)
}
average_angle = atan2(y, x)
我现在忽略指南针航向从北开始,顺时针方向移动,而“正常”笛卡尔坐标沿 X 轴从零开始,然后逆时针方向移动。无论如何,数学应该以相同的方式进行。
对于两个角度的特殊情况:
答案( (a + b) mod 360 ) / 2是错误的。对于角度 350 和 2,最近点是 356,而不是 176。
单位向量和三角解可能过于昂贵。
我从小修修补补中得到的是:
diff = ( ( a - b + 180 + 360 ) mod 360 ) - 180
angle = (360 + b + ( diff / 2 ) ) mod 360
ackb 是正确的,这些基于向量的解决方案不能被视为角度的真实平均值,它们只是单位向量对应物的平均值。但是,ackb 建议的解决方案在数学上似乎并不合理。
以下是从最小化 (angle[i] - avgAngle)^2 的目标(必要时校正差异)的数学方法得出的解决方案,这使其成为角度的真正算术平均值。
首先,我们需要准确地查看哪些情况下角度之间的差异与其正常数对应物之间的差异不同。考虑角度 x 和 y,如果 y >= x - 180 并且 y <= x + 180,那么我们可以直接使用差 (xy)。否则,如果第一个条件不满足,那么我们必须在计算中使用 (y+360) 而不是 y。相应地,如果不满足第二个条件,那么我们必须使用 (y-360) 代替 y。由于我们最小化的曲线方程仅在这些不等式从真变为假或反之亦然的点处发生变化,因此我们可以将整个 [0,360) 范围分成一组段,由这些点分隔。然后,我们只需要找到这些段的最小值,然后每个段的最小值中的最小值,也就是平均值。
这是一张图片,展示了计算角度差时出现问题的地方。如果 x 位于灰色区域,那么就会有问题。
为了最小化一个变量,根据曲线,我们可以取我们想要最小化的导数,然后我们找到转折点(导数 = 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 的子集,其中正确的(角度)差异 (a[i]-360)-x cn = c 的大小 d = a 的子集,其中正确(角度)差异 (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
仅此一项不足以获得最小值,而它适用于具有无限集合的正常值,因此结果肯定会在集合的范围内,因此是有效的。我们需要一个范围内的最小值(由段定义)。如果最小值小于我们分段的下限,则该分段的最小值必须在下限(因为二次曲线只有 1 个转折点),如果最小值大于我们分段的上限,则该分段的最小值在上限。在我们得到每个段的最小值之后,我们只需找到我们要最小化的值最低的那个 (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}。
这是该算法在 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,但是众所周知,算术平均值受边缘值的影响很大。您还应该记住,角度不是矢量,有时在处理角度时可能看起来很吸引人。
该算法当然也适用于所有服从模运算(调整最少)的量,例如一天中的时间。
我还要强调的是,尽管这是角度的真实平均值,但与矢量解决方案不同,这并不一定意味着它是您应该使用的解决方案,相应单位矢量的平均值很可能是您实际的值应该使用。
您必须更准确地定义平均值。对于两个角度的具体情况,我可以想到两种不同的场景:
不过,我看不出如何将第二种选择推广到两个以上角度的情况。
我想分享一个我在没有浮点或三角函数功能的微控制器上使用的方法。我仍然需要“平均”10 个原始轴承读数以消除变化。
这并不理想;它可以打破。在这种情况下,我侥幸逃脱,因为设备的旋转速度非常慢。我会把它放在那里,以防其他人发现自己在类似的限制下工作。
在 python 中,角度在 [-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 ) 和绿色 ( a )的平均值产生蓝绿色点:
角度“环绕”(例如 355 + 10 = 5),但标准算术将忽略此分支点。但是,如果角度b与分支点相反,则 ( b + g )/2 给出最接近的平均值:蓝绿色点。
对于任意两个角度,我们可以旋转问题,使其中一个角度与分支点相反,执行标准平均,然后旋转回去。
像所有平均值一样,答案取决于指标的选择。对于给定的度量 M,对于 [1,N] 中的 k,[-pi,pi] 中的某些角度 a_k 的平均值是使距离平方和 d^2_M(a_M,a_k) 最小化的角度 a_M。对于加权平均值,只需将权重 w_k 包含在总和中(使得 sum_k w_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/'.
这是完整的解决方案:(输入是以度为单位的方位数组(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;
}
用英语:
在蟒蛇中:
#numpy NX1 角度数组
if np.var(A) < np.var((A-180)%360):
average = np.average(A)
else:
average = (np.average((A-180)%360)+180)%360
我会使用复数采用矢量方式。我的示例是在 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不需要构建一个临时的新向量列表,以上所有操作都可以一步完成;我只是选择这种方式来近似适用于其他语言的伪代码。
这是一个完整的 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);
}
它采用双精度向量形式的角度,并将平均值简单地作为双精度返回。角度必须以度为单位,当然平均值也以度为单位。
基于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));
}
如果有人正在寻找JavaScript解决方案,我已经在mathjs 库的帮助下将维基百科页面Mean of circular quantity(在Nick 的回答中也提到过)中给出的示例翻译成 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;
}
希望这些解决方案对面临与我类似的编程挑战的人有所帮助。
这是一个想法:通过始终计算最接近的角度的平均值来迭代地构建平均值,并保持权重。
另一个想法:找到给定角度之间的最大间隙。找到平分它的点,然后选择圆上的相对点作为参考零来计算平均值。
让我们用圆圆周上的点来表示这些角度。
我们可以假设所有这些点都落在圆的同一半边吗?(否则,没有明显的方法来定义“平均角度”。想想直径上的两个点,例如 0 度和 180 度 --- 是平均 90 度还是 270 度?当我们有 3 个或更多时会发生什么均匀分布点?)
有了这个假设,我们选择该半圆上的任意点作为“原点”,并测量相对于该原点的给定角度集(称为“相对角度”)。请注意,相对角度的绝对值严格小于 180 度。最后,取这些相对角度的平均值以获得所需的平均角度(当然相对于我们的原点)。
没有单一的“正确答案”。我推荐阅读 KV Mardia 和 PE Jupp,“Directional Statistics”一书(Wiley,1999 年),以进行全面分析。
(只是想从估计理论或统计推断中分享我的观点)
Nimble 的试验是获得一组角度的 MMSE^ 估计,但它是找到“平均”方向的一种选择;人们还可以找到一个 MMAE^ 估计,或其他一些估计作为“平均”方向,这取决于您的度量量化方向误差;或更一般地,在估计理论中,成本函数的定义。
^ MMSE/MMAE 对应于最小均方/绝对误差。
ackb 说“平均角度 phi_avg 应该具有 sum_i|phi_avg-phi_i|^2 变得最小的属性......他们平均一些东西,但不是角度”
----您以均方的方式量化误差,这是最常见的方法之一,但不是唯一的方法。这里大多数人喜欢的答案(即单位向量之和并得到结果的角度)实际上是合理的解决方案之一。如果向量的方向被建模为 von Mises 分布,则(可以证明)ML 估计器充当我们想要的“平均”方向。这种分布并不花哨,只是从 2D 高斯分布中定期采样的分布。见方程式。(2.179) 在 Bishop 的《模式识别与机器学习》一书中。同样,它绝不是唯一代表“平均”方向的最佳方法,但是,它是一种非常合理的方法,具有良好的理论依据和简单的实现。
Nimble 说“ackb 是对的,这些基于矢量的解决方案不能被视为角度的真实平均值,它们只是单位矢量对应物的平均值”
- - 这不是真的。“单位向量对应物”揭示了向量的方向信息。角度是一个不考虑向量长度的量,单位向量是带有长度为 1 的附加信息的东西。您可以将“单位”向量定义为长度为 2,这并不重要。
这是一个使用移动平均线并注意标准化值的完全算术解决方案。如果所有角度都在圆的一侧(彼此在 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;
}
好吧,我参加聚会已经很晚了,但我想我会加上我的 2 美分,因为我真的找不到任何明确的答案。最后,我实现了以下 Java 版本的 Mitsuta 方法,我希望它能提供一个简单而健壮的解决方案。特别是因为标准偏差提供了测量分散,并且如果 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;
Alnitak 有正确的解决方案。Nick Fortescue 的解决方案在功能上是相同的。
对于 where 的特殊情况
( sum(x_component) = 0.0 && sum(y_component) = 0.0 ) // 例如 10. 和 190. 度的 2 个角度 ea。
使用 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)。上面的大多数建议都做了一些没有 最小属性的事情,即它们平均一些东西,但不是角度。
我在@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) 算法(再次使用弧度并避免使用单位向量):
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
要使用度数,只需将 pi 替换为 180。如果您打算使用更多维度,那么您很可能必须使用迭代方法来求解平均值。
问题非常简单。1. 确保所有角度都在 -180 到 180 度之间。2. a 将所有非负角相加,取它们的平均值,并计算有多少 2. b. 将所有负角相加,取它们的平均值并计算有多少。3. 取 pos_average 减去 neg_average 的差值 如果差值大于 180,则将差值更改为 360 减去差值。否则,只需更改差异的符号。请注意,差异始终是非负的。Average_Angle 等于 pos_average 加上差异乘以“权重”,负数除以负数和正数之和
这是一些平均角度的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;
}
我有一种与@Starblue 不同的方法,可以对上面给出的某些角度给出“正确”的答案。例如:
它使用对连续角度之间差异的总和。代码(在 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