programing tip

순환 데이터 세트의 평균을 어떻게 계산합니까?

itbloger 2020. 6. 21. 20:12
반응형

순환 데이터 세트의 평균을 어떻게 계산합니까?


순환 데이터 세트의 평균을 계산하고 싶습니다. 예를 들어 나침반을 읽은 결과가 몇 가지있을 수 있습니다. 물론 문제는 랩 어라운드를 처리하는 방법입니다. 동일한 알고리즘이 문자판에 유용 할 수 있습니다.

실제 질문은 더 복잡합니다. 통계는 구 또는 대수 공간에서 "주변"(예 : 첨가제 그룹 mod n)에서 무엇을 의미합니까? 답은 고유하지 않을 수 있습니다. 예를 들어 평균 359도, 1 도는 0도 또는 180도이지만 통계적으로 0이 더 좋습니다.

이것은 실제 프로그래밍 문제이며 수학 문제처럼 보이지 않도록 노력하고 있습니다.


각도에서 단위 벡터를 계산하고 평균 각도를 가져옵니다.


이 질문은 "Statistics On Spheres", 아칸소 대학교 게 프리 S. 왓슨, 수학 과학 강의 노트, 1983 John Wiley & Sons, Inc. ( http : //catless.ncl. Bruce Karsh의 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가 제공하는 방법은 계산 상 동일하지만 그의 이유는 명확하고 프로그래밍 방식으로 더 효율적이며 제로의 경우에도 잘 작동하므로 그에게 큰 도움이됩니다.

이 주제는 이제 Wikipedia 및 분수 부분과 같은 다른 용도 로 더 자세히 탐색됩니다 .


예를 들어 45 '각도와 315'각도를 가진 경우 "자연"평균은 180 '이지만 원하는 값은 실제로 0'입니다.

스타 블루가 뭔가 있다고 생각합니다. 각 각도의 (x, y) 직교 좌표를 계산하고 그 결과 벡터를 함께 추가하십시오. 최종 벡터의 각도 오프셋은 필요한 결과 여야합니다.

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

나침반 머리글이 북쪽에서 시작하여 시계 방향으로 이동하는 반면 "일반"직교 좌표는 X 축을 따라 0으로 시작한 다음 시계 반대 방향으로 이동한다는 점을 무시하고 있습니다. 수학은 상관없이 같은 방식으로 작동해야합니다.


두 각도의 특수한 경우 :

((a + b) mod 360) / 2WRONG 입니다. 각도 350과 2의 경우 가장 가까운 점은 176이 아니라 356입니다.

단위 벡터 및 삼각 솔루션이 너무 비쌀 수 있습니다.

내가 약간 어설프게 얻은 것은 다음과 같습니다.

diff = ( ( a - b + 180 + 360 ) mod 360 ) - 180
angle = (360 + b + ( diff / 2 ) ) mod 360
  • 0, 180-> 90 (이에 대한 두 가지 답변 :이 방정식은 a에서 시계 방향으로 응답합니다)
  • 180, 0-> 270 (위 참조)
  • 180, 1-> 90.5
  • 1, 180-> 90.5
  • 20, 350-> 5
  • 350, 20-> 5 (다음 예도 모두 올바르게 반전 됨)
  • 10, 20-> 15
  • 350, 2-> 356
  • 359, 0-> 359.5
  • 180, 180-> 180

ackb는 이러한 벡터 기반 솔루션이 실제 각도의 평균으로 간주 될 수 없으며 단위 벡터 대응의 평균 일 뿐이라는 점이 맞습니다. 그러나 ackb의 제안 된 솔루션은 수학적으로 들리는 것처럼 보이지 않습니다.

다음은 (angle [i]-avgAngle) ^ 2 (필요한 경우 차이가 수정되는 경우) 최소화라는 목표에서 수학적으로 도출 된 솔루션으로 각도의 실제 산술 평균이됩니다.

먼저 각도의 차이가 정상적인 수의 차이와 다른 경우를 정확히 찾아야합니다. y> = x-180 및 y <= x + 180 인 경우 각도 x와 y를 고려하면 차이 (xy)를 직접 사용할 수 있습니다. 그렇지 않으면 첫 번째 조건이 충족되지 않으면 y 대신 (y + 360)을 계산에 사용해야합니다. 해당하는 두 번째 조건이 충족되지 않으면 y 대신 (y-360)을 사용해야합니다. 곡선의 방정식은 이러한 불평등이 참에서 거짓으로 또는 그 반대로 변하는 지점에서만 변화를 최소화하기 때문에 전체 [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 [i] -xc = a의 부분 집합 정확한 (각도) 차이 (a [i] -360) -x cn = cd의 크기 = 정확한 (각도) 차이 (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

이것만으로는 최소값을 얻는 데 충분하지 않으며, 한정되지 않은 세트를 가진 정상 값에 대해 작동하므로 결과는 세트의 범위 내에 있으므로 유효합니다. 범위 내에서 최소값이 필요합니다 (세그먼트로 정의). 최소값이 세그먼트의 하한보다 작 으면 해당 세그먼트의 최소값이 하한에 있어야합니다 (2 차 곡선에 1 개의 회전 점이 있기 때문에) 최소값이 세그먼트의 상한값보다 큰 경우 세그먼트의 최소값은 상한. 각 세그먼트에 대한 최소값을 얻은 후 최소화하려는 값이 가장 낮은 세그먼트를 찾으면됩니다 (sum ((b [i] -x) ^ 2) + sum ((((c [i] -360)) ) -b) ^ 2) + 합 (((d [i] +360) -c) ^ 2)).

다음은 곡선에 대한 이미지입니다. x = (a [i] +180) % 360 지점에서 어떻게 변하는지를 보여줍니다. 해당 데이터 세트는 {65,92,230,320,250}입니다.

곡선

다음은 일부 최적화를 포함하여 Java로 알고리즘을 구현 한 것으로 복잡도는 O (nlogn)입니다. 비교 기반 정렬을 기수 정렬과 같은 비 비교 기반 정렬로 바꾸면 O (n)으로 줄일 수 있습니다.

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 일 것입니다. 그러나 산술 평균은 엣지 값에 큰 영향을받는 것으로 잘 알려져 있습니다. 또한 각도가 벡터가 아니라는 점을 기억해야합니다. 때로는 각도를 처리 할 때 보이는 것처럼 매력적입니다.

이 알고리즘은 물론 시간과 같이 모듈 식 산술 (최소 조정)을 준수하는 모든 수량에도 적용됩니다.

또한 이것이 벡터 솔루션과 달리 이것이 실제 각도의 평균 임에도 불구하고 반드시 사용해야하는 솔루션이라는 의미는 아니지만 해당 단위 벡터의 평균이 실제 값일 수 있음을 강조하고 싶습니다. 사용해야합니다.


평균을 보다 정확하게 정의해야 합니다. 두 각도의 특정 경우에 대해 두 가지 시나리오를 생각할 수 있습니다.

  1. "참"평균, 즉 (a + b) / 2 % 360.
  2. 같은 반원을 유지하면서 다른 두 사람 사이를 "예"(예 : 355 및 5)하는 각도는 180이 아닌 0입니다. 이렇게하려면 두 각도의 차이가 180보다 큰지 확인해야합니다. 또는 아닙니다. 그렇다면 위의 공식을 사용하기 전에 작은 각도를 360 씩 늘리십시오.

그래도 두 개 이상의 각도의 경우 두 번째 대안이 어떻게 일반화 될 수 있는지 알지 못합니다.


모든 평균과 마찬가지로 답은 측정 항목 선택에 따라 다릅니다. 주어진 메트릭 M에 대해, [1, N]의 k에 대한 [-pi, pi]의 일부 각도 a_k의 평균은 제곱 거리 d ^ 2_M (a_M, a_k)의 합을 최소화하는 각도 a_M입니다. 가중 평균의 경우, 단순히 가중치 w_k (sum_k w_k = 1)를 합계에 포함합니다. 그건,

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

두 가지 일반적인 메트릭 선택은 Frobenius 및 Riemann 메트릭입니다. Frobenius 지표의 경우 순환 통계에서 일반적인 평균 방위 개념에 해당하는 직접 공식이 존재합니다. 자세한 내용은 "회전 그룹의 평균 및 평균화", Maher Moakher, SIAM Journal of Matrix Analysis and Applications, Volume 24, Issue 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, robjohnson@alphawaveresearch.com
%     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 개의 원시 베어링 판독 값을 "평균화"해야했습니다.

  1. 첫 번째 베어링의 범위가 270-360 또는 0-90 도인 지 확인하십시오 (2 사분면).
  2. 만약 그렇다면이 값과 이후의 모든 판독 값을 180도 회전시켜 0 <= bearing <360 범위의 모든 값을 유지하십시오. 그렇지 않으면 판독 값을 그대로 가져갑니다.
  3. 10 회 판독 된 후 랩 어라운드가 없다고 가정하고 수치 평균을 계산하십시오.
  4. 180도 회전이 적용된 경우 계산 된 평균을 180도 회전하여 "진정한"베어링으로 ​​돌아갑니다.

이상적이지 않습니다. 깨질 수 있습니다. 이 경우 장치가 매우 느리게 회전하기 때문에이 문제를 해결했습니다. 다른 사람이 비슷한 제한을 받고 일하는 것을 발견 할 경우를 대비해 설명하겠습니다.


전체 솔루션은 다음과 같습니다. (입력은 각도 배열 (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;
}

파이썬에서 [-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)

세부:

의 평균 이 각도 가이 명 평균 180 ° 떨어져있다, 그러나 우리는 가까운 평균을 할 수 있습니다.

시각적으로, 청색 (평균 B ) 및 녹색 ( )를 맺는 틸 점 :

실물

각도는 '주변'(예 : 355 + 10 = 5)이지만 표준 산술에서는이 분기점을 무시합니다. 그러나 각도 b 가 분기점과 반대이면 ( b + g ) / 2가 가장 가까운 평균 인 청록 점을 제공합니다.

두 각도의 경우 문제를 회전시켜 각도 중 하나가 분기점과 반대가되도록 표준 평균화를 수행 한 다음 다시 회전 할 수 있습니다.

회전반환


나는 복잡한 숫자를 사용하여 벡터 길을 갈 것입니다. 내 예제는 복잡한 숫자가 내장 된 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)

파이썬은 임시로 새로운 벡터 목록을 만들 필요 가 없습니다 . 위의 모든 것을 한 단계로 수행 할 수 있습니다. 방금 다른 언어에도 적용 가능한 의사 코드를 추정하기 위해이 방법을 선택했습니다.


완전한 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);
}

복소수 벡터의 형태로 각도를 취하고 평균을 단순히 복소수로 반환합니다. 각도는도 단위 여야하며 물론 평균도 도입니다.


아이디어는 다음과 같습니다. 가중치를 유지하면서 서로 가장 가까운 각도의 평균을 항상 계산하여 평균을 반복적으로 만듭니다.

또 다른 아이디어 : 주어진 각도 사이에서 가장 큰 간격을 찾으십시오. 이등분하는 점을 찾은 다음 원의 반대쪽 점을 참조 0으로 선택하여 평균을 계산합니다.


원의 원주에있는 점으로이 각도를 표현해 봅시다.

이 모든 점들이 원의 같은 절반에 있다고 가정 할 수 있습니까? (그렇지 않으면, "평균 각도"를 정의하는 명백한 방법은 없습니다. 예를 들어 0도 및 180도에 대해 직경의 두 점을 생각해보십시오. 평균 90도 또는 270 도는 무엇입니까? 골고루 퍼지나요?)

이 가정을 통해 반원에서 임의의 점을 "원점"으로 선택하고이 원점과 관련하여 주어진 각도 세트를 측정합니다 ( "상대 각도"라고 함). 상대 각도는 절대 값이 엄격하게 180도 미만입니다. 마지막으로, 이러한 상대 각도의 평균을 취하여 원하는 평균 각도를 얻습니다 (물론 원점 기준).


하나의 "정답"은 없습니다. 철저한 분석을 위해 KV Mardia와 PE Jupp, "Directional Statistics"(Wiley, 1999) 책을 읽는 것이 좋습니다.


영어로:

  1. 모든 각도가 180만큼 시프트 된 두 번째 데이터 세트를 만듭니다.
  2. 두 데이터 세트의 분산을 취하십시오.
  3. 가장 작은 분산으로 데이터 세트의 평균을 취하십시오.
  4. 이 평균이 쉬프트 된 세트에서 나온 경우 답을 180만큼 다시 쉬프트하십시오.

파이썬에서 :

#numpy NX1 각도 배열

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

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

여기에는 이동 평균을 사용하고 값을 정규화하는 데주의를 기울이는 완전히 산술적 인 솔루션이 있습니다. 모든 각도가 원의 한쪽에있는 경우 (서로 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;
}

글쎄, 나는 파티에 크게 늦었지만 실제로 어떤 정답도 찾을 수 없으므로 2 센트 가치를 더할 것이라고 생각했습니다. 결국 나는 간단하고 강력한 솔루션을 제공하기 위해 다음 Java 버전의 Mitsuta 메소드를 구현했습니다. 특히 표준 편차는 측정 분산을 제공하고 sd == 90 인 경우 입력 각도로 인해 평균이 모호하다는 것을 나타냅니다.

편집 : 실제로 나는 원래의 구현이 훨씬 단순화 될 수 있음을 깨달았습니다 . 실제로 다른 답변에서 진행되는 모든 대화와 삼각법을 고려할 때 걱정스럽게 간단합니다.

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

... 그리고 당신 (자바) 괴짜 모두를 위해, 위의 방법을 사용하여 한 줄에 평균 각도를 얻을 수 있습니다.

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

Alnitak은 올바른 솔루션을 제공합니다. Nick Fortescue의 솔루션은 기능적으로 동일합니다.

특별한 경우에

(sum (x_component) = 0.0 && sum (y_component) = 0.0) // 예를 들어 10도 및 190 도의 2 개 각도.

합계로 0.0도 사용

atan2 (0., 0.)이 정의되지 않았고 오류가 발생하므로 계산적으로이 경우를 테스트해야합니다.


평균 각도 phi_avg는 sum_i | phi_avg-phi_i | ^ 2가 최소가되는 특성을 가져야하며, 그 차이는 [-Pi, Pi) 여야합니다 (반대로 돌아가는 것이 더 짧을 수 있기 때문입니다). 이것은 모든 입력 값을 [0, 2Pi)로 정규화하고 실행 평균 phi_run을 유지하고 정규화 | phi_i-phi_run | [-Pi, Pi)에 (2Pi를 더하거나 빼서). 대부분의 제안은 위 않는 다른 것을 할 수 없는 , 즉 그 최소한의 속성을 가지고, 그들은 평균 뭔가를 ,하지만 각도.


(추정 이론 또는 통계적 추론의 관점을 공유하고 싶습니다)

민첩한 시도는 각도의 MMSE ^ 추정값을 얻는 것이지만 "평균화 된"방향을 찾는 것 중 하나입니다. MMAE ^ 추정치 또는 "평균화 된"방향으로 추정되는 다른 추정치를 찾을 수 있으며 이는 방향의 미터법 정량화 오류에 따라 다릅니다. 또는 더 일반적으로 추정 이론에서 비용 함수의 정의.

^ MMSE / MMAE는 최소 평균 제곱 / 절대 오류에 해당합니다.

ackb는 "평균 각도 phi_avg는 sum_i | phi_avg-phi_i | ^ 2가 최소가되는 특성을 가져야한다.

---- 평균 제곱의 의미로 오류를 수량화하지만 유일한 방법은 아니지만 가장 일반적인 방법 중 하나입니다. 여기서 대부분의 사람들이 선호하는 대답 (즉, 단위 벡터의 합과 결과의 각도를 얻는 것)은 실제로 합리적인 해결책 중 하나입니다. 벡터의 방향이 von Mises 분포로 모델링되는 경우 원하는 "평균화 된"방향의 역할을하는 것은 ML 추정기입니다. 이 배포판은 환상적이지 않으며 2D Guassian에서 주기적으로 샘플링 된 배포판입니다. 식을 참조하십시오. (2.179) Bishop의 저서 "패턴 인식 및 기계 학습". 다시 말하지만 이것이 "평균"방향을 나타내는 유일한 방법은 아니지만 이론적으로 정당화되고 간단한 구현이 가능한 것은 합리적입니다.

Nimble은 "이 벡터 기반 솔루션은 각도의 실제 평균으로 간주 될 수 없으며 단위 벡터 대응의 평균에 불과하다"고 말했다.

---- 이것은 사실이 아닙니다. "단위 벡터 대응 물"은 벡터 방향의 정보를 나타낸다. 각도는 벡터의 길이를 고려하지 않은 수량이며, 단위 벡터는 길이가 1이라는 추가 정보가있는 것입니다. "단위"벡터를 길이 2로 정의 할 수 있습니다. 실제로 중요하지 않습니다.


@David_Hanak의 답변을 통해 문제를 해결했습니다. 그가 말한대로 :

같은 반원을 유지하면서 다른 두 사람 사이를 "예"(예 : 355 및 5)하는 각도는 180이 아닌 0입니다. 이렇게하려면 두 각도의 차이가 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;

완벽하게 작동합니다.


파이썬 함수 :

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 : 모든 프로그래밍 언어에 대한 다음 링크에서 솔루션과 약간의 설명을 볼 수 있습니다

예를 들어 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;
}

산출:

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

starblue의 답은 평균 단위 벡터의 각도를 제공하지만 0에서 2 * pi 범위 (또는 0 °에서 2 °까지의 범위에서 둘 이상의 응답이있을 수 있음을 인정하는 경우 산술 평균의 개념을 각도로 확장 할 수 있습니다. 360 °). 예를 들어, 평균 0 ° 및 180 °는 90 ° 또는 270 ° 일 수 있습니다.

산술 평균은 입력 값에 대한 제곱 거리의 최소 합을 갖는 단일 값의 속성을 갖습니다. 두 단위 벡터 사이의 단위 원을 따른 거리는 내적의 역 코사인으로 쉽게 계산할 수 있습니다. 벡터와 각 입력 단위 벡터의 내적의 제곱 역 코사인의 합을 최소화하여 단위 벡터를 선택하면 평균이 같습니다. 예외적 인 경우 최소 2 개 이상이있을 수 있습니다.

단위 구를 따른 거리는 단위 원을 따른 거리 (두 단위 벡터의 내적의 역 코사인)와 정확히 같은 방식으로 계산 될 수 있기 때문에이 개념은 여러 차원으로 확장 될 수 있습니다.

원의 경우이 평균을 여러 가지 방법으로 풀 수 있지만 다음 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으로 바꾸십시오. 더 많은 차원을 사용하려면 평균을 풀기 위해 반복적 인 방법을 사용해야 할 것입니다.


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 음이 아닌 모든 각도를 추가하고 평균을 계산 한 다음 개수를 계산합니다. b. 모든 음의 각도를 추가하고 평균을 계산하고 개수를 계산합니다. 3. pos_average에서 neg_average를 뺀 차이 차이가 180보다 크면 차이를 360에서 뺀 차이로 변경하십시오. 그렇지 않으면 차이의 부호 만 변경하십시오. 차이는 항상 음이 아닙니다. Average_Angle은 pos_average와 차이에 "weight"를 더한 음수, 음수를 음수와 양수의 합으로 나눈 값과 같습니다.


다음은 평균 각도에 대한 Java 코드입니다. 합리적으로 강력하다고 생각합니다.

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

위의 일부 각도에 대해 "올바른"답변을 제공하는 @Starblue와 다른 방법이 있습니다. 예를 들면 다음과 같습니다.

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

연속 각도의 차이에 대한 합을 사용합니다. 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

참고 URL : https://stackoverflow.com/questions/491738/how-do-you-calculate-the-average-of-a-set-of-circular-data

반응형