outlier 빼고 상관계수 구하기 : Mahalanobis 거리

   데이터를 분석할 때 거의 항상 문제가 되는 것이 '튀는 점', 즉 outlier 이다. 보통 1차원이고 정규분포를 따르는 데이터일 때는 아주 간단하게는 평균 ± d*표준편차 밖에 있는 점들을 outlier 라고, d 는 2나 1.96 등의 값을 사용한다. 그런데 만약 2차원 이상의 데이터라면 계산값이 조금 달라진다. 데이터 분포가 다음과 같을 때를 생각해 보자.

  데이터가 위와 같은 분포를 보일 때, 중심에서의 거리는 파란색 점보다는 녹색 점이 더 먼 것으로 생각할 수 있다. 왜냐 하면, 분산이 파란색 점이 있는 축이 더 작기 때문이다. 즉, 적게 퍼져 있는 분포에서 1 먼 것과, 많이 퍼져 있는 분포에서 1 먼 것은, 상대적으로 비교해 보면 많이 퍼져 있는 곳에서 1 먼 것이 좀 더 가깝다고 할 수 있기 때문이다. 이와 같은 개념을 Eucledian distance 에 적용시킨 것이 Mahalanobis distance 이다. 즉, 각 축(basis)에서의 거리를 구할 때, 원래의 값을 그 축에 해당하는 표준편차로 나누어 주는 것이다. Eucledian distance는 모든 축의 표준편차가 1 인 Mahalanobis 거리가 되는 것이다. 구체적인 식은, 위키에 있는 식을 가져 오면 다음과 같다.

이 때 S는 covariance matrix 이고, u 는 각 축의 평균이다. 다시 좀 더 쉽게 2차원 (x,y) 에서 생각을 해보면, u = (x 평균, y 평균), x = (x, y)이고, Mahalanobis distance D_M(x)는 다음과 같다.

보 통 이 단계까지 한 이후부터가 정말로 어려워지는 부분인데, Mahalanobis distance (MD) 가 일정한 분포를 이룰 때, 어느 값부터가 outlier 라 할 것인가? MD 값이 1,1,1,1,1,1,1,3,6 이런 식이면 MD 값의 평균+2*표준편차 밖에 있는 것으로 하면 6 이 쉽게 잡힐 것이다. 문제는 과연 [평균 - 2*표준편차, 평균 + 2*표준편차] 이 구간 이외의 값을 outlier 라 하는 것이 타당하느냐 하는 것인데, 때에 따라 그렇지 않은 경우가 많이 있다. 이 문서[1]에서 간략히 그러한 내용을 볼 수 있다. 문서 [1] 은 원래 [2]에 대한 간략한 예이다. [2]에 있는 원래 내용은 실제 관측된 MD 분포가 chi-square 분포를 따르고, 원래 데이터가 정규분포를 따른다면 관측된 MD 분포가 정규분포로 수렴할 것이므로, 분포의 가장 끝단만을 확인해서 데이터에 맞게끔 outlier 가 되는 p-value 를 정한다는 내용이다.

  나는 일단 [1]번과 [2]번에 나온 내용을 구현하기에는 배보다 배꼽이 더 커질 것 같아서, 우선 MD 거리를 큰 순서로 나열한 후 제일 큰 5%의 데이터를 outlier 라 가정했다 (percentile로 trimming. 앵겔지수도 이렇게 percentile regression 에서 나왔다고 한다). 결과는 다음과 같다.

빨 간색과 파란색 선은 각각의 네모에 대한 추세선이다. 점선은 모든 데이터를 고려했을 때의 추세선이고, 실선은 위에서 말한대로 튀는 데이터들(분홍색으로 표시한)을 제외한 후 계산한 추세선이다. 데이터들이 몰려 있는 곳에서 멀리 떨어져 있는 데이터들이 제거되는 것을 볼 수 있다. 그리고 그것 때문에 빨간색 데이터의 경우 상관계수 값이 양수에서 음수로 변한 것을 알 수 있다, 물론 이 경우 R² 값이 작아서 상관계수가 별로 의미가 없긴 하지만.


   구현의 경우, covariance matrix의 역행렬을 구해야 하는데, 일단 급한데로 중학교(고등학교?) 때 배운 역행렬 공식을 그냥 때려 넣어서 구했다. n-by-n matrix로까지 확장시키기 위해서는 많은 노력이 필요해서...

이에 대한 코드는 다음과 같다. 필요한 함수의 구현까지 모두 첨부하기 때문에 좀 길어졌는데, Mahalanobis 거리를 이용한 예제를 보기 위해서는, 아래 코드 중 제일 마지막에 있는 함수인 get_robust_correlation 부분을 보면 된다. 이 함수는, 데이터를 Mahalanobis 거리로 정렬한 후 거리가 제일 멀리 떨어진 특정 %의 데이터를 날리고 Pearson's correlation 을 구한다.


[code cpp]
#include <numeric>

#include <math.h>
bool get_median(std::vector<double> src, double* dest)
{
    if(src.empty() == true) return false;
    std::sort(src.begin(), src.end());
    int half_pos = (int)(((int)(src.size()) * 0.5 + 0.5));
    *dest = src[half_pos];
    return true;
}
bool get_mean_variance(const std::vector<double>& x, double *mean, double *var)
{
    if(x.empty() == true) return false;
    int n = (int)(x.size());
    std::vector<double>::const_iterator pos = x.begin();
    double sum = 0;
    double square_sum = 0;
    for(; pos != x.end(); pos++){
        sum += *pos;
        square_sum += (*pos * *pos);
    }

    double m = sum/n;
    if(mean != NULL) *mean = m;
    *var = square_sum/n - m*m;

    return true;
}
bool get_mean_variance(const std::set<double>& x, double *mean, double *var)
{
    if(x.empty() == true) return false;
    int n = (int)(x.size());
    std::set<double>::const_iterator pos = x.begin();
    double sum = 0;
    double square_sum = 0;
    for(; pos != x.end(); pos++){
        sum += *pos;
        square_sum += (*pos * *pos);
    }

    double m = sum/n;
    if(mean != NULL) *mean = m;
    if(var != NULL) *var = square_sum/n - m * m;

    return true;
}
bool get_covariance(const std::vector<double>& x, const std::vector<double>& y, double *cov)
{
    if(x.size() != y.size()) return false;
    std::vector<double>::const_iterator xpos = x.begin();
    std::vector<double>::const_iterator ypos = y.begin();
    double cross_sum = 0;
    double x_sum = 0;
    double y_sum = 0;
    for(; xpos != x.end(); xpos++,ypos++){
        x_sum += (*xpos);
        y_sum += (*ypos);
        cross_sum += (*xpos * *ypos);
    }
    int n = (int)(x.size());
    *cov = cross_sum/n - x_sum*y_sum/n/n;


    return true;
}

// y_i = _alpha * x_i + _beta + epsilon
bool get_correlation(std::vector<double> *x, std::vector<double>* y, double *corr, double *Rsquare, double *_alpha, double* _beta)
{
    if(x->size() != y->size()) return false;

    int n = (int)(x->size());
    if(n == 0) return false;

    std::vector<double>::const_iterator xpos = x->begin();
    std::vector<double>::const_iterator ypos = y->begin();
   
    double x_sum = 0;
    double y_sum = 0;
    double x_square_sum = 0;
    double y_square_sum = 0;
    double xy_sum = 0;

    for(; xpos != x->end(); xpos++,ypos++){
        x_sum += *xpos;
        y_sum += *ypos;
        x_square_sum += (*xpos * *xpos);
        y_square_sum += (*ypos * *ypos);
        xy_sum += (*xpos * *ypos);
    }

    if(corr != NULL){
        *corr = (n*xy_sum-x_sum*y_sum)/(sqrt(n*x_square_sum-x_sum*x_sum)*sqrt(n*y_square_sum-y_sum*y_sum));
    }
    double alpha = (x_sum*y_sum-n*xy_sum)/(x_sum*x_sum-n*x_square_sum);
    double beta = (x_sum*xy_sum - y_sum*x_square_sum)/(x_sum*x_sum-n*x_square_sum);
    if(_alpha != NULL){
        *_alpha = alpha;
    }
    if(_beta != NULL){
        *_beta = beta;
    }
    if(Rsquare == NULL){
        return true;
    }

    double y_mean = y_sum/n;
    double SStot = 0;
    double SSerr = 0;
    xpos = x->begin();
    ypos = y->begin();
    for(; xpos != x->end(); xpos++,ypos++){
        SStot += ((*ypos - y_mean)*(*ypos - y_mean));
        SSerr += ((*ypos - alpha*(*xpos) - beta)*(*ypos - alpha*(*xpos) - beta));
    }

    if(Rsquare != NULL){
        *Rsquare = 1 - SSerr/SStot;
    }


    return true;
}
bool sort_by_maha(const std::pair<int,double>& op1, const std::pair<int,double>& op2)
{
    return op1.second > op2.second;
}
bool get_robust_correlation(double trim_percent, const std::vector<double>& x, const std::vector<double>& y, double *corr, double *Rsquare, double *_alpha, double *_beta)
{
    //if(x == NULL || y == NULL) return false;
    if(x.empty() == true || y.empty() == true) return false;
    if(x.size() != y.size() ) return false;

    int n = (int)(x.size());
    double rho = 0;
    //get_correlation(x, y, &rho, NULL, NULL, NULL);

   

    //double x_sum = std::accumulate(x->begin(), x->end(), 0.0);
    //double y_sum = std::accumulate(y->begin(), y->end(), 0.0);
    double x_mean = 0;
    double y_mean = 0;
    double x_var = 0;
    double y_var = 0;

    get_mean_variance(x, &x_mean, &x_var);
    get_mean_variance(y, &y_mean, &y_var);

    double cov_xy = 0;
    get_covariance(x,y, &cov_xy);

    double denominator = 1.0/(x_var*y_var - cov_xy*cov_xy);
    double a = y_var/denominator;
    double b = -1.0 * cov_xy/denominator;
    double c = b;
    double d = x_var / denominator;

    // Now, S^-1 = | a     b|
    //             | c    d|
    // where S is covariance matrix of x and y.

    double maha_dist_sum = 0;
    double maha_dist_square_sum = 0;
    std::vector<std::pair<int,double> > idx2maha;
    std::vector<double>::const_iterator xpos = x.begin();
    std::vector<double>::const_iterator ypos = y.begin();
    int idx = 0;
    for(; xpos != x.end(); xpos++,ypos++,idx++ ){
        //double maha_dist = (*xpos - x_mean)*(*xpos - x_mean) + 2*rho*(*xpos - x_mean)*(*ypos - y_mean) + (*ypos - y_mean)*(*ypos - y_mean);
        double x_diff = *xpos - x_mean;
        double y_diff = *ypos - y_mean;
        double maha_dist = a*x_diff*x_diff + (b+c)*x_diff*y_diff + d*y_diff*y_diff;
        maha_dist_sum += maha_dist;
        maha_dist_square_sum += (maha_dist*maha_dist);
        idx2maha.push_back(std::make_pair<int,double>(idx,maha_dist));
    }

    std::sort(idx2maha.begin(), idx2maha.end(), sort_by_maha);

    int min_threshold_index = (int)(trim_percent*n+0.5);
    std::vector<double> new_x; std::vector<double> new_y;
    int i = 0;
    for(i = min_threshold_index; i<n; i++){
        new_x.push_back(x.operator[](idx2maha[i].first));
        new_y.push_back(y.operator[](idx2maha[i].first));
    }

    get_correlation(&new_x, &new_y, corr, Rsquare, _alpha, _beta);

    return true;
}

[/code]


[1] A multivariate outlier detection method (pdf임) by P. Filzmoser
[2] A robust and efficient adaptive reweighted estimator of multivariate location and scatter by Daniel Gervini, Journal of Multivariate Analysis, 2003, vol.84, pp-116-144