Pearson correlation between $ X $ and $ X^n $

2020-03-25

2020-11-27

matplotlib, person correlation, python, statistics, visualization

wikipedia: Pearson Correlation Coefficient
선형대수, 통계학, 기하적 직관

Pearson Correlation Coefficient는 두 데이터가 얼마나 연관성을 가지고 있는지 보여줍니다.
수치로는 -1에서 1 사이의 값으로 표현되고, 수식으로는 다음과 같이 표현됩니다.
$$ \rho_{X, Y} = \frac{Cov (X, Y)}{\sigma_X \sigma_Y} \text{, where } Cov \text{ is the covariance.}$$
이를 $X$와 $Y$의 평균($ \mu $)과 기댓값($ \mathbb{E}$)으로 표현하면 다음과 같습니다.
$$ \begin{aligned} \rho_{X, Y} &= \frac{\mathbb{E}[(X - \mu_X)(Y - \mu_Y)]}{\sigma_X \sigma_Y} \ &= \frac{\mathbb{E}[XY]-\mathbb{E}[X]\mathbb{E}[Y]}{\sqrt{\mathbb{E}[X^2]-[\mathbb{E}[X]]^2}\sqrt{\mathbb{E}[Y^2]-[\mathbb{E}[Y]]^2}} \end{aligned}$$

기하학적 직관을 바탕으로 $Y$의 $X$로의 projection이라는 통찰을 보이신 분도 있습니다만
간단하게 “$ Y=X $와 같은 직선에서 벗어나면 멀어진다” 정도는 알고 있기에,
$X$와 $X^n$의 상관계수가 궁금해졌습니다.

1. Pearson Correlation Coefficient

$ X^n $은 X의 범위에 따라 X에서 벗어나는 정도가 다릅니다.
X가 0~~1일 때, 그리고 1~~250일 때 두 가지 경우를 나누어 테스트해 보았습니다.
- X가 0~1일 때 데이터: test1.csv
  1
  2
  df1 = pd.read_csv('test1.csv')
  df1.head()
- X가 1~250일 때 데이터: test250.csv
  1
  2
  df250 = pd.read_csv('test250.csv')
  df250.head()

X1~X5는 각기 $ X^n$ $(n=1,2,3,4,5)$이고,
Y2~Y5는 각기 $\displaystyle \sum_{k=1}^n X^k$ $(n=2,3,4,5)$ 입니다.

두 데이터 셋에 $ X $와 더불어 $ X^n$ $(n=1,2,3,4)$를 담았습니다.
- 그래프를 그려보면 다음과 같습니다.

pandas의 DataFrame.corr()함수를 이용하여 PCC를 간단하게 구할 수 있습니다.
1
df1[['X', 'X2', 'X3', 'X4', 'X5']].corr(method='pearson')
- 실행결과 : $X$와 가장 먼 $X^5$도 0.82나(!) 됩니다.

이번에는 훨씬 차이가 나는 1~250까지를 봅시다.
1
2
df250 = pd.read_csv('test250.csv')
df250[['X', 'X2', 'X3', 'X4', 'X5']].corr(method='pearson')
- 실행결과 : 0~1까지와 차이가 거의 없습니다.
  소숫점 5째 자리는 가야 차이가 드러납니다.

일단 호기심이 해결되었습니다.

2. Feature Importance

matplotlib.axes.Axes.violon

그렇다면, $ Y = X + X^2 + \ldots + X^n $의 회귀식을 만들었을 때, 각각의 기여도는 얼마일까요.
- 회귀식은 무조건 1 : 1 : … : 1이 나올 것입니다.
- Random Forest의 Feature Importance를 통해서 판단해보겠습니다.

Random Forest는 데이터 분할에 따라 결과가 달라집니다.

다음과 같은 코드를 작성하여 100회 반복한 결과를 수집했습니다.

from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score

def calc_featimps(df):
    
    featimps = {}
    for i in range(2,6):
        y = df1[f'Y{i}']
        x = df1.iloc[:,:i]
        print(f'X features={x.columns.tolist()}, Y={y.name}')

        featimp = pd.DataFrame()

        for j in range(100):
            X_train, X_test, y_train, y_test = train_test_split(x, y)
            rf = RandomForestRegressor()
            rf.fit(X_train, y_train)
            y_pred = rf.predict(X_test)
            r2 = r2_score(y_test, y_pred)

            featimp_ = pd.Series(rf.feature_importances_, index=X_train.columns, name=f'{j}')
            featimp = pd.concat([featimp, featimp_], axis=1)

        featimps[f'{i}'] = featimp  
    
    return featimps

matplotlib의 violon plot을 이용해서 분포를 확인해 보겠습니다.

Feature별 Feature Importance의 평균값과 표준편차를 함께 표시하겠습니다.

plt.rcParams['mathtext.fontset'] = 'cm'
plt.rcParams['mathtext.default'] = 'it'  
font_title = {'fontsize': 15,
              'fontweight' : plt.rcParams['axes.titleweight'],
              'verticalalignment': 'center',
              'horizontalalignment': 'center'}
pad_title = 12

font_label = {'verticalalignment': 'bottom',
             'horizontalalignment': 'center'}
pad_label = 20
  
def plot_featimps():
    fig, axes = plt.subplots(ncols=4, nrows=1, figsize=(10,4))
    
    xlabels = [r'$ X $', r'$ X^2 $', r'$ X^3 $', r'$ X^4 $', r'$ X^5 $']
    for i, ax in enumerate(axes):
        labels = featimps[f'{i+2}'].index.tolist()
        print(labels)

        positions=np.arange(i+2).tolist()
        ax.violinplot([featimps[f'{i+2}'].loc[idx] for idx in labels], 
                      showextrema=False, positions=positions)
        ax.errorbar(labels, featimps[f'{i+2}'].mean(axis=1), yerr=featimps[f'{i+2}'].std(axis=1), 
                    fmt='o', mec=f'C{i}', mew=3, mfc='w', ms=15,
                    lw=2, c=f'C{i}')

        ax.set_xlim((-0.5, i+1.5))
        ax.set_xticks(np.arange(i+2))
        ax.set_xticklabels(xlabels[:i+2], fontdict=font_label)
        ax.tick_params(axis='x', which='major', pad=pad_label)
  
    axes[0].set_ylim((0.3, 0.7))
    axes[1].set_ylim((0.2, 0.5))
    axes[2].set_ylim((0.1, 0.4))
    axes[3].set_ylim((0.1, 0.3))

    axes[0].set_title(r'$Y=X+X^2$', fontdict=font_title, pad=pad_title)
    axes[1].set_title(r'$Y=X+X^2+X^3$', fontdict=font_title, pad=pad_title)
    axes[2].set_title(r'$Y=X+X^2 $      '+'\n' + r'$ +X^3+X^4$', fontdict=font_title, pad=pad_title)
    axes[3].set_title(r'$Y=X+X^2 $             '+'\n' + r'$ +X^3+X^4+X^5$', fontdict=font_title, pad=pad_title)

    plt.tight_layout()

2.1. X = 0 ~ 1

먼저 X가 0과 1 사이일 때를 살펴보겠습니다.
1
2
featimps1 = calc_featimps(df1)
plot_featimps()
- 실행결과: 모든 경우에서 Feature Importance가 균등하게 분배됩니다.
  Random Forest 수행에 따른 편차가 violin plot으로 드러나 있기는 하지만요.

2.2. X = 1 ~ 250

이제 X가 1과 250 사이일 때를 살펴보겠습니다.
1
2
featimps1 = calc_featimps(df250)
plot_featimps()
- 실행결과: 앞서와 같은 결과입니다.

2.3. 정리

Random Forest를 수행한 경우에도, 정상적으로 평균 Feature Importance는 $\frac{1}{n}$입니다.
수정되기 전의 글에서는 특정 차수가 높게 나왔는데, 데이터 오류 때문으로 확인되었습니다.

도움이 되셨나요? 카페인을 투입하시면 다음 포스팅으로 변환됩니다

PythonData Science

1. Pearson Correlation Coefficient

2. Feature Importance

2.1. X = 0 ~ 1

2.2. X = 1 ~ 250

2.3. 정리