我想计算一组循环数据的平均值。例如,我可能有从指南针读数中获取的多个样本。问题在于如何处理循环。同样的算法也可以用于时钟表盘。
实际问题更加复杂 - 统计学在球面或“环绕”的代数空间中的含义是什么,例如模n的加法群。答案可能不唯一,例如359度和1度的平均值可以是0度或180度,但统计上看起来0度更好。
这对我来说是一个真正的编程问题,我正在尝试使它不仅仅像一个数学问题。
我想计算一组循环数据的平均值。例如,我可能有从指南针读数中获取的多个样本。问题在于如何处理循环。同样的算法也可以用于时钟表盘。
实际问题更加复杂 - 统计学在球面或“环绕”的代数空间中的含义是什么,例如模n的加法群。答案可能不唯一,例如359度和1度的平均值可以是0度或180度,但统计上看起来0度更好。
这对我来说是一个真正的编程问题,我正在尝试使它不仅仅像一个数学问题。
从角度计算单位向量,然后取它们的平均角度。
sum_i_from_1_to_N sin(a[i])
a = arctangent ---------------------------
sum_i_from_1_to_N cos(a[i])
starblue提供的方法在计算上与原方法等效,但他的解释更加清晰,可能在编程上更有效,而且在零情况下也能很好地工作,所以要对他表示赞扬。
该主题现在在维基百科上有更详细的探讨,还包括其他用途,如小数部分。
我明白问题所在 - 例如,如果您有一个45度角和一个315度角,"自然"平均值将是180度,但您实际想要的值是0度。
我认为Starblue的想法很有用。只需计算每个角的(x,y)笛卡尔坐标,并将这些结果向量相加即可。最终向量的角偏移应该是您需要的结果。
x = y = 0
foreach angle {
x += cos(angle)
y += sin(angle)
}
average_angle = atan2(y, x)
我现在暂时忽略一个指南针方向从北开始,沿顺时针方向转动,而“普通”的笛卡尔坐标系则从X轴上的零度开始,然后逆时针旋转的事实。无论如何,数学计算应该得出相同的结果。
cos()
、sin()
和atan2()
提供了近似值(虽然很好,但仍然偏离1或2 ulps),因此您平均的次数越多,包含的错误就越多。 - Matthieu对于两个角度的特殊情况:
答案 ((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,则我们可以直接使用差(x-y)。否则,如果不满足第一个条件,则必须在计算中使用(y+360)。相应地,如果第二个条件未满足,则必须使用(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
仅有这些还不足以得到最小值,虽然对于普通值可以使用,但由于它具有无界集合,因此结果肯定在集合的范围内,因此是有效的。我们需要在一个范围内找到最小值(由该段定义)。如果最小值小于我们段的下限,则该段的最小值必须在下限处(因为二次曲线只有一个转折点),如果最小值大于我们段的上限,则该段的最小值在上限处。找到每个段的最小值后,我们只需找到使我们正在最小化的内容的值最低的那个。
下面是一张曲线图,显示了当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,但众所周知,算术平均值受边缘值影响很大。还应记住,角度不是向量,尽管在处理角度时这种理解非常吸引人。在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,在[-pi,pi]中对于k在[1,N]的某些角度a_k,其平均值是那个角度a_M,它最小化平方距离d^2_M(a_M,a_k)的总和。对于加权平均值,只需在总和中包括权重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的文章“Means and Averaging in the Group of Rotations”,SIAM Journal on Matrix Analysis and Applications,Volume 24,Issue 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;
}