Ответ 1
Вычислить единичные векторы из углов и взять угол их среднего.
Я хочу рассчитать среднее значение набора круговых данных. Например, у меня может быть несколько образцов из чтения компаса. Проблема, конечно, в том, как бороться с обманом. Тот же алгоритм может быть полезен для часов.
Актуальный вопрос более сложный - что означают статистические данные о сфере или в алгебраическом пространстве, которое "обертывается", например. аддитивная группа mod n. Ответ может быть не единственным, например. среднее значение 359 градусов и 1 градуса может быть 0 градусов или 180, но статистически 0 выглядит лучше.
Это реальная проблема программирования для меня, и я пытаюсь сделать ее не похожей на математическую проблему.
Вычислить единичные векторы из углов и взять угол их среднего.
Этот вопрос подробно рассматривается в книге: "Статистика по сферам", Джеффри С. Уотсон, Университет Арканзаса Лекция Примечания в математических науках, 1983 John Wiley and Sons, Inc., упомянутые в http://catless.ncl.ac.uk/Risks/7.44.html#subj4 Брюса Карша.
Хороший способ оценить средний угол A от набора угловых измерений a [i] 0 <= i
sum_i_from_1_to_N sin(a[i])
a = arctangent ---------------------------
sum_i_from_1_to_N cos(a[i])
Метод, заданный starblue, является эквивалентным по вычислительной мере, но его причины более ясны и, вероятно, программно более эффективны, а также хорошо работают в нулевом случае, так что он ему славен.
Теперь предмет более подробно рассматривается в Википедии, а также в других целях, таких как дробные части.
Я вижу проблему - например, если у вас есть угол 45 'и угол 315', "естественное" среднее значение будет равно 180, но значение, которое вы хотите, на самом деле равно 0 ".
Я думаю, что Starblue на что-то. Просто вычислите (x, y) декартовы координаты для каждого угла и добавьте эти полученные векторы вместе. Смещение angular конечного вектора должно быть вашим требуемым результатом.
x = y = 0
foreach angle {
x += cos(angle)
y += sin(angle)
}
average_angle = atan2(y, x)
Я сейчас игнорирую, что курсор компаса начинается с севера и идет по часовой стрелке, тогда как "нормальные" декартовы координаты начинаются с нуля вдоль оси X, а затем идут против часовой стрелки. Математика должна работать одинаково независимо.
ДЛЯ СПЕЦИАЛЬНОГО ДЕЛА ДВА УГЛЫ:
Ответ ((a + b) mod 360)/2 WRONG. Для углов 350 и 2 ближайшая точка равна 356, а не 176.
Единичный вектор и триггерные решения могут быть слишком дорогими.
То, что я получаю от небольшого мастерства, это:
diff = ( ( a - b + 180 + 360 ) mod 360 ) - 180
angle = (360 + b + ( diff / 2 ) ) mod 360
ackb прав, что эти векторные решения не могут считаться истинными средними по углам, они являются лишь средним числом единичных векторов. Однако предлагаемое решение ackb не кажется математически звуковым.
Следующее - это решение, которое математически выводится из цели минимизации (angle [i] - avgAngle) ^ 2 (где при необходимости исправляется разница), что делает его истинным средним арифметическим углов.
Во-первых, нам нужно посмотреть, в каких случаях разница между углами отличается от разницы между их нормальными номерами. Рассмотрим углы x и y, если y >= x - 180 и y <= x + 180, то мы можем непосредственно использовать разность (x-y). В противном случае, если первое условие не выполняется, мы должны использовать (y + 360) в вычислении вместо y. Соответственно, если второе условие не выполняется, мы должны использовать (y-360) вместо y. Поскольку уравнение кривой мы минимизируем только изменения в тех точках, где эти неравенства изменяются от истинных до ложных или наоборот, мы можем разделить полный диапазон [0,360) на множество сегментов, разделенных этими точками. Тогда нам нужно всего лишь найти минимум каждого из этих сегментов, а затем минимум минимума каждого сегмента, который является средним.
Здесь показано изображение, где возникают проблемы при расчете разностей углов. Если x лежит в серой области, тогда будет проблема.
Чтобы свести к минимуму переменную, зависящую от кривой, мы можем взять производную от того, что мы хотим минимизировать, а затем найдем точку поворота (где производная = 0).
Здесь мы применим идею минимизации квадрата разности для получения общей средней арифметической формулы: sum (a [i])/n. Кривая y = sum ((a [i] -x) ^ 2) может быть минимизирована таким образом:
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
Теперь применим его к кривым с нашими скорректированными разностями:
b = подмножество a, где правильная (angular) разность a [i] -x c = подмножество a, где правильная (angular) разность (a [i] -360) -x cn = размер c d = подмножество a, где правильная (angular) разность (a [i] +360) -x dn = размер 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
Этого недостаточно, чтобы получить минимум, в то время как он работает для нормальных значений, у которых есть неограниченный набор, поэтому результат определенно будет находиться в пределах заданного диапазона и, следовательно, действителен. Нам нужен минимум в пределах диапазона (определенного сегментом). Если минимум меньше нижней границы нашего сегмента, то минимум этого сегмента должен быть на нижней границе (поскольку квадратичные кривые имеют только одну точку поворота), и если минимум больше верхней границы сегмента, то минимальный отрезок находится на верхняя граница. После того, как у нас есть минимум для каждого сегмента, мы просто найдем тот, который имеет самое низкое значение для того, что мы минимизируем (sum ((b [i] -x) ^ 2) + sum (((c [i] -360 ) -b) ^ 2) + sum (((d [i] +360) -c) ^ 2)).
Вот изображение кривой, которое показывает, как оно изменяется в точках, где x = (a [i] +180)% 360. Набор данных находится под вопросом: {65,92,230,320,250}.
Вот реализация алгоритма в Java, включая некоторые оптимизации, его сложность - O (nlogn). Его можно свести к O (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_;
}
Среднее арифметическое набора углов может не совпадать с вашим интуитивным представлением о том, что должно быть в среднем. Например, среднее арифметическое набора {179,179,0,181,181} составляет 216 (и 144). Ответ, о котором вы сразу думаете, вероятно, равен 180, однако хорошо известно, что среднее арифметическое сильно зависит от значений краев. Вы также должны помнить, что углы не являются векторами, как это может показаться, когда они иногда имеют дело с углами.
Этот алгоритм, конечно, также применим ко всем величинам, которые подчиняются модулярной арифметике (с минимальной корректировкой), например, времени суток.
Я также хотел бы подчеркнуть, что, хотя это истинное среднее значение углов, в отличие от векторных решений, это не обязательно означает, что это решение, которое вы должны использовать, среднее значение соответствующих единичных векторов вполне может быть значение, которое вы действительно должны использовать.
Вы должны определить среднее значение более точно. Для конкретного случая с двумя углами я могу представить два разных сценария:
Я не вижу, как вторая альтернатива может быть обобщена для случая с более чем двумя углами.
Как и все средние значения, ответ зависит от выбора метрики. Для данной метрики M среднее значение некоторых углов a_k в [-pi, pi] для k в [1, N] - это угол a_M, который минимизирует сумму квадратов расстояний d ^ 2_M (a_M, a_k). Для взвешенного среднего просто включается в сумму веса w_k (такие, что sum_k w_k = 1). То есть
a_M = arg min_x sum_k w_k d ^ 2_M (x, a_k)
Два общих варианта метрики - это метрики Фробениуса и Римана. Для метрики Фробениуса существует прямая формула, соответствующая обычному понятию среднего значения в круговой статистике. Подробнее см. "Средства и усреднение в группе ротации", Махер Моакер, журнал SIAM по анализу и применению матриц, том 24, выпуск 1, 2002.
http://link.aip.org/link/?SJMAEL/24/1/1
Здесь функция GNU Octave 3.2.4, которая вычисляет:
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/'.
Я хотел бы поделиться методом, который я использовал с микроконтроллером, который не обладал возможностями с плавающей точкой или тригонометрикой. Я до сих пор нужно "средних" 10 сырых показаний подшипников для того, чтобы сгладить вариации.
Это не идеально; он может сломаться. В этом случае я сошел с рук, потому что устройство вращается очень медленно. Я поставлю его там, если кто-то другой найдет себя работающим с аналогичными ограничениями.
Я бы использовал векторный путь, используя сложные числа. Мой пример находится в Python, который имеет встроенные комплексные номера:
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)
Обратите внимание, что Python не нужно создавать временный новый список векторов, все это можно сделать за один шаг; Я просто выбрал этот способ, чтобы аппроксимировать псевдокод, применимый и к другим языкам.
Вот полное решение: (вход представляет собой массив подшипников в градусах (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;
}
Здесь полное решение на С++:
#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);
}
Он принимает углы в виде вектора удвоений и возвращает среднее просто как двойное. Углы должны быть в градусах, и, конечно, среднее значение также находится в градусах.
Здесь идея: постройте среднее значение итеративно, всегда вычисляя среднее значение углов, которые находятся ближе всего друг к другу, сохраняя вес.
Другая идея: найти наибольший разрыв между заданными углами. Найти точку, разделяющую его, а затем выбрать противоположную точку на окружности в качестве опорного нуля, чтобы вычислить среднее значение из.
Представим эти углы с точками на окружности круга.
Можно ли считать, что все эти точки падают на одну и ту же половину круга? (В противном случае нет очевидного способа определить "средний угол". Подумайте о двух точках диаметра, например, 0 град и 180 градусов --- это средний 90 градусов или 270 градусов? Что происходит, когда у нас есть 3 или более равномерно распределенные точки?)
В этом предположении мы выбираем произвольную точку полуокружности как "начало" и измеряем заданный набор углов относительно этого начала (назовем это "относительным углом" ). Обратите внимание, что относительный угол имеет абсолютное значение строго менее 180 град. Наконец, возьмите среднее значение этих относительных углов, чтобы получить желаемый средний угол (относительно нашего происхождения).
Нет единого "правильного ответа". Я рекомендую прочитать книгу, К. В. Мардия и П. Е. Юпп, "Направленная статистика", (Wiley, 1999), для тщательного анализа.
Alnitak имеет правильное решение. Решение Nic Fortescue функционально одинаково.
Для частного случая где
(sum (x_component) = 0.0 & sum (y_component) = 0.0)//например. 2 угла 10 и 190 градусов.
используйте 0,0 градуса в виде суммы
Вычислительно вы должны проверить этот случай, так как atan2 (0., 0) undefined и генерирует ошибку.
Средний угол phi_avg должен обладать тем свойством, что sum_i | phi_avg-phi_i | ^ 2 становится минимальным, где разница должна быть в [-Pi, Pi) (потому что это может быть короче, чтобы идти наоборот!), Это легко достигается путем нормализации всех входных значений до [0, 2Pi), сохраняя текущее среднее значение phi_run и выбрав нормализацию | phi_i-phi_run | к [-Pi, Pi) (путем добавления или вычитания 2Pi). Большинство предложений выше делают что-то еще, что не имеют такое минимальное свойство, т.е. что-то среднее, но не углы.
По-английски:
В python:
A # массивный массив NX1 углов
if np.var(A) < np.var((A-180)%360):
average = np.average(A)
else:
average = (np.average((A-180)%360)+180)%360
(Просто хочу поделиться своей точкой зрения с теорией оценки или статистическим выводом)
Проворное испытание - получить оценку MMSE ^ набора углов, но это один из вариантов выбора "усредненного" направления; можно также найти оценку MMAE ^ или некоторую другую оценку как "усредненное" направление, и это зависит от вашей метрики, определяющей ошибку направления; или в целом в теории оценки, определение функции стоимости.
^ MMSE/MMAE соответствует минимальной среднеквадратичной/абсолютной ошибке.
ackb сказал: "Средний угол phi_avg должен обладать свойством, что sum_i | phi_avg-phi_i | ^ 2 становится минимальным... они усредняют что-то, но не имеют углы"
---- вы количественно оцениваете ошибки в среднеквадратичном смысле, и это один из наиболее распространенных способов, однако, не единственный способ. Ответ, которым пользуются большинство людей здесь (т.е. Сумма единичных векторов и получение угла результата), на самом деле является одним из разумных решений. Это (может быть доказано) оценка ML, которая служит "усредненным" направлением, которое мы хотим, если направления векторов моделируются как распределение фон Мизеса. Это распределение не причудливо и представляет собой просто периодическое выборочное распределение от двумерного гуассиана. См. (2.179) в епископской книге "Распознавание образов и машинное обучение". Опять же, ни в коем случае это единственный лучший, чтобы представлять "среднее" направление, однако, это вполне разумный, который имеет как хорошее теоретическое обоснование, так и простую реализацию.
Nimble сказал, что "ackb прав, что эти векторные решения не могут считаться истинными средними по углам, они являются лишь средним числом единичных векторных копий"
---- это не так. "Единичные векторные сопоставления" раскрывают информацию о направлении вектора. Угол - это величина без учета длины вектора, а единичный вектор - с дополнительной информацией о том, что длина равна 1. Вы можете определить, что ваш "единичный" вектор имеет длину 2, это не имеет особого значения.
Я решил проблему с помощью ответа от @David_Hanak. Как он утверждает:
Угол, который указывает "между" двумя другими, оставаясь в том же полукруге, например. для 355 и 5 это будет 0, а не 180. Для этого вам нужно проверить, превышает ли разность между двумя углами 180 или нет. Если это так, увеличьте меньший угол на 360, прежде чем использовать приведенную выше формулу.
Так что я сделал, вычислил среднее значение всех углов. А затем все углы, которые меньше этого, увеличат их на 360. Затем пересчитайте среднее значение, добавив их все и разделив их по длине.
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;
Работает отлично.
Функция 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
Вы можете использовать эту функцию в 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
Вы можете увидеть решение и небольшое объяснение в следующей ссылке, для ЛЮБОГО языка программирования: https://rosettacode.org/wiki/Averages/Mean_angle
Например, С++-решение:
#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;
}
Вывод:
Mean Angle for 1st set : -0.000000 degrees
Mean Angle for 2nd set : -90.000000 degrees
Mean Angle for 3rd set : 20.000000 degrees
Или Решение 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
В то время как ответ звездной звезды дает угол среднего единичного вектора, можно расширить понятие среднего арифметического на углы, если вы согласны с тем, что может быть более одного ответа в диапазоне от 0 до 2 * pi (или От 0 ° до 360 °). Например, среднее значение 0 ° и 180 ° может составлять либо 90 °, либо 270 °.
Среднее арифметическое имеет свойство быть единственным значением с минимальной суммой квадратов расстояний до входных значений. Расстояние по единичному кругу между двумя единичными векторами можно легко вычислить как обратный косинус их точечного произведения. Если мы выберем единичный вектор, минимизируя сумму квадрата обратного косинуса точечного произведения нашего вектора и каждого входного единичного вектора, то мы имеем эквивалентное среднее. Опять же, имейте в виду, что в исключительных случаях могут быть два или более минимума.
Эта концепция может быть расширена до любого числа измерений, так как расстояние вдоль единичной сферы можно рассчитать точно так же, как расстояние по единице окружности - обратный косинус точечного произведения двух единичных векторов.
Для кругов мы могли бы решить для этого среднего несколькими способами, но я предлагаю следующий алгоритм O (n ^ 2) (углы находятся в радианах, и я не могу вычислить единичные векторы):
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
Если все углы находятся в пределах 180 ° друг от друга, мы могли бы использовать более простой алгоритм O (n) + O (сортировка) (снова используя радианы и избегая использования единичных векторов):
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
Чтобы использовать градусы, просто замените pi на 180. Если вы планируете использовать больше измерений, вам, скорее всего, придется использовать итеративный метод для решения для среднего.
Вот полностью арифметическое решение, использующее скользящие средние и заботящееся о нормализации значений. Он быстрый и дает правильные ответы, если все углы находятся на одной стороне круга (в пределах 180 ° друг от друга).
Это математически эквивалентно добавлению смещения, которое сдвигает значения в диапазон (0, 180), вычисляя среднее значение, а затем вычитая смещение.
В комментариях описывается, какой диапазон может принимать определенное значение в любой момент времени
// 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;
}
На основе ответа Alnitak я написал метод Java для вычисления среднего значения нескольких углов:
Если ваши углы находятся в радианах:
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);
}
Если ваши углы находятся в градусах:
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));
}
Проблема чрезвычайно проста. 1. Убедитесь, что все углы находятся между -180 и 180 градусов. 2. a Добавьте все неотрицательные углы, возьмите их среднее значение и COUNT сколько 2. b.Добавьте все отрицательные углы, возьмите их среднее значение и COUNT сколько. 3. Возьмите разницу pos_average минус neg_average Если разница больше 180, измените разницу на 360 минус разницу. В противном случае просто измените знак различия. Обратите внимание, что разница всегда неотрицательна. Среднее_Англе равно разнице pos_average plus раз "вес", отрицательное количество, деленное на сумму отрицательного и положительного count
В среднем два угла есть два средних значения на 180 ° друг от друга, но нам может понадобиться более близкое среднее.
Визуально среднее значение синего (b) и зеленого ( a) дает точку чика:
Обтекание углов вокруг (например, 355 + 10 = 5), но стандартная арифметика будет игнорировать эту точку ветвления. Однако, если угол b противоположный точке ветвления, то (b + g)/2 дает самое близкое среднее: точка чика.
Для любых двух углов мы можем вращать проблему, чтобы один из углов был противоположным точке ветвления, выполнял стандартное усреднение, а затем вращался назад.
В python с углами между [-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)
У меня есть другой метод, чем @Starblue, который дает "правильные" ответы на некоторые из приведенных выше углов. Например:
Он использует сумму по разности между последовательными углами. Код (в 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