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

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

A questão real é mais complicada - o que as estatísticas 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 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 dos ângulos e calcule o ângulo da média.

Esta questão é examinada em detalhes no livro: "Statistics On Spheres", Geoffrey S. Watson, Universidade de Arkansas, Lecture Notes in the Mathematics Sciences, 1983 John Wiley & Sons, Inc., como mencionado em http: //catless.ncl. ac.uk/Risks/7.44.html#subj4 por 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 seus motivos são mais claros e provavelmente programaticamente mais eficientes, e também funcionam bem no caso zero, então parabéns a ele.

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

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ê deseja é realmente 0'.

Eu acho que Starblue gosta de algo. Apenas calcule as coordenadas cartesianas (x, y) para cada ângulo e adicione 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)

Estou ignorando agora que uma direção 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 depois no sentido anti-horário. A matemática deve funcionar da mesma maneira, independentemente.

PARA O CASO ESPECIAL DE DOIS ANGLES:

A resposta ((a + b) mod 360) / 2 está 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: esta equação pega a resposta no sentido horário de a)
• 180, 0 -> 270 (veja acima)
• 180, 1 -> 90,5
• 1, 180 -> 90,5
• 20, 350 -> 5
• 350, 20 -> 5 (todos os exemplos a seguir também revertem corretamente)
• 10, 20 -> 15
• 350, 2 -> 356
• 359, 0 -> 359,5
• 180, 180 -> 180

A ackb está certa de que essas soluções baseadas em vetores não podem ser consideradas médias reais dos ângulos, elas são apenas uma média das contrapartes dos vetores unitários. No entanto, a solução sugerida pelo ackb não parece matematicamente correta.

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

Primeiro, precisamos examinar exatamente em quais casos a diferença entre os ângulos é diferente da diferença entre os números 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 estamos minimizando apenas as mudanças nos pontos em que essas desigualdades mudam de verdadeiro para falso ou vice-versa, podemos separar o intervalo [0,360) completo em um conjunto de segmentos, separados por esses pontos. Então, precisamos encontrar apenas 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.

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

Aqui aplicaremos a idéia de minimizar a diferença ao quadrado para derivar a fórmula da média aritmética comum: sum (a [i]) / n. A curva y = soma ((a [i] -x) ^ 2) pode ser minimizada desta 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 a curvas com nossas diferenças ajustadas:

b = subconjunto de a onde a diferença (angular) correta a [i] -xc = subconjunto de a onde a diferença (angular) correta (a [i] -360) -x cn = tamanho do cd = subconjunto de a onde o diferença (angular) correta (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 trabalha para valores normais, que tem um conjunto ilimitado, portanto o resultado definitivamente ficará 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 deverá estar no limite inferior (porque as curvas quadráticas têm apenas 1 ponto de viragem) e se o mínimo for maior que o limite superior do segmento, o mínimo do segmento estará no limite superior. Depois de termos o mínimo para cada segmento, simplesmente encontramos o que tem o menor valor para o que estamos minimizando (soma ((b [i] -x) ^ 2) + soma (((c [i] -360 ) -b) ^ 2) + soma (((d [i] +360) -c) ^ 2)).

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

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ê substituir a classificação baseada em comparação por uma classificação não baseada em comparação, como a classificação 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<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_;
}

A média aritmética de um conjunto de ângulos pode não concordar com sua ideia intuitiva de qual deve ser a média. Por exemplo, a média aritmética do conjunto {179,179,0,181,181} é 216 (e 144). A resposta em que você pensa imediatamente é provavelmente 180, no entanto, é sabido que a média aritmética é fortemente afetada pelos valores das arestas. Lembre-se também de que ângulos não são vetores, por mais atraente que possa parecer ao lidar com ângulos algumas vezes.

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

Eu também gostaria de enfatizar que, embora essa seja uma média real de ângulos, diferentemente das soluções vetoriais, isso não significa necessariamente que é a solução que você deve usar, a média dos vetores unitários correspondentes pode muito bem 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, e não 180. Para fazer isso, é necessário verificar se a diferença entre os dois ângulos é maior que 180 ou não. Nesse caso, aumente o ângulo menor em 360 antes de usar a fórmula acima.

Não vejo como a segunda alternativa possa ser generalizada para o caso de mais de dois ângulos.

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

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

Duas opções comuns de métricas são as métricas de Frobenius e Riemann. Para a métrica Frobenius, existe uma fórmula direta que corresponde à noção usual de rumo médio nas estatísticas circulares. Veja "Médias e médias no grupo de rotações", Maher Moakher, SIAM Journal on Matrix Analysis and Applications, Volume 24, Edição 1, 2002, para obter detalhes.
http://link.aip.org/link/?SJMAEL/24/1/1

Aqui está uma função 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=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, [email protected]
%     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 usei com um microcontrolador que não possui recursos de ponto flutuante ou trigonometria. Eu ainda precisava "mediar" 10 leituras brutas de rolamentos para suavizar as variações.

1. Verifique se o primeiro rolamento está no 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 na faixa de 0 <= rolamento <360. Caso contrário, faça as leituras como elas vêm.
3. Uma vez realizadas 10 leituras, 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 é o ideal; Isso pode quebrar. Eu me safei nesse caso porque o dispositivo gira apenas muito lentamente. Vou divulgá-lo caso alguém mais se veja trabalhando sob restrições semelhantes.

Em inglês:

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

Em python:

Uma matriz de ângulos NX1 #numpy

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

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

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

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 existem 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 verde-azulado:

Os ângulos 'envolvem' (por exemplo, 355 + 10 = 5), mas a aritmética padrão ignorará esse ponto de ramificação. No entanto, se o ângulo b for oposto ao ponto de ramificação, então ( b + g ) / 2 fornece a média mais próxima: o ponto de cerceta.

Para quaisquer dois ângulos, podemos rotacionar o problema para que um dos ângulos fique oposto ao ponto de ramificação, faça a média padrão e depois gire para trás.

Eu seguiria o caminho do vetor usando números complexos. Meu exemplo está no Python, que possui números complexos internos:

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)

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

Aqui está uma solução completa em 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);
}

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

Com base na resposta de Alnitak , escrevi um método Java para calcular a média de vários ângulos:

Se seus ângulos estão em radianos:

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);
}

Se seus ângulos estão em graus:

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));
}

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

Outra idéia: encontre 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 assumir que todos esses pontos caem na mesma metade do círculo? (Caso contrário, não há uma 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 distribuir uniformemente pontos?)

Com essa suposição, escolhemos um ponto arbitrário nesse semicírculo como a "origem" e medimos o conjunto de ângulos especificado em relação a essa origem (chamamos isso de "ângulo relativo"). Observe que o ângulo relativo tem um valor absoluto estritamente menor que 180 graus. Por fim, calcule a média desses ângulos relativos para obter o ângulo médio desejado (em relação à nossa origem, é claro).

Não existe uma única "resposta certa". Eu recomendo a leitura do livro, KV Mardia e PE Jupp, "Directional Statistics", (Wiley, 1999), para uma análise completa.

(Apenas quero compartilhar meu ponto de vista da teoria da estimativa ou inferência estatística)

O julgamento de Nimble é obter a estimativa MMSE ^ de um conjunto de ângulos, mas é uma das opções encontrar uma direção "média"; também é possível encontrar uma estimativa MMAE ^, ou alguma outra estimativa como a direção "média", e isso depende do seu erro de direção de quantificação métrica; ou, mais geralmente, na teoria das estimativas, a definição da função de custo.

^ MMSE / MMAE corresponde ao erro médio quadrático / absoluto mínimo.

ackb disse "O ângulo médio phi_avg deve ter a propriedade que sum_i | phi_avg-phi_i | ^ 2 se torna mínima ... eles calculam a média de algo, mas não de ângulos"

---- você quantifica erros no sentido quadrático médio e é uma das maneiras mais comuns, no entanto, não é a única. A resposta favorecida pela maioria das pessoas aqui (ou seja, soma dos vetores unitários e obter o ângulo do resultado) é na verdade uma das soluções razoáveis. É (pode ser provado) o estimador de ML que serve como a direção "média" que queremos, se as direções dos vetores forem modeladas como distribuição de von Mises. Essa distribuição não é sofisticada e é apenas uma distribuição periodicamente amostrada de um guassiano 2D. Veja Eqn. (2.179) no livro de Bishop "Reconhecimento de padrões e aprendizado de máquina". 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 "a ackb está certa de que essas soluções baseadas em vetores não podem ser consideradas médias reais de ângulos, elas são apenas uma média das contrapartes em vetores unitários"

----isso não é verdade. O "vetor unitário equivalente" revela as informações 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 de que o comprimento é 1. Você pode definir o vetor "unitário" como sendo o comprimento 2, isso realmente não importa.

Aqui está uma solução completamente aritmética usando médias móveis e tomando o cuidado de normalizar valores. É rápido e fornece respostas corretas se todos os ângulos estiverem em um lado do círculo (a 180 ° um do outro).

É matematicamente equivalente a adicionar o deslocamento que altera os valores para o intervalo (0, 180), calcular a média e subtrair o deslocamento.

Os comentários descrevem qual intervalo um valor específico pode assumir a qualquer momento

// 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;
}

Bem, estou muito atrasado para a festa, mas pensei em adicionar meus 2 centavos, pois não consegui encontrar nenhuma resposta definitiva. No final, implementei a seguinte versão Java do método Mitsuta que, espero, forneça uma solução simples e robusta. Especialmente porque o desvio padrão fornece uma dispersão de medida e, se sd == 90, indica que os ângulos de entrada resultam em uma média ambígua.

Edição: Na verdade, percebi que minha implementação original pode ser ainda mais simplificada, de fato preocupantemente simples, considerando toda a conversa e trigonometria acontecendo nas outras respostas.

/**
 * 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))};
}

... e para todos os seus geeks (Java), você pode usar a abordagem acima para obter o ângulo médio em uma linha.

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

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

Para o caso especial de onde

(soma (x_componente) = 0,0 && soma (y_componente) = 0,0) // por exemplo, 2 ângulos de 10. e 190. graus ea.

use 0,0 graus como a soma

Computacionalmente, você deve testar esse caso, pois atan2 (0, 0) é indefinido e gera 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 deve estar em [-Pi, Pi) (porque pode ser mais curto para o contrário!). Isso é facilmente alcançado normalizando todos os valores de entrada para [0, 2Pi), mantendo um phi_run médio em execução e escolhendo normalizando | phi_i-phi_run | para [-Pi, Pi) (adicionando ou subtraindo 2Pi). A maioria das sugestões acima faz outra coisa que não tem essa propriedade mínima, ou seja, elas têm uma média de alguma coisa , mas não ângulos.

Resolvi o problema com a ajuda da resposta de @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, e não 180. Para fazer isso, é necessário verificar se a diferença entre os dois ângulos é maior que 180 ou não. Nesse caso, aumente 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 menores que isso, aumentam em 360. Em seguida, recalcule a média adicionando todos eles e dividindo-os pelo 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 função 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 

Você pode ver uma solução e uma pequena explicação no link a seguir, para QUALQUER linguagem de programação: https://rosettacode.org/wiki/Averages/Mean_angle

Por exemplo, solução 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;
}

Resultado:

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

Ou solução 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

Enquanto a resposta do starblue fornece o ângulo do vetor unitário médio, é possível estender o conceito da média aritmética para ângulos se você aceitar que pode haver mais de uma resposta no intervalo de 0 a 2 * pi (ou 0 ° a 360 °). Por exemplo, a média de 0 ° e 180 ° pode ser 90 ° ou 270 °.

A média aritmética tem a propriedade de ser o valor único com a soma mínima de distâncias ao quadrado dos valores de entrada. A distância ao longo do círculo unitário entre dois vetores unitários pode ser facilmente calculada como o cosseno inverso do seu produto escalar. Se escolhermos um vetor unitário minimizando a soma do cosseno inverso ao quadrado do produto escalar de nosso vetor e cada vetor unitário de entrada, teremos uma média equivalente. Novamente, lembre-se de que pode haver dois ou mais mínimos em casos excepcionais.

Esse conceito pode ser estendido a qualquer número de dimensões, uma vez que a distância ao longo da esfera unitária pode ser calculada exatamente da mesma maneira que a distância ao longo do círculo unitário - o cosseno inverso do produto escalar de dois vetores unitários.

Para círculos, poderíamos resolver essa média de várias maneiras, mas proponho o seguinte algoritmo O (n ^ 2) (os ângulos estão em radianos e evito calcular os vetores unitários):

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

Se todos os ângulos estiverem a 180 ° um do outro, poderíamos usar um algoritmo O (n) + O (classificação) mais simples (novamente usando radianos e evitando o uso de vetores unitários):

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

Para usar graus, basta substituir pi por 180. Se você planeja usar mais dimensões, provavelmente precisará usar um método iterativo para resolver a média.

O problema é extremamente simples. 1. Verifique se todos os ângulos estão entre -180 e 180 graus. 2. a Adicione todos os ângulos não negativos, calcule a média e COUNT quantos 2. b.Adicione todos os ângulos negativos, calcule a média e COUNT quantos. 3. Considere a diferença de pos_average menos neg_average Se a diferença for maior que 180, altere a diferença para 360 menos diferença. Caso contrário, basta alterar o sinal da diferença. Observe que a diferença é sempre não negativa. O Average_Angle é igual a pos_average mais diferença vezes o "peso", contagem negativa dividida pela soma da contagem negativa e positiva

Aqui está um código java para ângulos médios, acho que é razoavelmente robusto.

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

Eu tenho um método diferente do @Starblue que fornece respostas "corretas" para alguns dos ângulos dados acima. Por exemplo:

• angle_avg ([350,10]) = 0
• angle_avg ([- 90,90,40]) = 13,333
• angle_avg ([350,2]) = 356

Ele usa uma soma sobre as diferenças entre ângulos consecutivos. O código (no 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