Como você calcula a média de um dataset circulares?

Eu quero calcular a média de um dataset circulares. Por exemplo, posso ter várias amostras da leitura de uma bússola. O problema, claro, é como lidar com o embrulho. O mesmo algoritmo pode ser útil para um mostrador de relógio.

A questão atual é mais complicada – o que as statistics significam em uma esfera ou em um espaço algébrico que “envolve”, por exemplo, o grupo aditivo mod n. A resposta pode não ser única, por exemplo, a média de 359 graus e 1 grau pode ser 0 graus ou 180, mas estatisticamente 0 parece melhor.

Este é um problema de programação real para mim e estou tentando fazer com que não pareça apenas um problema de matemática.

Calcule vetores unitários a partir dos ângulos e tome o ângulo da sua média.

Esta questão é examinada em detalhes no livro: “Estatísticas sobre Esferas”, Geoffrey S. Watson, Universidade de Arkansas, Notas de Aula em Ciências Matemáticas, 1983, John Wiley & Sons, Inc. como mencionado em http: // catless.ncl. ac.uk/Risks/7.44.html#subj4 de Bruce Karsh.

Uma boa maneira de estimar um ângulo médio, A, a partir de um conjunto de medidas de ângulo a [i] 0 < = i

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

O método dado por starblue é computacionalmente equivalente, mas suas razões são mais claras e provavelmente programaticamente mais eficientes, e também funcionam bem no caso zero, então muitos elogios para ele.

O assunto agora é explorado com mais detalhes na Wikipédia e com outros usos, como partes fracionárias.

Eu vejo o problema – por exemplo, se você tem um ângulo de 45 ‘e um ângulo de 315’, a média “natural” seria 180 ‘, mas o valor que você quer é na verdade 0’.

Eu acho que o Starblue está em algo. Basta calcular as coordenadas cartesianas (x, y) para cada ângulo e adicionar os vetores resultantes juntos. O deslocamento angular do vetor final deve ser o resultado desejado.

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

Eu estou ignorando por enquanto que um rumo da bússola começa no norte, e vai no sentido horário, enquanto as coordenadas cartesianas “normais” começam com zero ao longo do eixo X, e então seguem no sentido anti-horário. A matemática deve funcionar da mesma maneira, independentemente disso.

PARA O CASO ESPECIAL DE DOIS ÂNGULOS:

A resposta ((a + b) mod 360) / 2 é ERRADA . Para os ângulos 350 e 2, o ponto mais próximo é 356, não 176.

O vetor unitário e as soluções trigonométricas podem ser muito caros.

O que eu tenho de um pouco de mexer é:

 diff = ( ( a - b + 180 + 360 ) mod 360 ) - 180 angle = (360 + b + ( diff / 2 ) ) mod 360 
  • 0, 180 -> 90 (duas respostas para isso: essa equação pega a resposta no sentido horário de a)
  • 180, 0 -> 270 (ver acima)
  • 180, 1 -> 90,5
  • 1, 180 -> 90,5
  • 20, 350 -> 5
  • 350, 20 -> 5 (todos os exemplos seguintes também são invertidos corretamente)
  • 10, 20 -> 15
  • 350, 2 -> 356
  • 359, 0 -> 359,5
  • 180, 180 -> 180

ackb está certo de que essas soluções baseadas em vetores não podem ser consideradas verdadeiras médias de ângulos, elas são apenas uma média das contrapartes vetoriais unitárias. No entanto, a solução sugerida pelo ackb não parece ter som matematicamente.

A seguir, uma solução que é derivada matematicamente da meta de minimizar (ângulo [i] – avgAngle) ^ 2 (onde a diferença é corrigida, se necessário), o que a torna uma verdadeira média aritmética dos ângulos.

Primeiro, precisamos ver exatamente quais casos a diferença entre os ângulos é diferente da diferença entre os seus correspondentes numéricos normais. Considere os ângulos x e y, se y> = x – 180 e y < = x + 180, então podemos usar a diferença (xy) diretamente. Caso contrário, se a primeira condição não for atendida, devemos usar (y + 360) no cálculo, em vez de y. Correspondente, se a segunda condição não for atendida, devemos usar (y-360) em vez de y. Como a equação da curva minimiza apenas as mudanças nos pontos em que essas desigualdades mudam de verdade para falso ou vice-versa, podemos separar o intervalo total [0,360] em um conjunto de segmentos, separados por esses pontos. Então, só precisamos encontrar o mínimo de cada um desses segmentos e, em seguida, o mínimo do mínimo de cada segmento, que é a média.

Aqui está uma imagem demonstrando onde os problemas ocorrem no cálculo das diferenças de ângulo. Se x estiver na área cinza, haverá um problema.

Comparações de ângulo

Para minimizar uma variável, dependendo da curva, podemos pegar a derivada do que queremos minimizar e então encontramos o ponto de virada (que é onde a derivada = 0).

Aqui, aplicaremos a idéia de minimizar a diferença quadrática para derivar a fórmula média aritmética comum: sum (a [i]) / n. A curva y = sum ((a [i] -x) ^ 2) pode ser minimizada da seguinte maneira:

 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 

Agora aplicando-o às curvas com nossas diferenças ajustadas:

b = subconjunto de a onde a diferença correta (angular) a [i] -xc = subconjunto de a onde a diferença correta (angular) (a [i] -360) -x cn = tamanho de cd = subconjunto de a onde o diferença correta (angular) (a [i] +360) -x dn = tamanho de 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 

Isso por si só não é suficiente para obter o mínimo, enquanto ele trabalha para valores normais, que tem um conjunto ilimitado, então o resultado ficará definitivamente dentro do intervalo do conjunto e é, portanto, válido. Precisamos do mínimo dentro de um intervalo (definido pelo segmento). Se o mínimo for menor que o limite inferior do nosso segmento, o mínimo desse segmento deve estar no limite inferior (porque as curvas quadráticas têm apenas 1 ponto de virada) e se o mínimo for maior que o limite superior do nosso segmento, o mínimo do segmento é limite superior. Depois que tivermos o mínimo para cada segmento, simplesmente encontraremos aquele que tem o menor valor para o que estamos minimizando (sum ((b [i] -x) ^ 2) + sum (((c [i] -360 ) -b) ^ 2) + sum (((d [i] +360) -c) ^ 2)).

Aqui está uma imagem para a curva, que mostra como ela muda nos pontos em que x = (a [i] +180)% 360. O dataset em questão é {65,92,230,320,250}.

Curva

Aqui está uma implementação do algoritmo em Java, incluindo algumas otimizações, sua complexidade é O (nlogn). Ele pode ser reduzido para O (n) se você replace a sorting baseada em comparação por uma sorting não baseada em comparação, como sorting radix.

 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 lowerAngles_ = new LinkedList(); List upperAngles_ = new LinkedList(); 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 lowerAnglesIter = lowerAngles_.iterator(); for(int i = 0; i < lowerAngles_.size(); i++) lowerAngles[i] = lowerAnglesIter.next(); upperAngles = new double[upperAngles_.size()]; Iterator upperAnglesIter = upperAngles_.iterator(); for(int i = 0; i < upperAngles_.size(); i++) upperAngles[i] = upperAnglesIter.next(); } List averageAngles = new LinkedList(); 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 averageAnglesIter = averageAngles.iterator(); for(int i = 0; i < averageAngles_.length; i++) averageAngles_[i] = averageAnglesIter.next(); return averageAngles_; } 

A média aritmética de um conjunto de ângulos pode não concordar com a sua ideia intuitiva do que deve ser a média. Por exemplo, a média aritmética do conjunto {179,179,0,181,181} é 216 (e 144). A resposta que você pensa imediatamente é provavelmente de 180, no entanto, é bem sabido que a média aritmética é fortemente afetada pelos valores de borda. Você também deve lembrar que os ângulos não são vetores, por mais atraentes que possam parecer quando lidamos com ângulos às vezes.

Esse algoritmo também se aplica a todas as grandezas que obedecem à aritmética modular (com ajuste mínimo), como a hora do dia.

Eu também gostaria de enfatizar que mesmo que esta seja uma verdadeira média de ângulos, ao contrário das soluções vetoriais, isso não significa necessariamente que é a solução que você deveria estar usando, a média dos vetores unitários correspondentes pode ser o valor que você realmente deveria estar usando.

Você precisa definir a média com mais precisão. Para o caso específico de dois ângulos, posso pensar em dois cenários diferentes:

  1. A média “verdadeira”, ou seja, (a + b) / 2% 360.
  2. O ângulo que aponta “entre” os outros dois enquanto permanece no mesmo semicírculo, por exemplo, para 355 e 5, seria 0, não 180. Para fazer isso, você precisa verificar se a diferença entre os dois ângulos é maior que 180. ou não. Em caso afirmativo, incremente o ângulo menor em 360 antes de usar a fórmula acima.

Eu não vejo como a segunda alternativa pode ser generalizada para o caso de mais de dois ângulos, no entanto.

Como todas as médias, a resposta depende da escolha da métrica. Para uma dada métrica M, a média de alguns ângulos a_k em [-pi, pi] para k em [1, N] é aquele ângulo a_M que minimiza a sum das distâncias quadradas d ^ 2_M (a_M, a_k). Para uma média ponderada, basta include na sum os pesos w_k (tal que sum_k w_k = 1). Isso é,

a_M = arg min_x sum_k w_k d ^ 2_M (x, a_k)

Duas escolhas comuns de métrica são as métricas Frobenius e Riemann. Para a métrica de Frobenius, existe uma fórmula direta que corresponde à noção usual de porte médio em estatística circular. Ver “Meios e Médias no Grupo de Rotações”, Maher Moakher, Jornal SIAM de Análise e Aplicações de Matrizes, Volume 24, Edição 1, 2002, para detalhes.
http://link.aip.org/link/?SJMAEL/24/1/1

Aqui está uma function para o GNU Octave 3.2.4 que faz o cálculo:

 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/'. 

Eu gostaria de compartilhar um método que eu usei com um microcontrolador que não tinha resources de ponto flutuante ou trigonometria. Eu ainda precisava “mediar” 10 leituras de rolamentos para suavizar as variações.

  1. Verifique se o primeiro rolamento é o intervalo 270-360 ou 0-90 graus (dois quadrantes do norte)
  2. Se estiver, gire esta e todas as leituras subseqüentes em 180 graus, mantendo todos os valores no intervalo 0 < = rolamento <360. Caso contrário, faça as leituras como elas são fornecidas.
  3. Uma vez que 10 leituras foram feitas calcule a média numérica assumindo que não houve
  4. Se a rotação de 180 graus estiver em vigor, gire a média calculada em 180 graus para voltar a um rolamento “verdadeiro”.

Não é ideal; pode quebrar. Eu escapei com isso neste caso porque o dispositivo gira apenas muito devagar. Vou colocá-lo lá para o caso de alguém se encontrar trabalhando sob restrições semelhantes.

Aqui está a solução completa: (a input é uma matriz de rolamento em graus (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; } 

Eu usaria o vetor usando números complexos. Meu exemplo é em Python, que possui números complexos incorporados:

 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) 

Note que o Python não precisa construir uma nova lista temporária de vetores, tudo isso pode ser feito em uma única etapa; Eu escolhi este caminho para aproximar o pseudo-código aplicável a outros idiomas também.

Aqui está uma solução completa em C ++:

 #include  #include  double dAngleAvg(const vector& 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); } 

Ele pega os ângulos na forma de um vetor de duplos e retorna a média simplesmente como um duplo. Os ângulos devem estar em graus e, é claro, a média também é em graus.

Em python, com ângulos entre [-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) 

Detalhes:

Para a média de dois ângulos, há duas médias separadas por 180 °, mas podemos querer a média mais próxima.

Visualmente, a média do azul ( b ) e verde ( a ) produz o ponto da cerceta:

Original

Angles ‘wrap around’ (por exemplo, 355 + 10 = 5), mas a aritmética padrão irá ignorar este ponto de ramificação. No entanto, se o ângulo b é oposto ao ponto de ramificação, então ( b + g ) / 2 dá a média mais próxima: o ponto de cerceta.

Para quaisquer dois ângulos, podemos girar o problema de modo que um dos ângulos fique em frente ao ponto de ramificação, faça a média padrão e depois gire de volta.

rodado devolvida

Aqui está uma ideia: construir a média iterativamente calculando sempre a média dos ângulos mais próximos, mantendo um peso.

Outra ideia: encontrar o maior espaço entre os ângulos dados. Encontre o ponto que o divide e, em seguida, escolha o ponto oposto no círculo como o zero de referência para calcular a média.

Vamos representar esses ângulos com pontos na circunferência do círculo.

Podemos presumir que todos esses pontos caem na mesma metade do círculo? (Caso contrário, não há nenhuma maneira óbvia de definir o “ângulo médio”. Pense em dois pontos no diâmetro, por exemplo, 0 graus e 180 graus — é a média de 90 graus ou 270 graus? O que acontece quando temos 3 ou mais espalhar uniformemente pontos?)

Com essa suposição, escolhemos um ponto arbitrário nesse semicírculo como a “origem” e medimos o conjunto de ângulos dado em relação a essa origem (chame isso de “ângulo relativo”). Note que o ângulo relativo tem um valor absoluto estritamente menor que 180 graus. Finalmente, tome a média desses ângulos relativos para obter o ângulo médio desejado (em relação à nossa origem, é claro).

Não há uma única “resposta certa”. Eu recomendo ler o livro, KV Mardia e PE Jupp, “Directional Statistics”, (Wiley, 1999), para uma análise completa.

Alnitak tem a solução certa. A solução de Nick Fortescue é funcionalmente a mesma.

Para o caso especial de onde

(sum (x_component) = 0.0 && sum (y_component) = 0.0) // por exemplo, 2 ângulos de 10. e 190. graus ea.

use 0,0 graus como a sum

Computationally você tem que testar para este caso desde atan2 (0., 0) é indefinido e irá gerar um erro.

O ângulo médio phi_avg deve ter a propriedade que sum_i | phi_avg-phi_i | ^ 2 se torna mínima, onde a diferença tem que ser em [-Pi, Pi) (porque pode ser mais curto fazer o caminho inverso!). Isso é facilmente alcançado normalizando todos os valores de input para [0, 2Pi), mantendo uma média de execução phi_run e escolhendo normalizar | phi_i-phi_run | para [-Pi, Pi) (adicionando ou subtraindo 2Pi). A maioria das sugestões acima faz outra coisa que não tem aquela propriedade mínima, isto é, eles medem alguma coisa , mas não os ângulos.

Em inglês:

  1. Faça um segundo dataset com todos os ângulos deslocados em 180.
  2. Tome a variação de ambos os conjuntos de dados.
  3. Tome a média do dataset com a menor variação.
  4. Se esta média for do conjunto alterado, mude a resposta novamente por 180.

Em python:

Uma matriz #numpy NX1 de ângulos

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

(Só quero compartilhar meu ponto de vista da Teoria da Estimação ou Inferência Estatística)

O teste de Nimble é obter a estimativa do MMSE de um conjunto de ângulos, mas é uma das opções para encontrar uma direção “média”; pode-se também encontrar uma estimativa do MMAE, ou alguma outra estimativa para ser a direção “média”, e isso depende do seu erro de quantificação métrica; ou mais geralmente na teoria de estimativa, a definição da function de custo.

^ MMSE / MMAE corresponde ao erro médio mínimo ao quadrado / absoluto.

ackb disse: “O phi_avg de ângulo médio deve ter a propriedade que sum_i | phi_avg-phi_i | ^ 2 se torna mínima … eles medem algo, mas não ângulos”

—- você quantifica erros no sentido médio-quadrado e é um dos meios mais comuns, no entanto, não é o único caminho. A resposta favorecida pela maioria das pessoas aqui (isto é, sum dos vetores unitários e obter o ângulo do resultado) é, na verdade, uma das soluções razoáveis. É (pode ser provado) o estimador ML que serve como a direção “média” desejada, se as direções dos vetores forem modeladas como distribuição de von Mises. Esta distribuição não é extravagante, e é apenas uma distribuição amostrada periodicamente de um Guassiano 2D. Veja Eqn. (2.179) no livro de Bishop “Pattern Recognition and Machine Learning”. Novamente, de modo algum é o único melhor para representar a direção “média”, no entanto, é bastante razoável que tenha boa justificativa teórica e implementação simples.

Nimble disse que “ackb está certo de que essas soluções baseadas em vetores não podem ser consideradas verdadeiras médias de ângulos, elas são apenas uma média das contrapartes unitárias de vetores”.

—-isso não é verdade. As “contrapartes vetoriais unitárias” revelam a informação da direção de um vetor. O ângulo é uma quantidade sem considerar o comprimento do vetor, e o vetor unitário é algo com informações adicionais que o comprimento é 1. Você pode definir o seu vetor “unitário” como sendo de comprimento 2, isso realmente não importa.

Eu resolvi o problema com a ajuda da resposta do @David_Hanak. Como ele afirma:

O ângulo que aponta “entre” os outros dois enquanto permanece no mesmo semicírculo, por exemplo, para 355 e 5, seria 0, não 180. Para fazer isso, você precisa verificar se a diferença entre os dois ângulos é maior que 180. ou não. Em caso afirmativo, incremente o ângulo menor em 360 antes de usar a fórmula acima.

Então o que eu fiz foi calcular a média de todos os ângulos. E então todos os ângulos que são menores que isso, aumentam-nos em 360. Em seguida, recalcule a média adicionando-os todos e dividindo-os por seu comprimento.

  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; 

Funciona perfeitamente.

Função 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 

Você pode usar esta function no 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 

You can see a solution and a little explanation in the following link, for ANY programming language: https://rosettacode.org/wiki/Averages/Mean_angle

For instance, C++ solution :

 #include #include 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; } 

Saída:

 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= 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

Here is some java code to average angles, I think it’s reasonably robust.

 public static double getAverageAngle(List 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; } 

I have a different method than @Starblue that gives “correct” answers to some of the angles given above. Por exemplo:

  • 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