Comment calculer la moyenne d'un ensemble de données circulaires ?

Question

Je veux calculer la moyenne d'un ensemble de données circulaires. Par exemple, je peux avoir plusieurs échantillons de la lecture d'une boussole. Le problème, bien sûr, est de savoir comment traiter l'enroulement. Le même algorithme pourrait être utile pour un cadran d'horloge.

La question réelle est plus compliquée - que signifient les statistiques sur une sphère ou dans un espace algébrique qui "s'enroule", par exemple le groupe additif mod n. La réponse peut ne pas être unique, par exemple la moyenne de 359 degrés et 1 degré pourrait être 0 degré ou 180, mais statistiquement 0 semble mieux.

C'est un vrai problème de programmation pour moi et j'essaie de faire en sorte qu'il ne ressemble pas à un simple problème de mathématiques.

Answer 1

Calculez les vecteurs unitaires à partir des angles et prenez l'angle de leur moyenne.

Answer 2

Cette question est examinée en détail dans le livre : "Statistics On Spheres", Geoffrey S. Watson, University of Arkansas Lecture Notes in the Mathematical Sciences, 1983 John Wiley & Sons, Inc. comme mentionné à l'adresse http://catless.ncl.ac.uk/Risks/7.44.html#subj4 par Bruce Karsh.

Une bonne façon d'estimer un angle moyen, A, à partir d'un ensemble de mesures d'angles a[i] 0<=i

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

La méthode donnée par starblue est équivalente du point de vue informatique, mais ses raisons sont plus claires et probablement plus efficaces du point de vue programmatique, et fonctionnent également bien dans le cas zéro, donc bravo à lui.

Le sujet est maintenant exploré plus en détail sur Wikipédia et d'autres usages, comme les parties fractionnées.

Answer 3

Je vois le problème - par exemple, si vous avez un angle de 45' et un angle de 315', la moyenne "naturelle" serait de 180', mais la valeur que vous voulez est en fait 0'.

Je pense que Starblue a trouvé quelque chose. Il suffit de calculer les coordonnées cartésiennes (x, y) pour chaque angle, et d'additionner les vecteurs résultants. Le décalage angulaire du vecteur final devrait être le résultat souhaité.

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

Pour l'instant, je ne tiens pas compte du fait que le cap de la boussole commence au nord et va dans le sens des aiguilles d'une montre, alors que les coordonnées cartésiennes "normales" commencent à zéro sur l'axe X et vont dans le sens inverse des aiguilles d'une montre. Les mathématiques devraient néanmoins donner les mêmes résultats.

Answer 4

POUR LE CAS PARTICULIER DE DEUX ANGLES :

La réponse ( (a + b) mod 360 ) / 2 es MAUVAIS . Pour les angles 350 et 2, le point le plus proche est 356, et non 176.

Les solutions de vecteur unitaire et de trigonométrie peuvent être trop coûteuses.

Ce que j'ai obtenu en bricolant un peu est :

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

• 0, 180 -> 90 (deux réponses pour cela : cette équation prend la réponse de a dans le sens des aiguilles d'une montre)
• 180, 0 -> 270 (voir ci-dessus)
• 180, 1 -> 90.5
• 1, 180 -> 90.5
• 20, 350 -> 5
• 350, 20 -> 5 (tous les exemples suivants s'inversent correctement aussi)
• 10, 20 -> 15
• 350, 2 -> 356
• 359, 0 -> 359.5
• 180, 180 -> 180

Answer 5

Ackb a raison de dire que ces solutions vectorielles ne peuvent pas être considérées comme de véritables moyennes d'angles, elles ne sont qu'une moyenne de leurs contreparties vectorielles unitaires. Cependant, la solution proposée par ackb ne semble pas mathématiquement valable.

La solution suivante est mathématiquement dérivée de l'objectif de minimiser (angle[i] - avgAngle)^2 (où la différence est corrigée si nécessaire), ce qui en fait une véritable moyenne arithmétique des angles.

Tout d'abord, nous devons examiner dans quels cas exactement la différence entre les angles est différente de la différence entre leurs homologues en nombres normaux. Considérons les angles x et y, si y >= x - 180 et y <= x + 180, alors nous pouvons utiliser directement la différence (x-y). Sinon, si la première condition n'est pas remplie, nous devons utiliser (y+360) dans le calcul au lieu de y. De même, si la deuxième condition n'est pas remplie, nous devons utiliser (y-360) au lieu de y. Étant donné que l'équation de la courbe que nous minimisons ne change qu'aux points où ces inégalités passent de vrai à faux ou vice versa, nous pouvons séparer la plage complète [0,360] en un ensemble de segments, séparés par ces points. Il nous suffit alors de trouver le minimum de chacun de ces segments, puis le minimum du minimum de chaque segment, qui est la moyenne.

Voici une image qui montre où se situent les problèmes de calcul des différences d'angle. Si x se trouve dans la zone grise, il y aura un problème.

Pour minimiser une variable, selon la courbe, on peut prendre la dérivée de ce que l'on veut minimiser et ensuite on trouve le point d'inflexion (qui est celui où la dérivée = 0).

Nous appliquerons ici l'idée de minimiser la différence au carré pour dériver la formule commune de la moyenne arithmétique : somme(a[i])/n. La courbe y = somme((a[i]-x)^2) peut être minimisée de cette manière :

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

Maintenant, on l'applique aux courbes avec nos différences ajustées :

b = sous-ensemble de a où la différence (angulaire) correcte a[i]-x c = sous-ensemble de a où la différence (angulaire) correcte (a[i]-360)-x cn = taille de c d = sous-ensemble de a où la différence (angulaire) correcte (a[i]+360)-x dn = taille 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

Cela ne suffit pas pour obtenir le minimum, alors que cela fonctionne pour les valeurs normales, qui ont un ensemble non limité, donc le résultat sera certainement dans la gamme de l'ensemble et est donc valide. Nous avons besoin du minimum dans un intervalle (défini par le segment). Si le minimum est inférieur à la limite inférieure de notre segment, alors le minimum de ce segment doit être à la limite inférieure (car les courbes quadratiques n'ont qu'un seul point de retournement) et si le minimum est supérieur à la limite supérieure de notre segment, alors le minimum du segment est à la limite supérieure. Une fois que nous avons le minimum de chaque segment, nous trouvons simplement celui qui a la valeur la plus basse pour ce que nous minimisons (somme((b[i]-x)^2) + somme(((c[i]-360)-b)^2) + somme(((d[i]+360)-c)^2)).

Voici une image de la courbe, qui montre comment elle change aux points où x=(a[i]+180)%360. L'ensemble de données en question est {65,92,230,320,250}.

Voici une implémentation de l'algorithme en Java, incluant quelques optimisations, sa complexité est O(nlogn). Elle peut être réduite à O(n) si vous remplacez le tri basé sur la comparaison par un tri non basé sur la comparaison, tel que le tri 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_;
}

La moyenne arithmétique d'un ensemble d'angles peut ne pas correspondre à votre idée intuitive de ce que devrait être la moyenne. Par exemple, la moyenne arithmétique de l'ensemble {179,179,0,181,181} est 216 (et 144). La réponse à laquelle vous pensez immédiatement est probablement 180, mais il est bien connu que la moyenne arithmétique est fortement affectée par les valeurs des angles. Vous devez également vous rappeler que les angles ne sont pas des vecteurs, aussi séduisant que cela puisse paraître lorsqu'on traite parfois des angles.

Cet algorithme s'applique bien sûr aussi à toutes les quantités qui obéissent à l'arithmétique modulaire (avec un ajustement minimal), comme l'heure du jour.

Je tiens également à souligner que même s'il s'agit d'une véritable moyenne d'angles, contrairement aux solutions vectorielles, cela ne signifie pas nécessairement que c'est la solution que vous devez utiliser, la moyenne des vecteurs unitaires correspondants peut très bien être la valeur que vous devez utiliser.