시계열 데이터의 상관성 구하기 (time-series correlation)

서로 다른 시계열 데이터의 상관성을 어떻게 알 수 있을까?

그리고 두 데이터의 길이가 다르다면??

 

공통적인 주의사항으로, 상관관계를 인과관계로 해석해서는 안된다는 것임을 통계학 수업에서 많이 들어봤을 것이다.

1. Pearson Correlation Coefficient (PCC, Pearson's r)

통계 시간에 배우는 그 피어슨-상관계수 맞다. 

 

\[ \rho_{X, \ Y} = \cfrac{\text{cov}(X, \ Y)}{\sigma_X \sigma_Y} = \cfrac{\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} }\]

$\text{cov}(X, \ Y)$는 공분산이고, $\sigma$는 표준편차, $\mathbb{E}$는 기댓값이다.

 

표본의 경우 아래처럼 구하고 sample correlation coefficient (or sample PCC)

\[ r_{xy} = \cfrac{\sum x_i y_i - n \bar{x} \bar{y}}{(n-1) s_x s_y} \]

이때 $n$은 표본의 크기, $\bar{x}$는 표본평균, $s_x$는 표본표준편차(표본분산의 루트)이다.

 

가장 간단하면서 단점(과 함정)이 많은 상관계수이다.

우선 두 데이터간의 선형적 관계만 파악할 수 있다는 것이 단점이다. 다른 관계는 찾지 못한다.

 

그리고 2가지 함정이 있는데 첫번째로는 이상치(outlier)가 상관계수의 질을 크게 악화(degraded correlation values)시킨다는 것이고, 두번째로는 PCC는 특정 조건을 만족하는 이변량 분포에서만 정의된다. (모집단의 공분산이 정의되고 그 값이 0이 아니어야 한다. 코시 분포의 경우 분산이 정의되지 않는다.)

 

그리고 애초에, 데이터의 길이가 다르면 사용할 수 없다.

data1과 data2는 상관성이 높아보이지만, 길이가 달라 계산할 수 없다.

2. Cross Correlation (교차상관, 상호상관)

신호처리에서 사용되는 방법으로, 두 series가 선후 관계에 있는지 유사도를 측정하는 방법이다. 직관적으로 설명하면, 두 수열 $x_t, y_t$에 대하여 $x$가 $y$보다 선행하는지, 후행하는지 알 수 있다는 것이다.

 

연속적인 분포에서 아래와 같이 정의된다.

\[ f \star g(\tau) = \int_{-\infty}^{\infty} f(t)g(t + \tau) dt \]

 

이산적인 분포는 아래와 같이 정의된다.

\[ \cfrac{1}{n} \sum_{x, y} f(x, y)t(x, y) \]

 

(convolution은 * 를, cross-correlation은 $\star$를 연산기호로 사용하고, convolution은 $f(t)$와 $g(-t)$의 cross-correlation과 연산이 같다고 한다.)

 

여담으로, 이미지처리분야에서 어떤 subimage $f(x, y)$에서 template image $t(x, y)$를 찾는 template-matching task에서도 cross-correlation이 사용되는데, 평균 또는 표준편차로 정규화할때마다 수식이 달라진다. (수식 형태를 보면 scaling issue만 있을 뿐 본질은 같다는 것을 알 수 있다)

\[ ZNCC = \cfrac{1}{n} \sum_{x, y}\cfrac{1}{\sigma_f \sigma_t}(f(x, y) - \bar{f})(t(x, y) - \bar{t}) \]

\[ NCC = \cfrac{1}{n} \sum_{x, y} \cfrac{1}{\sigma_f \sigma_t}f(x, y)t(x,y) \]

 

다시 시계열로 돌아와서, Time delay analysis에서 유용하게 사용된다. (앞서 언급한 선행/동행/후행 판단) 두 신호(시계열로 보아도 무방)의 cross correlation을 계산하여 cross-correlation function의 최댓값(또는 음의 상관의 경우는 최솟값)은 두 신호가 가장 잘 정렬된 시점을 나타낸다. 따라서

\[ \tau_{\text{delay}} = \text{argmax}_{t \in \mathbb{R}} (f \star g) (t) \]

 

그러나 cross-correlation도 해석할 때 주의할 점이 있다.

우선, 비선형계(non-linear system)에서 적용할 때 시스템의 입력/출력의 cross-correlation은 가려져 나타나지 않을 수 있다. 통계적 의미인 상관성이 수치로 드러나지 않을 수 있다는 것.

두번째는로 주의할 점은 자기상관관계(autocorrelation)가 없다는 가정이 필요하다는 것이다.

참고로 자기상관(autocorrelation)은 서로 다른 신호가 아니라 같은 신호가 correlate하다는 것이다. 자기상관의 정의는 다음과 같다. (신호처리에 문외한인 나로써는 재귀적인 구조처럼 보인다)

\[ f \star f( \tau ) = \int_{-\infty}^{\infty} f(t) f(t + \tau) dt \]

 

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(0, 6.28, num=50)
y1 = np.sin(x)
y2 = np.sin(x - np.pi / 2)

# plot the data
fig, ax = plt.subplots()
ax.plot(x, y1, label='data1', marker='o')
ax.plot(x, y2, label='data2', marker='o')
ax.legend()

# get Cross-Correlation(y1, y2)
corr = np.correlate(y1, y2, mode='full')
max_idx = np.argmax(corr)
min_idx = np.argmin(corr)
plt.plot(range(len(corr)), corr)
plt.vlines(x=max_idx, color='r', ymin=0, ymax=corr[max_idx], label='max-cross-correlation')
plt.vlines(x=min_idx, color='b', ymin=corr[min_idx], ymax=0.0, label='min-cross-correlation')
plt.hlines(y=0, color='black', xmin=0, xmax=len(corr), linestyle='dashed')
plt.legend()

두 사인파의 cross-correlation

반응형

3. Dynamic Time Warping (DTW)

DTW는 길이가 다른 시계열 데이터에 대하여 상관성을 도출한다. 그리고 비-유클리드 거리를 이용한 것이기 때문에 피어슨-상관계수와 달리 선형적 관계를 지니지 않아도 두 데이터의 상관성을 얻어낼 수 있다.

DTW는 비유클리드 거리를 이용하여 두 데이터의 거리를 계산한다.

 

DWT는 컴퓨터 알고리즘 중에서 Dynamic Programming을 이용한 방법이다. DTW의 distance를 구하는 방법은 마치 학부 수업 중에서 sequence alignment와 아이디어가 거의 같다.

 

두 배열 $s[1 \cdots n]$과 $t[1 \cdots m]$의 DTW 거리를 구하는 공식은 다음과 같다. 

(1) 초기에 $DTW(i, \ j) = \infty$로 초기화하고 $DTW(0, 0) = 0$으로 초기화한다. 

(2) $i = 1, \cdots n, \ j = 1, \cdots m$ 동안 아래를 반복한다

\[ DTW(i,\ j) = \lvert s_i - t_j \rvert + \min\left(DTW(i- 1, \ j), \ DTW(i, \ j-1), \ DTW(i-1, \ j-1) \right) \]

 

 

이렇게 생긴 DTW 행렬을 costMatrix로 부른다. 보다시피 시간복잡도는 $O(nm)$이다.

 

시계열 데이터를 어떤 미분가능한 함수의 이산형으로 간주하여 빠르게 DTW를 구하는 방법도 존재한다.

오픈소스 중 FastDTW는 다양한 방법과 재귀적인 방법으로 $O(N)$ 만에 DTW를 구한다고 (주장)한다. 

 

# generate data
start = 0
end = 2 * np.pi
step = 0.1
k = 2
 
x1 = np.arange(start, end, k * step)
x2 = np.arange(start, end / k, step)
 
noise = np.random.uniform(start, end, len(x2)) / 10
 
y1 = np.sin(x1) + np.sin(2 * x1) + noise
y2 = np.sin(k * x2) + np.sin(k * 2 * x2)

fig, ax = plt.subplots()
ax.plot(x1, y1, label='data1', marker='o')
ax.plot(x2, y2, label='data2', marker='o')
ax.legend()

# apply DTW algorithm
import dtw

alignment = dtw.dtw(y1, y2, keep_internals=True)

# costMatrix
plt.imshow(alignment.costMatrix, interpolation='nearest', cmap='gray_r')
plt.colorbar()

# print DTW distance
print(alignment.distance, alignment.normalizedDistance)

# plot how to match the two data
alignment.plot(type="twoway")
alignment.plot(type="threeway")

(좌) 랜덤 생성한 데이터. (우) DTW의 cost matrix를 시각화. 밝을 수록 cost가 작다.

(좌) two way로 어떻게 최소 거리를 매칭했는지 시각화. (우) three way로 거리 시각화

 

출처 및 참고

https://en.wikipedia.org/wiki/Pearson_correlation_coefficient

Pearson correlation coefficient - Wikipedia

https://en.wikipedia.org/wiki/Cross-correlation

Cross-correlation - Wikipedia

https://en.wikipedia.org/wiki/Dynamic_time_warping

Dynamic time warping - Wikipedia

https://github.com/rmaestre/FastDTW

GitHub - rmaestre/FastDTW: FastDTW- Dynamic Time Warping (DTW) with a linear time and memory complexity