Python での相互相関

相互相関は、一方の数値列をもう一方に対してずらしながら、2つの実数列を比較します。相対位置ごとにスコアが得られるため、パターンがどの位置で一致するかを調べられます。既定のスコアは積の合計であり、自動的に正規化された相関係数になるわけではありません。

相互相関で測定するもの

実数列 x と y の相互相関では、重なり合う値同士を掛け、その積を各シフト位置で合計します。正のスコアが高いほど、その位置では値の符号やパターンが似ている傾向があります。負のスコアは、符号が反対であることを示します。スコアは列の尺度と、重なり合う値の個数にも左右されます。

full モードでは、少なくとも1組の値が重なるすべてのシフトを含めます。x の長さが N、y の長さが M の場合、スコアは N + M - 1 個になります。valid モードでは、短い方の列が長い方の列の範囲内に完全に収まるシフトだけを残します。same モードでは、最初の入力を中心とする、最初の入力と同じ長さの結果を返します。長さが偶数の場合、中心位置の厳密な規則はライブラリによって異なります。これらのモードが選ぶのはシフトであり、正規化方法ではありません。

純粋なPythonで相互相関を計算する

直接実装すると、値の重なり方とラグの規則が明確になります。実数列 x と y に対して、次のコードは重なり合う値の積を合計した、正規化されていないスコアを返します。実数入力では、ラグの順序と結果はNumPyの correlate(x, y, mode="full") の規則と一致します。

次の例を実行してください。

このセクションのコードと表示されるラベルや出力は英語のままです。

x = [1, 2, 3]
y = [4, 5]

scores = []
for lag in range(-len(y) + 1, len(x)):
    score = sum(
        x[index + lag] * y[index]
        for index in range(len(y))
        if 0 <= index + lag < len(x)
    )
    scores.append(score)

print(scores)

[5, 14, 23, 12]

この列にはラグごとに1つのスコアが含まれ、-(len(y) - 1) から len(x) - 1 までの順に並びます。各スコアには、インデックスが重なる組み合わせだけが含まれます。この単純な方法は、小さな入力や計算の仕組みを理解するのに適していますが、入れ子になった処理の量は2つの列の長さの積に応じて増えます。

NumPyで相互相関を計算する

NumPyは python -m pip install numpy でインストールできます。numpy.correlate は1次元の実数値数値配列を処理し、full、same、valid の各モードを利用できます。既定のモードは valid なので、必要なシフト範囲に合わせてモードを指定してください。

このセクションのコードと表示されるラベルや出力は英語のままです。

import numpy as np

x = np.array([1, 2, 3])
y = np.array([4, 5])

scores = np.correlate(x, y, mode="full")
print(scores.tolist())

[5, 14, 23, 12]

この実数入力では、結果のラグ順序と正規化されていない積の合計は、純粋なPythonの例と同じです。複素数データではNumPyは2番目の入力を複素共役にしますが、上の純粋なPythonコードはそうしません。そのため、両者の一致を主張できるのは実数入力の場合です。通常のPython形式が必要な場合に限り、配列をリストに変換してください。その後も数値処理を続けるなら、一般にNumPy配列の方が便利です。この比較では、入力が実数値の1次元配列である必要があります。

SciPyで相互相関を計算する

SciPyは python -m pip install scipy でインストールできます。SciPyにはNumPyも必要です。scipy.signal.correlate は同じ3つの出力モードを提供し、SciPyのほかの信号処理ツールと併用できます。correlation_lags を使うと、結果の各位置に対応するラグが返されるため、インデックスを手作業で計算する必要がありません。

このセクションのコードと表示されるラベルや出力は英語のままです。

import numpy as np
from scipy.signal import correlate, correlation_lags

x = np.array([1, 2, 3])
y = np.array([4, 5])

scores = correlate(x, y, mode="full", method="auto")
lags = correlation_lags(len(x), len(y), mode="full")
print(scores.tolist())
print(lags.tolist())

[5, 14, 23, 12]
[-1, 0, 1, 2]

既定では、SciPyの method オプションが計算方法を選びますが、指定した相互相関の定義は変わりません。異なるライブラリの結果を比較するときは、同じ実数入力と同じモードを使ってください。NumPyと同様、SciPyは複素数データでは2番目の入力を複素共役にしますが、上の純粋なPythonコードはそうしません。この例では実数値の1次元配列を使います。

ラグを解釈して方法を選ぶ

実数値の x を第1引数、実数値の y を第2引数として渡す場合、NumPyとSciPyの full モードでは、ラグは -(len(y) - 1) から len(x) - 1 の順に並びます。この規則では、正のラグは x の後方部分と y を比較します。ピークを時間差に読み替える前に、必ずライブラリの引数順とラグの規則を確認してください。入力を入れ替えると、ラグの向きも逆になります。

依存パッケージを使わない基準計算が必要な場合や、列が短い場合は純粋なPythonを選びます。数値処理コードで簡潔な配列演算が必要ならNumPyを選びます。明示的なラグ値やほかの信号処理関数も必要ならSciPyを選びます。サンプル間隔が一定なら、ラグのインデックスにサンプリング間隔を掛けて時間に換算できます。正規化係数が必要な場合は、正規化方法を別途定義し、重なり合うサンプル数も考慮してください。生の相互相関スコアは、重なりの長さや尺度が異なる場合、そのままでは比較できません。

相互相関とPandasの相関

pandas.Series.corr は、観測値の位置を揃えて単一のペアワイズ相関係数を計算します。一方の系列をもう一方に対してシフトしたり、ラグごとのスコアを返したりはしないため、相互相関関数ではありません。ラグを順に調べるにはNumPyまたはSciPyを使います。Pandasのデータを使う場合は、まず意図した方法で位置を揃えて数値を取り出し、選んだ実装に配列として渡してください。インデックスの不一致や欠損値をどう扱うかも、明示的に決める必要があります。