Оригинальный адрес:Four ways to quantify synchrony between time series data

Оригинальный автор:Jin Hyun Cheong

Перевод с:Программа перевода самородков

Постоянная ссылка на эту статью:GitHub.com/rare earth/gold-no…

Переводчик:EmilyQiRabbit

Корректор:zhmhhu

Четыре метода количественной синхронизации между данными временных рядов

Примеры кода и данных для расчета метрик синхронизации включают: корреляцию Пирсона, взаимную корреляцию с временным запаздыванием, динамическую деформацию времени и мгновенную фазовую синхронизацию.

Airplanes flying in synchrony, photo by [Gabriel Gusmao](https://unsplash.com/@gcsgpp?utm_source=medium&utm_medium=referral) on [Unsplash](https://unsplash.com?utm_source=medium&utm_medium=referral)

В психологии межличностная синхрония является важным источником информации о социальной динамике и потенциальных социальных результатах. Он может проявляться во многих областях, в том числе в движениях тела (Ramseyer & Tschacher, 2011),Выражение лица(Riehle, Kempkensteffen, & Lincoln, 2017), расширение зрачка (Kang & Wheatley, 2015) и нейронные сигналы (Stephens, Silbert, & Hasson, 2010). так или иначе,синхронностьМогут быть предоставлены различные значения, и существует множество способов количественной оценки синхронизации двух сигналов.

В этой статье я рассматриваю плюсы и минусы некоторых из наиболее часто используемых показателей синхронизации, включая: корреляцию Пирсона, взаимную корреляцию с временной задержкой (TLCC) и оконную TLCC, динамическую деформацию времени и мгновенную фазовую синхронизацию. Для лучшей иллюстрации метрики синхронизации рассчитываются с использованием выборочных данных, представляющих собой трехминутное видео разговора двух участников, из которого мы извлечем улыбки на лицах (изображение ниже — скриншот). Чтобы вы могли лучше ориентироваться в содержании статьи, вы можете скачать ее бесплатноданные лица, извлеченные из образцови включает в себя весь пример кодаЗаметки Юпитера.

содержание

Корреляции Пирсона
Взаимная корреляция с временной задержкой (TLCC) и оконная TLCC
Динамическое искажение времени (DTW)
Мгновенная фазовая синхронизация

Sample data is the smiling facial expression between two participants having a conversation.

1. Корреляция Пирсона — самый простой и лучший метод

Корреляции ПирсонаМера того, как два непрерывных сигнала изменяются вместе во времени, и линейная зависимость между ними может быть выражена как числа -1 (отрицательно коррелированные), 0 (некоррелированные) и 1 (идеально коррелированные). Это интуитивно понятно, легко понять и легко объяснить. Но при использовании корреляции Пирсона необходимо помнить о двух вещах: во-первых, аномальные данные могут мешать результатам оценки корреляции; во-вторых, предполагается, что все данныегомоскедастичность, чтобы дисперсия данных была однородной по всему диапазону данных. Как правило, корреляция является моментальной мерой глобальной синхронности. Следовательно, он не может предоставить информацию о направленности между двумя сигналами, например, какой сигнал является ведущим сигналом, а какой является последующим сигналом.

Многие пакеты используют корреляцию Пирсона, включая Numpy, Scipy и Pandas. Если ваши данные содержат нулевые или отсутствующие значения, метод корреляции в Pandas отбросит эти строки перед расчетом, и если вы хотите использовать приложения корреляции Numpy или Scipy для Pearson, вам придется вручную очистить их. Эти данные.

Следующий код загружает образцы данных (они находятся в той же папке, что и код), использует 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")

Sample data on top, moment-to-moment synchrony from moving window correlation on bottom.

В целом, корреляция Пирсона является хорошим вводным учебным пособием, она предоставляет очень простой способ расчета глобальной и локальной синхронизации. Однако он не предоставляет информацию о динамике сигнала, например, какой сигнал появился первым, что можно измерить с помощью взаимной корреляции.

2. Взаимная корреляция с временным запаздыванием — оценка динамики сигнала

Взаимная корреляция с запаздыванием по времени (TLCC) может определять направленность между двумя сигналами, например, отношение ведущего-следующего, в котором сигнал ведущего инициирует ответ, а сигнал ведомого повторяет его. Существуют и другие способы исследовать этот тип отношений, в том числеПричинность Грейнджера, которые часто используются в экономике, но имейте в виду, что они не обязательно отражают истинную причинно-следственную связь. Однако, глядя на взаимную корреляцию, мы все же можем извлечь информацию о том, какой сигнал появился первым.

[http://robosub.eecs.wsu.edu/wiki/ee/hydrophones/start](http://robosub.eecs.wsu.edu/wiki/ee/hydrophones/start)

Как показано на рисунке выше, TLCC измеряется путем постепенного смещения вектора временного ряда (красная линия) и многократного вычисления корреляции между двумя сигналами. Если пик корреляции находится в центре (смещение = 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()

Peak synchrony is not at the center, suggesting a leader-follower signal dynamic.

На приведенном выше рисунке мы можем сделать вывод из отрицательных координат, что сигнал Субъекта 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()]);

Windowed time lagged cross correlation for discrete windows

Как показано на рисунке выше, временной ряд делится на 20 временных периодов одинаковой длины, а затем вычисляется взаимная корреляция каждого временного окна. Это дает нам более детальное представление о взаимодействии сигналов. Например, в первом окне (верхний ряд) красный пик справа говорит нам о том, что S2 запускает взаимодействие в самом начале. Однако в третьем или четвертом окне (строке) мы можем обнаружить, что S1 запускает больше самозагрузочных взаимодействий. Мы также можем продолжить вычисления, тогда мы сможем получить гладкое изображение, подобное изображенному ниже.

Rolling window time lagged cross correlation for continuous windows

Взаимная корреляция с временной задержкой и оконная взаимная корреляция с временной задержкой — отличные способы увидеть более детализированные динамические взаимодействия между двумя сигналами, такие как отношения опережения-следования и их изменение во времени. Однако такие вычисления сигналов предполагают, что события происходят одновременно и имеют одинаковую продолжительность, что будет рассмотрено в следующем разделе.

3. Динамическое искажение времени — сигналы с разной длительностью синхронизации.

Dynamic Time Warping (DTW) — это метод расчета пути между двумя сигналами, который минимизирует расстояние между ними. Самым большим преимуществом этого метода является то, что он может обрабатывать сигналы разной длины. Первоначально он был изобретен для лингвистического анализа (вэто видеоВы можете узнать больше в ), DTW вычисляет наименьшее расстояние, которое соответствует двум сигналам, вычисляя евклидово расстояние каждого кадра для всех остальных кадров. Одним из недостатков является то, что он не может обрабатывать пропущенные значения, поэтому, если ваши точки данных отсутствуют, вам нужно будет интерполировать данные заранее.

XantaCross [CC BY-SA 3.0 ([https://creativecommons.org/licenses/by-sa/3.0](https://creativecommons.org/licenses/by-sa/3.0))]

Для расчета DTW мы будем использовать Pythondtwпакет, он сможет ускорить работу.

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. Мгновенная фазовая синхронизация

Наконец, если у вас есть данные временного ряда, которые, по вашему мнению, могут иметь колебательные свойства (например, ЭЭГ и фМРТ), вы также можете измерить мгновенную фазовую синхронизацию в этот момент. Он также может вычислять синхронизацию между двумя сигналами в каждый момент времени. Этот результат может варьироваться от человека к человеку, поскольку вам необходимо отфильтровать данные, чтобы получить сигнал на интересующих вас длинах волн, но у вас могут быть только непрактичные причины для идентификации этих диапазонов. Чтобы вычислить фазовую синхронизацию, нам нужно извлечь фазу сигнала, что можно сделать с помощью преобразования Гильберта, которое разделяет фазу и энергию сигнала (Вы можете узнать больше о преобразовании Гильберта здесь). Это позволяет нам оценить, находятся ли два сигнала в фазе (оба сигнала увеличиваются или уменьшаются вместе).

Gonfer at English Wikipedia [CC BY-SA 3.0 ([https://creativecommons.org/licenses/by-sa/3.0](https://creativecommons.org/licenses/by-sa/3.0))]

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()

Filtered time series (top), angle of each signal at each moment in time (middle row), and instantaneous phase synchrony measure (bottom).

Мгновенная оценка фазовой синхронизации — это хороший способ вычислить моментальную синхронизацию двух сигналов, и он не требует от нас произвольного указания ширины окна, как мы делаем для корреляций со скользящим окном. Если вы хотите узнать сравнение мгновенной фазовой синхронизации и корреляции окна,Проверьте мой предыдущий блог здесь.

Суммировать

Мы описываем четыре метода расчета корреляции данных временных рядов: корреляция Пирсона, взаимная корреляция с временным запаздыванием, динамическая деформация времени и мгновенная фазовая синхронизация. Решите, какую меру корреляции использовать на основе вашего типа сигнала, предположений, которые вы делаете о своем сигнале, и какие цели синхронности данных вы хотите искать в своих данных, не стесняйтесь задавать мне любые вопросы и добро пожаловать Оставьте комментарий ниже.

Полный код находится в Jupyter Notes, в котором используетсяПример данных здесь.

Если вы обнаружите ошибки в переводе или в других областях, требующих доработки, добро пожаловать наПрограмма перевода самородковВы также можете получить соответствующие бонусные баллы за доработку перевода и PR. начало статьиПостоянная ссылка на эту статьюЭто ссылка MarkDown этой статьи на GitHub.

Программа перевода самородковэто сообщество, которое переводит высококачественные технические статьи из Интернета сНаггетсДелитесь статьями на английском языке на . Охват контентаAndroid,iOS,внешний интерфейс,задняя часть,блокчейн,продукт,дизайн,искусственный интеллектЕсли вы хотите видеть более качественные переводы, пожалуйста, продолжайте обращать вниманиеПрограмма перевода самородков,официальный Вейбо,Знай колонку.