你如何计算一组循环数据的平均值?
我想计算一组循环数据的平均值。 例如,我可能会从指南针的阅读中得到几个样本。 当然问题是如何处理这个概括。 相同的algorithm可能对钟面有用。
实际的问题更复杂 – 统计数据在一个球体上或者在一个代表性的“代替”的代数空间中意味着什么,比如加法群mod n。 答案可能不是唯一的,例如,359度和1度的平均值可能是0度或180度,但统计上0看起来更好。
这是一个真正的编程问题,我试图让它看起来不仅仅是一个math问题。
从angular度计算单位vector,并取平均angular度。
这个问题在“统计领域”,杰弗里·沃森,阿肯色大学math科学讲义,1983年John Wiley&Sons公司,在http://catless.ncl。 ac.uk/Risks/7.44.html#subj4 by Bruce Karsh。
从一组angular度测量值a [i] 0 <= i估计平均angular度A的一种好方法
sum_i_from_1_to_N sin(a[i]) a = arctangent --------------------------- sum_i_from_1_to_N cos(a[i])
starblue给出的方法在计算上是等价的,但是他的理由更清楚,可能在程序上更有效率,而且在零情况下也能很好地工作,所以对他很有好处。
这个主题现在更详细地在维基百科上进行了探讨,还有其他用途,如小数部分。
我看到了这个问题 – 例如,如果你有45度angular和315度angular,那么“自然”平均值就是180度,但是你想要的值实际上是0度。
我认为Starblue是一些东西。 只需计算每个angular度的(x,y)笛卡尔坐标,并将这些结果vector添加到一起。 最终vector的angular度偏移量应该是您所需的结果。
x = y = 0 foreach angle { x += cos(angle) y += sin(angle) } average_angle = atan2(y, x)
现在我忽略了罗盘航向从北开始顺时针旋转,而“正常”笛卡尔坐标从X轴开始为零,然后逆时针旋转。 math应该以相同的方式工作。
针对两个angular度的特殊情况:
答案((a + b)mod 360)/ 2是错误的 。 对于angular度350和2,最近的点是356,而不是176。
单位vector和触发解决scheme可能太昂贵了。
我从一点修补中得到的是:
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
ackb是正确的,这些基于vector的解决scheme不能被视为真正的angular度平均,他们只是单位vector对应的平均值。 然而,ackbbuild议的解决scheme似乎没有math上的声音。
下面是从最小化(angular度[i] – avgAngle)^ 2(在必要时校正差异)的目标math导出的解决scheme,这使得它成为angular度的真算术平均值。
首先,我们需要看看哪些情况下,angular度之间的差异与正常数字之间的差异是不同的。 考虑angular度x和y,如果y> = x – 180和y <= x + 180,那么我们可以直接使用差异(xy)。 否则,如果不满足第一个条件,则我们必须在计算中使用(y + 360)而不是y。 相应地,如果第二个条件不符合,那么我们必须使用(y-360)而不是y。 由于曲线的方程式中我们只是在这些不等式从真到错或反过来的点上的变化最小化,所以我们可以将完整的[0,360)范围分成由这些点分开的一组段。 然后,我们只需要找出每个分段的最小值,然后每个分段的最小值即平均值。
下面是一张图像,演示计算angular度差异时出现的问题。 如果x位于灰色区域,那么会出现问题。
为了最小化一个variables,根据曲线,我们可以取我们想要最小化的导数,然后find转折点(这是导数= 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 =正确的(angular度)差异的子集a [i] -xc =其中正确(angular度)差异(a [i] -360)的子集-x cn = cd的大小=其中正确(angular度)差异(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
仅仅这一点不足以达到最小值,而它适用于正常值,这是一个无限的集合,所以结果肯定会在集合的范围内,因此是有效的。 我们需要一个范围内的最小值(由段定义)。 如果最小值小于我们线段的下限,那么线段的最小值必须在下限(因为二次曲线只有一个转折点),如果最小值大于我们线段的上限,则线段的最小值在上限。 在我们有每个段的最小值之后,我们只需find最小值的那一个(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中algorithm的实现,包括一些优化,它的复杂度是O(nlogn)。 如果使用非基于比较的sorting(例如基数sorting)replace基于比较的sorting,则它可以简化为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_; }
一组angular度的算术平均值可能不符合你对平均值应该是什么的直觉。 例如,集合{179,179,0,181,181}的算术平均值是216(和144)。 你马上想到的答案可能是180,但是众所周知,算术平均值很大程度上受边缘值的影响。 你还应该记住,angular度不是vector,有时候处理angular度时可能看起来很吸引人。
当然这个algorithm也适用于所有服从模运算(具有最小调整)的量,如一天中的时间。
我还要强调,尽pipe这是一个真正的angular度平均值,与vector解法不同,这不一定意味着它是你应该使用的解决scheme,相应单位vector的平均值可能是你实际上的值应该被使用。
你必须更准确地定义平均值 。 对于两个angular度的具体情况,我可以想到两种不同的情况:
- “真实”的平均值,即(a + b)/ 2%360。
- 在这两个angular度之间的angular度,在两个人之间,而停留在同一个半圆,如355和5,这将是0,而不是180.要做到这一点,你需要检查两个angular度之间的差异是否大于180或不。 如果是这样,在使用上述公式之前,将小angular度增加360°。
不过,我不明白第二种select如何在两个以上的angular度上得到推广。
像所有的平均值一样,答案取决于度量的select。 对于给定的度量M,[1,N]中的[-pi,pi]中的一些angular度a_k的平均值是angular度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,SIAM Journal on Matrix Analysis and Applications,卷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=xz; 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(ab) < pi/2 for some value b % % % na=3, a=sort(mod(randn(1,na)+1,2)-1)*pi; % da=diff([aa(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/'.
我想分享一个我没有使用浮点或三angular函数的微控制器的方法。 我仍然需要“平均”10个原始轴承读数以平滑变化。
- 检查第一个方位是在270-360度还是0-90度(北方两象限)
- 如果是的话,将这个和所有随后的读数旋转180度,保持所有值在0 <=方位<360的范围内。否则,读取数据。
- 一旦读取了10个读数,假定没有环绕,计算数字平均值
- 如果180度旋转已经生效,则将计算的平均值旋转180度,以回到“真实”的方位。
这不是理想的; 它可以打破。 在这种情况下,我摆脱了这种情况,因为设备只能非常缓慢地旋转。 我会把它放在那里,以防其他人发现自己在类似的限制下工作。
我会去使用复数的vector方式。 我的例子是在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 不需要build立一个临时的新的向量列表,以上所有内容都可以一步完成。 我只是select这种方式来近似适用于其他语言的伪代码。
这里是完整的解决scheme:(input是一个轴承的度数(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 ++解决scheme:
#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); }
它以双打vectorforms呈现angular度,并简单地将平均值作为双精度返回。 angular度必须是度数,当然平均也是度数。
这里有一个想法:通过总是计算最接近的angular度的平均值,保持一个重量来迭代地build立平均值。
另一个想法:find给定angular度之间的最大差距。 find平分它的点,然后选取圆上的相反点作为参考零点,以计算平均值。
我们用圆周上的点来表示这些angular度。
我们可以假设所有这些点落在同一半的圈子? (否则,没有明显的方法来定义“平均angular度”,想想直径上的两个点,例如0度和180度—是平均90度还是270度?当我们有3个或更多均匀分布点?)
在这个假设下,我们选取半圆上的一个任意点作为“原点”,并测量相对于这个原点的给定angular度集合(称之为“相对angular度”)。 请注意,相对angular度的绝对值严格小于180度。 最后,取这些相对angular度的平均值,以得到所需的平均angular度(相对于我们当然的起源)。
没有一个“正确答案”。 我build议您阅读KV Mardia和PE Jupp的“方向统计”(Wiley,1999)一书,进行彻底的分析。
Alnitak有正确的解决scheme。 Nick Fortescue的解决scheme在function上是一样的。
对于哪里的特殊情况
(sum(x_component)= 0.0 && sum(y_component)= 0.0)//例如2个angular度为10°和190°ea。
使用0.0度作为总和
计算上你必须testing这种情况,因为atan2(0,,0)是未定义的,将会产生一个错误。
平均angular度phi_avg应该具有sum_i | phi_avg-phi_i | ^ 2变得最小的特性,其中差异必须在[-Pi,Pi)中(因为反过来可能更短)。 通过将所有input值归一化为[0,2Pi),保持移动平均值phi_run并select归一化| phi_i-phi_run | 到[-Pi,Pi)(通过join或减less2Pi)。 上面的大多数build议都做了一些没有这个最小属性的东西 ,也就是说,它们是平均的东西 ,而不是angular度。
用英语:
- 做所有angular度移动180的第二个数据集。
- 取两个数据集的方差。
- 取最小方差的数据集的平均值。
- 如果这个平均值来自移位集合,那么将这个答案再次移动180。
在python中:
#numpy NX1angular度arrays
if np.var(A) < np.var((A-180)%360): average = np.average(A) else: average = (np.average((A-180)%360)+180)%360
(只是想分享我的观点从估计理论或统计推断)
Nimble的试验是获得一组angular度的MMSE ^估计值,但是find一个“平均”方向是一种select。 人们也可以find一个MMAE ^估计值,或者其他一些估计值作为“平均值”的方向,这取决于你的度量方向误差的量度。 或更一般地说在估计理论中,成本函数的定义。
^ MMSE / MMAE对应最小均方差/绝对误差。
ackb说:“平均angular度phi_avg应该具有sum_i | phi_avg-phi_i | ^ 2变得最小的属性…他们平均的东西,而不是angular度”
—-你用均方意义量化误差,这是最常用的方法之一,但不是唯一的方法。 这里大多数人喜欢的答案(即单位向量的总和,并得到结果的angular度)实际上是合理的解决scheme之一。 如果向量的方向被build模为von Mises分布,则(可以certificate)ML估计器作为我们想要的“平均”方向。 这种分布不是花哨的,只是一个二维高斯阶的周期性采样分布。 参见方程 (2.179)在主教的书“模式识别和机器学习”。 同样,它绝不是代表“平均”方向的唯一最好的方法,但是,理论上的合理性和实施的简单性是相当合理的。
Nimble说:“ackb是正确的,这些基于向量的解决scheme不能被视为真正的angular度平均值,它们只是单位向量对应的平均值”
– – 这不是真的。 “单位vector对象”揭示vector方向的信息。 angular度是一个没有考虑vector长度的数量,单位vector是长度为1的附加信息。你可以定义你的“单位”vector长度为2,这并不重要。
我在@David_Hanak的帮助下解决了这个问题。 正如他所说:
在这两个angular度之间的angular度,在两个人之间,而停留在同一个半圆,如355和5,这将是0,而不是180.要做到这一点,你需要检查两个angular度之间的差异是否大于180或不。 如果是这样,在使用上述公式之前,将小angular度增加360°。
所以我做的是计算所有angular度的平均值。 然后所有小于这个的angular度增加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
对于任何编程语言,您可以在以下链接中看到解决scheme和一些解释: https : //rosettacode.org/wiki/Averages/Mean_angle
例如, C ++解决scheme :
#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
Or Matlab solution :
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
While starblue's answer gives the angle of the average unit vector, it is possible to extend the concept of the arithmetic mean to angles if you accept that there may be more than one answer in the range of 0 to 2*pi (or 0° to 360°). For example, the average of 0° and 180° may be either 90° or 270°.
The arithmetic mean has the property of being the single value with the minimum sum of squared distances to the input values. The distance along the unit circle between two unit vectors can be easily calculated as the inverse cosine of their dot product. If we choose a unit vector by minimizing the sum of the squared inverse cosine of the dot product of our vector and each input unit vector then we have an equivalent average. Again, keep in mind that there may be two or more minimums in exceptional cases.
This concept could be extended to any number of dimensions, since the distance along the unit sphere can be calculated in the exact same way as the distance along the unit circle–the inverse cosine of the dot product of two unit vectors.
For circles we could solve for this average in a number of ways, but I propose the following O(n^2) algorithm (angles are in radians, and I avoid calculating the unit vectors):
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
If all the angles are within 180° of each other, then we could use a simpler O(n)+O(sort) algorithm (again using radians and avoiding use of unit vectors):
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
To use degrees, simply replace pi with 180. If you plan to use more dimensions then you will most likely have to use an iterative method to solve for the average.
Here is a completely arithmetic solution using moving averages and taking care to normalize values. It is fast and delivers correct answers if all angles are on one side of the circle (within 180° of each other).
It is mathimatically equivalent to adding the offset which shifts the values into the range (0, 180), calulating the mean and then subtracting the offset.
The comments describe what range a specific value can take on at any given time
// 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; }
Based on Alnitak's answer , I've written a Java method for calculating the average of multiple angles:
If your angles are in radians:
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); }
If your angles are in degrees:
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)); }
The problem is extremely simple. 1. Make sure all angles are between -180 and 180 degrees. 2. a Add all non-negative angles, take their average, and COUNT how many 2. b.Add all negative angles, take their average and COUNT how many. 3. Take the difference of pos_average minus neg_average If difference is greater than 180 then change difference to 360 minus difference. Otherwise just change the sign of difference. Note that difference is always non-negative. The Average_Angle equals the pos_average plus difference times the "weight", negative count divided by the sum of negative and positive count
I have a different method than @Starblue that gives "correct" answers to some of the angles given above. 例如:
- angle_avg([350,10])=0
- angle_avg([-90,90,40])=13.333
- angle_avg([350,2])=356
It uses a sum over the differences between consecutive angles. The code (in 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