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

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

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

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

답변

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