Jak obliczyć średnią zbioru danych cyklicznych?

Chcę obliczyć średnią zbioru danych cyklicznych. Na przykład mogę mieć kilka próbek z odczytu kompasu. Problem polega oczywiście na tym, jak radzić sobie z zawijaniem. Ten sam algorytm może być przydatny dla tarczy zegara.

Rzeczywiste pytanie jest bardziej skomplikowane - co oznaczają statystyki na kuli lub w przestrzeni algebraicznej, która „otacza”, np. Grupa addytywna mod n. Odpowiedź może nie być jednoznaczna, np. Średnia 359 stopni i 1 stopień może wynosić 0 stopni lub 180, ale statystycznie 0 wygląda lepiej.

To jest dla mnie prawdziwy problem programistyczny i staram się, aby nie wyglądał jak zwykły problem matematyczny.

Oblicz wektory jednostkowe z kątów i weź kąt ich średniej.

Kwestia ta została szczegółowo zbadana w książce: „Statistics On Spheres”, Geoffrey S. Watson, University of Arkansas Lecture Notes in the Mathematical Sciences, 1983 John Wiley & Sons, Inc., jak wspomniano na stronie http: //catless.ncl. ac.uk/Risks/7.44.html#subj4 autor: Bruce Karsh.

Dobry sposób na oszacowanie średniego kąta A na podstawie zestawu pomiarów kątów a [i] 0 <= i

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

Metoda podana przez starblue jest równoważna obliczeniowo, ale jego powody są jaśniejsze i prawdopodobnie bardziej wydajne programowo, a także działają dobrze w przypadku zerowym, więc chwała mu.

Temat jest teraz bardziej szczegółowo omawiany w Wikipedii i przy innych zastosowaniach, takich jak części ułamkowe.

Widzę problem - na przykład, jeśli masz kąt 45 'i kąt 315', „naturalna” średnia będzie równa 180 ', ale żądana wartość to w rzeczywistości 0'.

Myślę, że Starblue jest na czymś. Po prostu oblicz współrzędne kartezjańskie (x, y) dla każdego kąta i dodaj razem otrzymane wektory. Wymaganym wynikiem powinno być przesunięcie kątowe końcowego wektora.

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

Na razie ignoruję, że kurs kompasu zaczyna się na północy i biegnie zgodnie z ruchem wskazówek zegara, podczas gdy „normalne” współrzędne kartezjańskie zaczynają się od zera wzdłuż osi X, a następnie idą w kierunku przeciwnym do ruchu wskazówek zegara. Niezależnie od tego matematyka powinna działać w ten sam sposób.

DLA SZCZEGÓLNEGO PRZYPADKU DWÓCH KĄTÓW:

Odpowiedź ((a + b) mod 360) / 2 jest ZŁA . Dla kątów 350 i 2 najbliższy punkt to 356, a nie 176.

Wektor jednostkowy i rozwiązania trygonometryczne mogą być zbyt drogie.

To, co otrzymałem z małego majsterkowania, to:

diff = ( ( a - b + 180 + 360 ) mod 360 ) - 180
angle = (360 + b + ( diff / 2 ) ) mod 360

• 0, 180 -> 90 (dwie odpowiedzi: to równanie przyjmuje odpowiedź zgodną z ruchem wskazówek zegara z a)
• 180, 0 -> 270 (patrz wyżej)
• 180, 1 -> 90,5
• 1, 180 -> 90,5
• 20, 350 -> 5
• 350, 20 -> 5 (wszystkie poniższe przykłady również odwracają się prawidłowo)
• 10, 20 -> 15
• 350, 2 -> 356
• 359, 0 -> 359,5
• 180, 180 -> 180

Ackb ma rację, że tych rozwiązań opartych na wektorach nie można uznać za prawdziwe średnie kątów, są one jedynie średnią odpowiedników wektorów jednostkowych. Jednak rozwiązanie sugerowane przez ACKB nie wydaje się matematycznie poprawne.

Poniższe rozwiązanie jest matematycznie wyprowadzone z celu minimalizacji (kąt [i] - avgAngle) ^ 2 (gdzie różnica jest korygowana, jeśli to konieczne), co czyni go prawdziwą średnią arytmetyczną kątów.

Po pierwsze, musimy dokładnie przyjrzeć się, w których przypadkach różnica między kątami jest inna niż różnica między ich odpowiednikami z normalną liczbą. Rozważ kąty x i y, jeśli y> = x - 180 i y <= x + 180, to możemy bezpośrednio użyć różnicy (xy). W przeciwnym razie, jeśli pierwszy warunek nie zostanie spełniony, w obliczeniach musimy użyć (y + 360) zamiast y. Odpowiednio, jeśli drugi warunek nie jest spełniony, musimy użyć (y-360) zamiast y. Ponieważ równanie krzywej minimalizujemy tylko zmiany w punktach, w których te nierówności zmieniają się z prawdziwego na fałszywy lub odwrotnie, możemy rozdzielić pełny [0,360) zakres na zestaw segmentów oddzielonych tymi punktami. Następnie musimy tylko znaleźć minimum każdego z tych segmentów, a następnie minimum minimum każdego segmentu, które jest średnią.

Oto obraz pokazujący, gdzie występują problemy przy obliczaniu różnic kątów. Jeśli x znajduje się w szarym obszarze, wystąpi problem.

Aby zminimalizować zmienną, w zależności od krzywej, możemy wziąć pochodną tego, co chcemy zminimalizować, a następnie znaleźć punkt zwrotny (gdzie pochodna = 0).

Tutaj zastosujemy ideę zminimalizowania kwadratowej różnicy w celu wyprowadzenia wspólnego wzoru na średnią arytmetyczną: suma (a [i]) / n. Krzywą y = sum ((a [i] -x) ^ 2) można zminimalizować w ten sposób:

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

Teraz zastosuj to do krzywych z naszymi dostosowanymi różnicami:

b = podzbiór a gdzie poprawna (kątowa) różnica a [i] -xc = podzbiór a gdzie poprawna (kątowa) różnica (a [i] -360) -x cn = wielkość cd = podzbiór a gdzie poprawna (kątowa) różnica (a [i] +360) -x dn = rozmiar 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

Samo to nie wystarczy, aby uzyskać minimum, podczas gdy działa dla normalnych wartości, które mają nieograniczony zbiór, więc wynik z pewnością będzie mieścił się w zakresie zestawu i dlatego jest ważny. Potrzebujemy minimum mieszczącego się w zakresie (określonym przez segment). Jeśli minimum jest mniejsze niż dolna granica naszego segmentu, to minimum tego segmentu musi znajdować się w dolnej granicy (ponieważ krzywe kwadratowe mają tylko 1 punkt zwrotny), a jeśli minimum jest większe niż górna granica naszego segmentu, to minimum segmentu jest na Górna granica. Gdy mamy minimum dla każdego segmentu, po prostu znajdujemy ten, który ma najniższą wartość tego, co minimalizujemy (suma ((b [i] -x) ^ 2) + suma (((c [i] -360 ) -b) ^ 2) + suma (((d [i] +360) -c) ^ 2)).

Oto obraz krzywej, który pokazuje, jak zmienia się w punktach, w których x = (a [i] +180)% 360. Zbiór danych, o którym mowa, to {65,92,230,320,250}.

Oto implementacja algorytmu w Javie, w tym kilka optymalizacji, jego złożoność wynosi O (nlogn). Można go zredukować do O (n), jeśli zastąpisz sortowanie oparte na porównaniu sortowaniem bez porównania, takim jak sortowanie 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_;
}

Średnia arytmetyczna zbioru kątów może nie zgadzać się z intuicyjnym wyobrażeniem o tym, jaka powinna być średnia. Na przykład średnia arytmetyczna zbioru {179,179,0,181,181} wynosi 216 (i 144). Odpowiedź, o której natychmiast pomyślisz, to prawdopodobnie 180, jednak dobrze wiadomo, że na średnią arytmetyczną duży wpływ mają wartości krawędzi. Należy również pamiętać, że kąty nie są wektorami, tak atrakcyjne, jak może się to czasami wydawać, gdy mamy do czynienia z kątami.

Algorytm ten oczywiście ma również zastosowanie do wszystkich wielkości, które są zgodne z arytmetyką modularną (z minimalną korektą), takich jak pora dnia.

Chciałbym również podkreślić, że chociaż jest to prawdziwa średnia kątów, w przeciwieństwie do rozwiązań wektorowych, to niekoniecznie oznacza, że ​​jest to rozwiązanie, którego powinieneś użyć, średnia odpowiednich wektorów jednostkowych może być wartością, którą faktycznie używasz powinien być używany.

Musisz dokładniej zdefiniować średnią . W konkretnym przypadku dwóch kątów mogę wymyślić dwa różne scenariusze:

1. „Prawdziwa” średnia, tj. (A + b) / 2% 360.
2. Kąt, który wskazuje "między" dwoma pozostałymi, pozostając w tym samym półkolu, np. Dla 355 i 5, to będzie 0, a nie 180. Aby to zrobić, musisz sprawdzić, czy różnica między dwoma kątami jest większa niż 180 albo nie. Jeśli tak, zwiększ mniejszy kąt o 360 przed użyciem powyższego wzoru.

Nie rozumiem jednak, jak można uogólnić drugą alternatywę dla przypadku więcej niż dwóch kątów.

Podobnie jak w przypadku wszystkich średnich, odpowiedź zależy od wyboru miernika. Dla danej miary M, średnia niektórych kątów a_k w [-pi, pi] dla k w [1, N] jest tym kątem a_M, który minimalizuje sumę kwadratów odległości d ^ 2_M (a_M, a_k). W przypadku średniej ważonej po prostu uwzględnia się w sumie wagi w_k (takie, że sum_k w_k = 1). To jest,

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

Dwie powszechne metryki to Frobenius i Riemann. W przypadku miernika Frobeniusa istnieje bezpośredni wzór, który odpowiada zwykłemu pojęciu średniego namiaru w statystykach cyklicznych. Szczegółowe informacje można znaleźć w „Means and Averaging in the Group of Rotations”, Maher Moakher, SIAM Journal on Matrix Analysis and Applications, tom 24, wydanie 1, 2002.
http://link.aip.org/link/?SJMAEL/24/1/1

Oto funkcja dla GNU Octave 3.2.4, która wykonuje obliczenia:

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

Chciałbym udostępnić metodę, której użyłem, z mikrokontrolerem, który nie miał możliwości zmiennoprzecinkowych ani trygonometrycznych. Nadal musiałem „uśrednić” 10 nieprzetworzonych odczytów łożysk, aby wyrównać zmiany.

1. Sprawdź, czy pierwszy namiar ma zakres 270-360 lub 0-90 stopni (dwie ćwiartki północne)
2. Jeśli tak, należy obrócić ten i wszystkie kolejne odczyty o 180 stopni, utrzymując wszystkie wartości w zakresie 0 <= namiar <360. W przeciwnym razie odczyty należy odczytywać na bieżąco.
3. Po wykonaniu 10 odczytów obliczyć średnią liczbową, zakładając, że nie nastąpiło zawijanie
4. Jeśli miał miejsce obrót o 180 stopni, należy obrócić obliczoną średnią o 180 stopni, aby powrócić do „prawdziwego” łożyska.

To nie jest idealne; może się zepsuć. W tym przypadku uszło mi to na sucho, ponieważ urządzenie obraca się bardzo wolno. Opublikuję to na wypadek, gdyby ktoś inny pracował z podobnymi ograniczeniami.

Po angielsku:

1. Utwórz drugi zestaw danych ze wszystkimi kątami przesuniętymi o 180.
2. Weź wariancję obu zestawów danych.
3. Weź średnią ze zbioru danych z najmniejszą wariancją.
4. Jeśli ta średnia pochodzi z przesuniętego zestawu, przesuń odpowiedź ponownie o 180.

W Pythonie:

Tablica kątów #numpy NX1

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

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

Oto pełne rozwiązanie: (dane wejściowe to tablica namiarów w stopniach (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;
}

W Pythonie, z kątami pomiędzy [-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)

Detale:

Dla średniej z dwóch kątów istnieją dwie średnie oddalone od siebie o 180 °, ale możemy chcieć bliższej średniej.

Wizualnie, średnia z niebieskiego ( b ) i zielonego ( a ) daje turkusowy punkt:

Kąty „zawijają się” (np. 355 + 10 = 5), ale standardowa arytmetyka ignoruje ten punkt rozgałęzienia. Jeśli jednak kąt b jest przeciwny do punktu rozgałęzienia, to ( b + g ) / 2 daje najbliższą średnią: punkt w kolorze turkusowym.

Dla dowolnych dwóch kątów możemy obrócić zadanie tak, aby jeden z kątów był przeciwny do punktu rozgałęzienia, wykonać standardowe uśrednienie, a następnie obrócić z powrotem.

Poszedłbym drogą wektorową, używając liczb zespolonych. Mój przykład jest w Pythonie, który ma wbudowane liczby zespolone:

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)

Zwróć uwagę, że Python nie musi tworzyć tymczasowej nowej listy wektorów, wszystkie powyższe można wykonać w jednym kroku; Po prostu wybrałem ten sposób, aby przybliżyć pseudokod mający zastosowanie również w innych językach.

Oto kompletne rozwiązanie 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);
}

Przyjmuje kąty w postaci wektorów podwójnych i zwraca średnią po prostu jako podwójną. Kąty muszą być podane w stopniach i oczywiście średnia jest również w stopniach.

Opierając się na odpowiedzi Alnitaka , napisałem metodę Java do obliczania średniej z wielu kątów:

Jeśli twoje kąty są w radianach:

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

Jeśli twoje kąty są w stopniach:

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

Oto pomysł: utwórz średnią iteracyjnie, zawsze obliczając średnią z kątów, które są najbliżej siebie, zachowując wagę.

Inny pomysł: znajdź największą lukę między podanymi kątami. Znajdź punkt, który dzieli ją na pół, a następnie wybierz przeciwległy punkt na okręgu jako zero odniesienia, na podstawie którego obliczysz średnią.

Przedstawmy te kąty punktami na obwodzie koła.

Czy możemy założyć, że wszystkie te punkty leżą na tej samej połowie koła? (W przeciwnym razie nie ma oczywistego sposobu zdefiniowania „średniego kąta”. Pomyśl o dwóch punktach na średnicy, np. 0 stopni i 180 stopni - czy średnia wynosi 90 stopni czy 270 stopni? Co się stanie, gdy mamy 3 lub więcej równomiernie rozłożyć punkty?)

Przy takim założeniu jako „początek” wybieramy dowolny punkt na tym półkolu i mierzymy podany zbiór kątów względem tego początku (nazywamy to „kątem względnym”). Zauważ, że względny kąt ma wartość bezwzględną dokładnie mniejszą niż 180 stopni. Na koniec, weź średnią tych kątów względnych, aby uzyskać pożądany średni kąt (w stosunku do naszego pochodzenia oczywiście).

Nie ma jednej „właściwej odpowiedzi”. Zalecam przeczytanie książki KV Mardia i PE Jupp, „Directional Statistics”, (Wiley, 1999), w celu dokładnej analizy.

(Chcę tylko podzielić się moim punktem widzenia z teorii szacowania lub wnioskowania statystycznego)

Próba Nimble'a polega na uzyskaniu oszacowania MMSE ^ zestawu kątów, ale jest to jeden z wyborów, aby znaleźć „uśredniony” kierunek; można również znaleźć oszacowanie MMAE ^ lub inne oszacowanie jako kierunek „uśredniony”, i zależy to od twojego kwantyfikującego błędu kierunku; lub bardziej ogólnie w teorii estymacji, definicja funkcji kosztu.

^ MMSE / MMAE odpowiada minimalnemu średniemu kwadratowi / błędowi bezwzględnemu.

ACKB powiedział: "Średni kąt phi_avg powinien mieć taką właściwość, że sum_i | phi_avg-phi_i | ^ 2 staje się minimalny ... coś uśrednia, ale nie kąty"

---- kwantyfikujesz błędy w sensie średniokwadratowym i jest to jeden z najczęściej stosowanych sposobów, ale nie jedyny. Odpowiedź preferowana przez większość ludzi tutaj (tj. Suma wektorów jednostkowych i uzyskanie kąta wyniku) jest w rzeczywistości jednym z rozsądnych rozwiązań. Jest to (można udowodnić) estymator ML, który służy jako "uśredniony" kierunek, jaki chcemy, jeśli kierunki wektorów są modelowane jako rozkład von Misesa. Ta dystrybucja nie jest wyszukana i jest po prostu okresowo próbkowaną dystrybucją z Guassian 2D. Zobacz Eqn. (2.179) w książce Bishopa „Rozpoznawanie wzorców i uczenie maszynowe”. Ponownie, w żadnym wypadku nie jest to jedyny najlepszy kierunek do reprezentowania "średniego" kierunku, jednak jest całkiem rozsądny, który ma zarówno dobre uzasadnienie teoretyczne, jak i prostą implementację.

Zwinny powiedział: „potwierdzenie ma rację, że tych rozwiązań opartych na wektorach nie można uznać za prawdziwe średnie kątów, są one jedynie średnią odpowiedników wektorów jednostkowych”

----to nie jest prawda. „Odpowiedniki wektora jednostkowego” ujawniają informacje o kierunku wektora. Kąt jest wielkością bez uwzględnienia długości wektora, a wektor jednostkowy to coś z dodatkową informacją, że długość wynosi 1. Możesz zdefiniować wektor „jednostkowy” tak, aby miał długość 2, to nie ma znaczenia.

Oto całkowicie arytmetyczne rozwiązanie wykorzystujące średnie kroczące i zwracające uwagę na normalizację wartości. Jest szybki i dostarcza poprawnych odpowiedzi, jeśli wszystkie kąty są po jednej stronie koła (w granicach 180 ° od siebie).

Jest to równoważne matematycznie z dodaniem przesunięcia, które przesuwa wartości do zakresu (0, 180), obliczeniem średniej, a następnie odjęciem przesunięcia.

Komentarze opisują, jaki zakres może przyjąć określona wartość w danym momencie

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

Cóż, jestem ogromnie spóźniony na imprezę, ale pomyślałem, że dodam moje 2 centy, ponieważ nie mogłem znaleźć żadnej ostatecznej odpowiedzi. W końcu zaimplementowałem następującą wersję metody Mitsuta w Javie, która, mam nadzieję, zapewnia proste i solidne rozwiązanie. Zwłaszcza, że ​​odchylenie standardowe zapewnia zarówno dyspersję miary, jak i, jeśli sd == 90, wskazuje, że kąty wejściowe dają niejednoznaczną średnią.

EDYCJA: Właściwie zdałem sobie sprawę, że moja oryginalna implementacja może być jeszcze bardziej uproszczona, w rzeczywistości niepokojąco prosta, biorąc pod uwagę całą rozmowę i trygonometrię, które mają miejsce w innych odpowiedziach.

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

... a dla wszystkich maniaków (Java), możesz użyć powyższego podejścia, aby uzyskać średni kąt w jednej linii.

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

Alnitak ma właściwe rozwiązanie. Rozwiązanie Nicka Fortescue jest funkcjonalnie takie samo.

W szczególnym przypadku gdzie

(sum (x_component) = 0,0 && sum (y_component) = 0,0) // np. 2 kąty po 10 i 190 stopni.

użyj sumy 0,0 stopni

Obliczeniowo musisz przetestować ten przypadek, ponieważ atan2 (0., 0) jest niezdefiniowane i spowoduje błąd.

Średni kąt phi_avg powinien mieć taką właściwość, że sum_i | phi_avg-phi_i | ^ 2 staje się minimalny, gdzie różnica musi być w [-Pi, Pi) (ponieważ może być krótsza, aby przejść na drugą stronę!). Można to łatwo osiągnąć poprzez normalizację wszystkich wartości wejściowych do [0, 2Pi), utrzymanie bieżącej średniej phi_run i wybranie normalizacji | phi_i-phi_run | do [-Pi, Pi) (dodając lub odejmując 2Pi). Większość powyższych sugestii robi coś innego, co nie ma tej minimalnej właściwości, tj. Uśrednia coś , ale nie kąty.

Rozwiązałem problem z pomocą odpowiedzi od @David_Hanak. Jak stwierdza:

Kąt, który wskazuje "między" dwoma pozostałymi, pozostając w tym samym półkolu, np. Dla 355 i 5, to będzie 0, a nie 180. Aby to zrobić, musisz sprawdzić, czy różnica między dwoma kątami jest większa niż 180 albo nie. Jeśli tak, zwiększ mniejszy kąt o 360 przed użyciem powyższego wzoru.

Więc obliczyłem średnią wszystkich kątów. A potem wszystkie kąty, które są mniejsze, zwiększ je o 360. Następnie oblicz ponownie średnią, dodając je wszystkie i dzieląc przez ich długość.

        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;

Działa świetnie.

Funkcja Pythona:

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

Możesz użyć tej funkcji w Matlabie:

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 

Możesz zobaczyć rozwiązanie i małe wyjaśnienie w poniższym linku, dla DOWOLNEGO języka programowania: https://rosettacode.org/wiki/Aagues/Mean_angle

Na przykład rozwiązanie 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;
}

Wynik:

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

Lub rozwiązanie 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

Chociaż odpowiedź Starblue'a podaje kąt średniego wektora jednostkowego, możliwe jest rozszerzenie pojęcia średniej arytmetycznej na kąty, jeśli przyjmiesz, że może istnieć więcej niż jedna odpowiedź w zakresie od 0 do 2 * pi (lub od 0 ° do 360 °). Na przykład, średnia 0 ° i 180 ° może wynosić 90 ° lub 270 °.

Średnia arytmetyczna ma tę właściwość, że jest pojedynczą wartością z minimalną sumą kwadratów odległości do wartości wejściowych. Odległość wzdłuż okręgu jednostkowego między dwoma wektorami jednostkowymi można łatwo obliczyć jako odwrotność cosinusa ich iloczynu skalarnego. Jeśli wybierzemy wektor jednostkowy, minimalizując sumę kwadratu odwrotnego cosinusa iloczynu skalarnego naszego wektora i każdego wejściowego wektora jednostkowego, otrzymamy równoważną średnią. Pamiętaj, że w wyjątkowych przypadkach mogą obowiązywać dwa lub więcej minimum.

Pojęcie to można rozszerzyć na dowolną liczbę wymiarów, ponieważ odległość wzdłuż sfery jednostkowej można obliczyć dokładnie w taki sam sposób, jak odległość wzdłuż okręgu jednostkowego - odwrotność cosinusa iloczynu skalarnego dwóch wektorów jednostkowych.

Dla okręgów możemy rozwiązać tę średnią na wiele sposobów, ale proponuję następujący algorytm O (n ^ 2) (kąty są w radianach, a ja unikam obliczania wektorów jednostkowych):

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

Jeśli wszystkie kąty znajdują się w granicach 180 ° od siebie, możemy użyć prostszego algorytmu O (n) + O (sort) (ponownie używając radianów i unikając wektorów jednostkowych):

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

Aby użyć stopni, po prostu zamień pi na 180. Jeśli planujesz użyć większej liczby wymiarów, najprawdopodobniej będziesz musiał użyć metody iteracyjnej w celu obliczenia średniej.

Problem jest niezwykle prosty. 1.Upewnij się, że wszystkie kąty mieszczą się w przedziale od -180 do 180 stopni. 2. a Dodaj wszystkie kąty nieujemne, weź ich średnią i LICZ ile 2. b. Dodaj wszystkie kąty ujemne, weź ich średnią i LICZ ile. 3. Weź różnicę pos_average minus neg_average Jeśli różnica jest większa niż 180, zmień różnicę na 360 minus różnica. W przeciwnym razie po prostu zmień znak różnicy. Zauważ, że różnica jest zawsze nieujemna. Average_Angle równa się pos_average plus różnica razy „waga”, ujemna liczba podzielona przez sumę ujemnej i dodatniej liczby

Oto trochę kodu Java do średnich kątów, myślę, że jest dość solidny.

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

Mam inną metodę niż @Starblue, która daje „poprawne” odpowiedzi na niektóre z powyższych kątów. Na przykład:

• kąt_śr ([350,10]) = 0
• kąt_śr ([- 90,90,40]) = 13,333
• kąt_śr ([350,2]) = 356

Wykorzystuje sumę różnic między kolejnymi kątami. Kod (w Matlabie):

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