[譯] 時間序列數據間量化同步的四種方法

• 原文地址：Four ways to quantify synchrony between time series data
• 原文做者：Jin Hyun Cheong
• 譯文出自：掘金翻譯計劃
• 本文永久連接：github.com/xitu/gold-m…
• 譯者：EmilyQiRabbit
• 校對者：zhmhhu

時間序列數據間量化同步的四種方法

用於計算同步指標的示例代碼和數據包括：皮爾遜相關，時間滯後互相關，動態時間扭曲和瞬時相位同步。

在心理學當中，人與人之間的同步性是能提供社會動態和潛在社交產出的重要信息。它能夠體現於衆多領域，包括肢體動做（Ramseyer & Tschacher, 2011）、面部表情（Riehle, Kempkensteffen, & Lincoln, 2017）、瞳孔的擴張（Kang & Wheatley, 2015）以及神經信號（Stephens, Silbert, & Hasson, 2010）。不管如何，同步性能夠提供多種的意義，同時，量化兩個信號的同步性也有不少方法。

在本篇文章中，我調研了一些最經常使用的同步指標的利弊，並衡量了包括：皮爾遜相關，時間滯後互相關（TLCC）以及加窗的 TLCC，動態時間扭曲和瞬時相位同步這幾種技術。爲了更好的說明，同步指標會使用樣本數據計算，樣本數據是一段三分鐘且包含兩個參與者對話的視頻，咱們將從中提取面部微笑表情（下圖是一個截屏）。爲了讓你能更好的跟上文章的內容，能夠免費下載從樣本中提取的面部數據以及包括了全部示例代碼的Jupyter 筆記。

目錄

1. 皮爾遜相關
2. 時間滯後互相關（TLCC）以及加窗的 TLCC
3. 動態時間扭曲（DTW）
4. 瞬時相位同步

---

1. 皮爾遜相關 —— 最簡單也是最好的方法

皮爾遜相關能夠衡量兩個連續信號如何隨時間共同變化，而且能夠以數字 -1（負相關）、0（不相關）和 1（徹底相關）表示出它們之間的線性關係。它很直觀，容易理解，也很好解釋。可是當使用皮爾遜相關的時候，有兩件事情須要注意，它們分別是：第一，異常數據可能會干擾相關評估的結果；第二，它假設數據都是同方差的，這樣的話，數據方差在整個數據範圍內都是同質的。一般狀況下，相關性是全局同步性的快照測量法。因此，它不能提供關於兩個信號間方向性的信息，例如，哪一個信號是引導信號，哪一個信號是跟隨信號。

不少包都應用了皮爾遜相關，包括 Numpy、Scipy 和 Pandas。若是你的數據中包含了空值或者缺失值，Pandas 中的相關性方法將會在計算前把這些行丟棄，而若是你想使用 Numpy 或者 Scipy 對於皮爾遜相關的應用，你則必須手動清除掉這些數據。

以下的代碼加載的就是樣本數據（它和代碼位於同一個文件夾下），並使用 Pandas 和 Scipy 計算皮爾遜相關，而後繪製出了中值濾波的數據。

import pandas as pd
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
import scipy.stats as stats

df = pd.read_csv('synchrony_sample.csv')
overall_pearson_r = df.corr().iloc[0,1]
print(f"Pandas computed Pearson r: {overall_pearson_r}")
# 輸出：使用 Pandas 計算皮爾遜相關結果的 r 值：0.2058774513561943

r, p = stats.pearsonr(df.dropna()['S1_Joy'], df.dropna()['S2_Joy'])
print(f"Scipy computed Pearson r: {r} and p-value: {p}")
# 輸出：使用 Scipy 計算皮爾遜相關結果的 r 值：0.20587745135619354，以及 p-value：3.7902989479463397e-51

# 計算滑動窗口同步性
f,ax=plt.subplots(figsize=(7,3))
df.rolling(window=30,center=True).median().plot(ax=ax)
ax.set(xlabel='Time',ylabel='Pearson r')
ax.set(title=f"Overall Pearson r = {np.round(overall_pearson_r,2)}");
複製代碼

再次重申，全部的皮爾遜 r 值都是用來衡量全局同步的，它將兩個信號的關係精簡到了一個值當中。儘管如此，使用皮爾遜相關也有辦法觀察每一刻的狀態，即局部同步性。計算的方法之一就是測量信號局部的皮爾遜相關，而後在全部滑動窗口重複該過程，直到全部的信號都被窗口覆蓋過。因爲能夠根據你想要重複的次數任意定義窗口的寬度，這個結果會因人而異。在下面的代碼中，咱們使用 120 幀做爲窗口寬度（4 秒左右），而後在下圖展現出咱們繪製的每一刻的同步結果。

# 設置窗口寬度，以計算滑動窗口同步性
r_window_size = 120
# 插入缺失值
df_interpolated = df.interpolate()
# 計算滑動窗口同步性
rolling_r = df_interpolated['S1_Joy'].rolling(window=r_window_size, center=True).corr(df_interpolated['S2_Joy'])
f,ax=plt.subplots(2,1,figsize=(14,6),sharex=True)
df.rolling(window=30,center=True).median().plot(ax=ax[0])
ax[0].set(xlabel='Frame',ylabel='Smiling Evidence')
rolling_r.plot(ax=ax[1])
ax[1].set(xlabel='Frame',ylabel='Pearson r')
plt.suptitle("Smiling data and rolling window correlation")
複製代碼

總的來講，皮爾遜相關是很好的入門學習教程，它提供了一個計算全局和局部同步性的很簡單的方法。可是，它不能提供信號動態信息，例如哪一個信號先出現，而這個能夠用互相關來衡量。

2. 時間滯後互相關 —— 評估信號動態性

時間滯後互相關（TLCC）能夠定義兩個信號之間的方向性，例如引導-追隨關係，在這種關係中，引導信號會初始化一個響應，追隨信號則重複它。還有一些其餘方法能夠探查這類關係，包括格蘭傑因果關係，它經常使用於經濟學，可是要注意這些仍然不必定能反映真正的因果關係。可是，經過查看互相關，咱們仍是能夠提取出哪一個信號首先出現的信息。

如上圖所示，TLCC 是經過逐步移動一個時間序列向量（紅色線）並反覆計算兩個信號間的相關性而測量獲得的。若是相關性的峯值位於中心（offset=0），那就意味着兩個時間序列在此時相關性最高。可是，若是一個信號在引導另外一個信號，相關性的峯值就可能位於不一樣的座標值上。下面這段代碼應用了一個使用了 pandas 提供功能的互相關函數。同時它也能夠將數據打包，這樣相關性邊界值也能經過添加信號另外一邊的數據而計算出來。

def crosscorr(datax, datay, lag=0, wrap=False):
    """ Lag-N cross correlation. Shifted data filled with NaNs Parameters ---------- lag : int, default 0 datax, datay : pandas.Series objects of equal length Returns ---------- crosscorr : float """
    if wrap:
        shiftedy = datay.shift(lag)
        shiftedy.iloc[:lag] = datay.iloc[-lag:].values
        return datax.corr(shiftedy)
    else: 
        return datax.corr(datay.shift(lag))

d1 = df['S1_Joy']
d2 = df['S2_Joy']
seconds = 5
fps = 30
rs = [crosscorr(d1,d2, lag) for lag in range(-int(seconds*fps-1),int(seconds*fps))]
offset = np.ceil(len(rs)/2)-np.argmax(rs)
f,ax=plt.subplots(figsize=(14,3))
ax.plot(rs)
ax.axvline(np.ceil(len(rs)/2),color='k',linestyle='--',label='Center')
ax.axvline(np.argmax(rs),color='r',linestyle='--',label='Peak synchrony')
ax.set(title=f'Offset = {offset} frames\nS1 leads <> S2 leads',ylim=[.1,.31],xlim=[0,300], xlabel='Offset',ylabel='Pearson r')
ax.set_xticklabels([int(item-150) for item in ax.get_xticks()])
plt.legend()
複製代碼

上圖中，咱們能夠從負座標推斷出，Subject 1（S1）信號在引導信號間的相互做用（當 S2 被推動了 47 幀的時候相關性最高）。可是，這個評估信號在全局層面會動態變化，例如在這三分鐘內做爲引導信號的信號就會如此。另外一方面，咱們認爲信號之間的相互做用也許會波動得更加明顯，信號是引導仍是跟隨，會隨着時間而轉換。

爲了評估粒度更細的動態變化，咱們能夠計算加窗的時間滯後互相關（WTLCC）。這個過程會在信號的多個時間窗內反覆計算時間滯後互相關。而後咱們能夠分析每一個窗口或者取窗口上的總和，來提供比較二者之間領導者跟隨者互動性差別的評分。

# 加窗的時間滯後互相關
seconds = 5
fps = 30
no_splits = 20
samples_per_split = df.shape[0]/no_splits
rss=[]
for t in range(0, no_splits):
    d1 = df['S1_Joy'].loc[(t)*samples_per_split:(t+1)*samples_per_split]
    d2 = df['S2_Joy'].loc[(t)*samples_per_split:(t+1)*samples_per_split]
    rs = [crosscorr(d1,d2, lag) for lag in range(-int(seconds*fps-1),int(seconds*fps))]
    rss.append(rs)
rss = pd.DataFrame(rss)
f,ax = plt.subplots(figsize=(10,5))
sns.heatmap(rss,cmap='RdBu_r',ax=ax)
ax.set(title=f'Windowed Time Lagged Cross Correlation',xlim=[0,300], xlabel='Offset',ylabel='Window epochs')
ax.set_xticklabels([int(item-150) for item in ax.get_xticks()]);

# 滑動窗口時間滯後互相關
seconds = 5
fps = 30
window_size = 300 #樣本
t_start = 0
t_end = t_start + window_size
step_size = 30
rss=[]
while t_end < 5400:
    d1 = df['S1_Joy'].iloc[t_start:t_end]
    d2 = df['S2_Joy'].iloc[t_start:t_end]
    rs = [crosscorr(d1,d2, lag, wrap=False) for lag in range(-int(seconds*fps-1),int(seconds*fps))]
    rss.append(rs)
    t_start = t_start + step_size
    t_end = t_end + step_size
rss = pd.DataFrame(rss)

f,ax = plt.subplots(figsize=(10,10))
sns.heatmap(rss,cmap='RdBu_r',ax=ax)
ax.set(title=f'Rolling Windowed Time Lagged Cross Correlation',xlim=[0,300], xlabel='Offset',ylabel='Epochs')
ax.set_xticklabels([int(item-150) for item in ax.get_xticks()]);
複製代碼

如上圖所示，是將時間序列分割成了 20 個等長的時間段，而後計算每一個時間窗口的互相關。這給了咱們更細粒度的視角來觀察信號的相互做用。例如，在第一個窗口內（第一行），右側的紅色峯值告訴咱們 S2 開始的時候在引導相互做用。可是，在第三或者第四窗口（行），咱們能夠發現 S1 開始更多的引導相互做用。咱們也能夠繼續計算下去，那麼就能夠得出下圖這樣平滑的圖像。

時間滯後互相關和加窗時間滯後互相關是查看兩信號之間更細粒度動態相互做用的很好的方法，例如引導-追隨關係以及它們如何隨時間改變。可是，這樣的對信號的計算的前提是假設事件是同時發生的，而且具備類似的長度，這些內容將會在下一部分涵蓋。

3. 動態時間扭曲 —— 同步長度不一樣的信號

動態時間扭曲（DTW）是一種計算兩信號間路徑的方法，它能最小化兩信號之間的距離。這種方法最大的優點就是他能處理不一樣長度的信號。最初它是爲了進行語言分析而被髮明出來（在這段視頻中你能夠了解更多），DTW 經過計算每一幀對於其餘全部幀的歐幾里得距離，計算出能匹配兩個信號的最小距離。一個缺點就是它沒法處理缺失值，因此若是你的數據點有缺失，你須要提早插入數據。

爲了計算 DTW，咱們將會使用 Python 的 dtw 包，它將可以加速運算。

from dtw import dtw,accelerated_dtw

d1 = df['S1_Joy'].interpolate().values
d2 = df['S2_Joy'].interpolate().values
d, cost_matrix, acc_cost_matrix, path = accelerated_dtw(d1,d2, dist='euclidean')

plt.imshow(acc_cost_matrix.T, origin='lower', cmap='gray', interpolation='nearest')
plt.plot(path[0], path[1], 'w')
plt.xlabel('Subject1')
plt.ylabel('Subject2')
plt.title(f'DTW Minimum Path with minimum distance: {np.round(d,2)}')
plt.show()
複製代碼

如圖所示咱們能夠看到白色凸形線繪製出的最短距離。換句話說，較早的 Subject2 數據和較晚的 Subject1 數據的同步性可以匹配。最短路徑代價是 d=.33，能夠用來和其餘信號的該值作比較。

4. 瞬時相位同步

最後，若是你有一段時間序列數據，你認爲它可能有振盪特性（例如 EEG 和 fMRI），此時你也能夠測量瞬時相位同步。它也能夠計算兩個信號間每一時刻的同步性。這個結果可能會因人而異由於你須要過濾數據以得到你感興趣的波長信號，可是你可能只有未經實踐的某些緣由來肯定這些波段。爲了計算相位同步性，咱們須要提取信號的相位，這能夠經過使用希爾伯特變換來完成，希爾波特變換會將信號的相位和能量拆分開（你能夠在這裏學習更多關於希爾伯特變換的知識）。這讓咱們可以評估兩個信號是否同相位（兩個信號一塊兒加強或減弱）。

from scipy.signal import hilbert, butter, filtfilt
from scipy.fftpack import fft,fftfreq,rfft,irfft,ifft
import numpy as np
import seaborn as sns
import pandas as pd
import scipy.stats as stats
def butter_bandpass(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a


def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    y = filtfilt(b, a, data)
    return y

lowcut  = .01
highcut = .5
fs = 30.
order = 1
d1 = df['S1_Joy'].interpolate().values
d2 = df['S2_Joy'].interpolate().values
y1 = butter_bandpass_filter(d1,lowcut=lowcut,highcut=highcut,fs=fs,order=order)
y2 = butter_bandpass_filter(d2,lowcut=lowcut,highcut=highcut,fs=fs,order=order)

al1 = np.angle(hilbert(y1),deg=False)
al2 = np.angle(hilbert(y2),deg=False)
phase_synchrony = 1-np.sin(np.abs(al1-al2)/2)
N = len(al1)

# 繪製結果
f,ax = plt.subplots(3,1,figsize=(14,7),sharex=True)
ax[0].plot(y1,color='r',label='y1')
ax[0].plot(y2,color='b',label='y2')
ax[0].legend(bbox_to_anchor=(0., 1.02, 1., .102),ncol=2)
ax[0].set(xlim=[0,N], title='Filtered Timeseries Data')
ax[1].plot(al1,color='r')
ax[1].plot(al2,color='b')
ax[1].set(ylabel='Angle',title='Angle at each Timepoint',xlim=[0,N])
phase_synchrony = 1-np.sin(np.abs(al1-al2)/2)
ax[2].plot(phase_synchrony)
ax[2].set(ylim=[0,1.1],xlim=[0,N],title='Instantaneous Phase Synchrony',xlabel='Time',ylabel='Phase Synchrony')
plt.tight_layout()
plt.show()
複製代碼

瞬時相位同步測算是計算兩個信號每一刻同步性的很好的方法，而且它不須要咱們像計算滑動窗口相關性那樣任意規定窗口寬度。若是你想要知道瞬時相位同步和窗口相關性的比對，能夠在這裏查看我更早些的博客。

---

總結

咱們講解了四種計算時間序列數據相關性的方法：皮爾遜相關，時間滯後互相關，動態時間扭曲及瞬時相位同步。基於你的信號類型，你對信號做出的假設，以及你想要從數據中尋找什麼樣的同步性數據的目標，來決定使用那種相關性測量，有任何問題均可以向我提出，並歡迎在下方留言。

完整的代碼在 Jupyter 筆記上，它使用的樣本數據在這裏。

若是發現譯文存在錯誤或其餘須要改進的地方，歡迎到 掘金翻譯計劃 對譯文進行修改並 PR，也可得到相應獎勵積分。文章開頭的 本文永久連接 即爲本文在 GitHub 上的 MarkDown 連接。

---

掘金翻譯計劃 是一個翻譯優質互聯網技術文章的社區，文章來源爲 掘金 上的英文分享文章。內容覆蓋 Android、iOS、前端、後端、區塊鏈、產品、設計、人工智能等領域，想要查看更多優質譯文請持續關注 掘金翻譯計劃、官方微博、知乎專欄。