¿Cómo se calcula el promedio de un conjunto de datos circulares?

Esta pregunta se examina en detalle en el libro: "Estadísticas en esferas", Geoffrey S. Watson, Universidad de Arkansas Lecciones Notas en las Ciencias Matemáticas, 1983 John Wiley & Sons, Inc. como se menciona en http://catless.ncl.ac.uk/Risks/7.44.html#subj4 por Bruce Karsh.

Una buena manera de estimar un ángulo promedio, A, a partir de un conjunto de medidas de ángulos a[i] 0<=i

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

El método dado por starblue es computacionalmente equivalente, pero sus razones son más claras y probablemente más eficientes programáticamente, y también funcionan bien en el caso cero, así que felicitaciones para él.

El tema se explora ahora con más detalle en Wikipedia, junto con otros usos, como partes fraccionarias.

Veo el problema: por ejemplo, si tienes un ángulo de 45° y un ángulo de 315°, el promedio "natural" sería 180°, pero el valor que quieres es en realidad 0°.

Creo que Starblue está en algo. Simplemente calcula las coordenadas cartesianas (x, y) para cada ángulo y suma esos vectores resultantes juntos. El desplazamiento angular del vector final debería ser tu resultado requerido.

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

Por ahora estoy ignorando que un rumbo de una brújula comienza en el norte y va en sentido horario, mientras que las coordenadas cartesianas "normales" comienzan con cero a lo largo del eje X y luego van en sentido contrario a las agujas del reloj. Las matemáticas deberían funcionar de la misma manera independientemente.

PARA EL CASO ESPECIAL DE DOS ÁNGULOS:

La respuesta ( (a + b) mod 360 ) / 2 está INCORRECTA. Para los ángulos 350 y 2, el punto más cercano es 356, no 176.

Las soluciones de vector unitario y trigonométricas pueden ser demasiado costosas.

Lo que he obtenido de un poco de experimentación es:

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

• 0, 180 -> 90 (dos respuestas para esto: esta ecuación toma la respuesta en sentido horario desde a)
• 180, 0 -> 270 (ver arriba)
• 180, 1 -> 90.5
• 1, 180 -> 90.5
• 20, 350 -> 5
• 350, 20 -> 5 (todos los ejemplos siguientes se invierten correctamente también)
• 10, 20 -> 15
• 350, 2 -> 356
• 359, 0 -> 359.5
• 180, 180 -> 180

ackb tiene razón en que estas soluciones basadas en vectores no pueden considerarse verdaderos promedios de ángulos, sino que son sólo un promedio de las contrapartes vectoriales unitarias. Sin embargo, la solución sugerida por ackb no parece tener sentido matemáticamente.

La siguiente es una solución que se deriva matemáticamente de la meta de minimizar (ángulo[i] - avgAngulo)^2 (donde la diferencia se corrige si es necesario), lo que la convierte en una verdadera media aritmética de los ángulos.

En primer lugar, tenemos que mirar exactamente en qué casos la diferencia entre los ángulos es diferente a la diferencia entre sus contrapartes de número normal. Consideremos los ángulos x y y, si y >= x - 180 y y <= x + 180, entonces podemos usar la diferencia (x-y) directamente. De lo contrario, si la primera condición no se cumple, entonces debemos usar (y+360) en el cálculo en lugar de y. Correspondientemente, si la segunda condición no se cumple, entonces debemos usar (y-360) en lugar de y. Como la ecuación de la curva que estamos minimizando sólo cambia en los puntos donde estas desigualdades cambian de verdadero a falso o viceversa, podemos separar el rango completo [0,360) en un conjunto de segmentos, separados por estos puntos. Entonces, sólo necesitamos encontrar el mínimo de cada uno de estos segmentos, y luego el mínimo del mínimo de cada segmento, que es el promedio.

Aquí hay una imagen que demuestra dónde se producen los problemas en el cálculo de las diferencias de ángulo. Si x se encuentra en el área gris, entonces habrá un problema.

http://i3.photobucket.com/albums/y69/02williamsa6/angle_comparisons.png

Para minimizar una variable, dependiendo de la curva, podemos tomar la derivada de lo que queremos minimizar y luego encontramos el punto de giro (que es donde la derivada = 0).

Aquí aplicaremos la idea de minimizar la diferencia al cuadrado para derivar la fórmula de la media aritmética común: suma(a[i])/n. La curva y = suma((a[i]-x)^2) puede ser minimizada de esta manera:

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

Ahora lo aplicamos a las curvas con nuestras diferencias ajustadas:

b = subconjunto de a donde la diferencia (angular) correcta a[i]-x c = subconjunto de a donde la diferencia (angular) correcta (a[i]-360)-x cn = tamaño de c d = subconjunto de a donde la diferencia (angular) correcta (a[i]+360)-x dn = tamaño 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

Esto por sí solo no es suficiente para obtener el mínimo, mientras que funciona para los valores normales, que tiene un conjunto sin límites, por lo que el resultado se encontrará definitivamente dentro del rango del conjunto y por lo tanto es válido. Necesitamos el mínimo dentro de un rango (definido por el segmento). Si el mínimo es menor que el límite inferior de nuestro segmento, entonces el mínimo de ese segmento debe estar en el límite inferior (porque las curvas cuadráticas sólo tienen un punto de giro) y si el mínimo es mayor que el límite superior de nuestro segmento, entonces el mínimo del segmento está en el límite superior. Después de tener el mínimo de cada segmento, simplemente encontramos el que tiene el valor más bajo de lo que estamos minimizando (suma((b[i]-x)^2) + suma(((c[i]-360)-b)^2) + suma(((d[i]+360)-c)^2)).

Aquí hay una imagen de la curva, que muestra cómo cambia en los puntos donde x=(a[i]+180)%360. El conjunto de datos en cuestión es {65,92,230,320,250}.

http://i3.photobucket.com/albums/y69/02williamsa6/curve.png

Aquí hay una implementación del algoritmo en Java, incluyendo algunas optimizaciones, su complejidad es O(nlogn). Puede reducirse a O(n) si se reemplaza la clasificación basada en la comparación con una clasificación no basada en la comparación, como la clasificación 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 media aritmética de un conjunto de ángulos puede no estar de acuerdo con su idea intuitiva de lo que debería ser la media. Por ejemplo, la media aritmética del conjunto {179,179,0,181,181} es 216 (y 144). La respuesta en la que piensas inmediatamente es probablemente 180, sin embargo es bien sabido que la media aritmética está fuertemente afectada por los valores de los bordes. También debes recordar que los ángulos no son vectores, por muy atractivo que parezca a veces cuando se trata de ángulos.

Por supuesto, este algoritmo también se aplica a todas las cantidades que obedecen a la aritmética modular (con un ajuste mínimo), como la hora del día.

También quisiera subrayar que, aunque se trata de una verdadera media de ángulos, a diferencia de las soluciones vectoriales, eso no significa necesariamente que sea la solución que debe utilizarse, la media de los vectores unitarios correspondientes bien puede ser el valor que realmente debe utilizarse.