FUNC-RETURN

システム、インフラ、AI、デザイン

【Python】EEG(脳波)をカオス時系列解析、takensの埋め込み定理「2章」





ターケンスの埋め込み定理(Takens' Embedding Theorem)


m 次元の力学系の1つの変数 y(t) の時系列データから
適当な時間隔 ∆t をとり作られた k 次元の時間遅れベクトル

\vec{v}(t) = (y(t), y(t + \Delta t), ..., y(t + (k − 1) \Delta t))

が k > 2m であるならば時間遅れベクトルから再構成されたアトラクターは k 次元の空間
に埋め込まれる



ターケンスの埋め込み定理は、力学系における時間遅れ座標を使った埋め込みです。
分かりやすくに言うと、系の一部の情報だけからその系と同じ状態を持った新しい系を再構成できると言うことです。

ちゃんとした説明は、参考を貼ることで省略させてもらいます。


【参考】

多変量解析を用いた力学系の埋め込み次元解析

「埋め込む」ことで大規模複雑システムを理解する

カオス時系列解析とコンピュータ



これを今回は、Pythonを使って実装していきます。




完成したコード

まずは、コード全体をみてみます。

// Takens.py

import numpy as np
from sklearn import metrics


class Takens:

	'''
	constant
	'''
	tau_max = 100


	'''
	initializer
	'''
	def __init__(self, data):
		self.data           = data
		self.tau, self.nmi  = self.__search_tau()


	'''
	reconstruct data by using searched tau
	'''
	def reconstruct(self):
		_data1 = self.data[:-2]
		_data2 = np.roll(self.data, -1 * self.tau)[:-2]
		_data3 = np.roll(self.data, -2 * self.tau)[:-2]
		return _data1, _data2, _data3


	'''
	find tau to use Takens' Embedding Theorem
	'''
	def __search_tau(self):

		# Create a discrete signal from the continunous dynamics
		hist, bin_edges = np.histogram(self.data, bins=200, density=True)
		bin_indices = np.digitize(self.data, bin_edges)
		data_discrete = self.data[bin_indices]

		# find usable time delay via mutual information
		before     = 1
		nmi        = []
		res        = None

		for tau in range(1, self.tau_max):
			unlagged = data_discrete[:-tau]
			lagged = np.roll(data_discrete, -tau)[:-tau]
			nmi.append(metrics.normalized_mutual_info_score(unlagged, lagged))

			if res is None and len(nmi) > 1 and nmi[-2] < nmi[-1]:
				res = tau - 1

		if res is None:
			res = 50

		return res, nmi

参考:node99




再構築に最適な遅延時間\Delta tを見つける

今回は上記の参考先でも使われている、FraserとSwinneyが提案した方法を使っています。

その論文はこちら
Independent coordinates for strange attractors from mutual information

具体的には、遅延時間と非遅延時間の相互情報量が最初に極小となる値を「遅れ」として採用する方法です。
総合情報量は sklearnmetrics.normalized_mutual_info_score() で求めることができます。


コードでは、この部分。

def __search_tau(self):

		# 適度な周波数ビンに離散化
		hist, bin_edges = np.histogram(self.data, bins=200, density=True)
		bin_indices = np.digitize(self.data, bin_edges)
		data_discrete = self.data[bin_indices]

		before     = 1
		nmi        = []
		res        = None

                # self.tau_maxは100
                # 0から100まで相互情報量を出していき、最初の極小値をresに代入。
		for tau in range(1, self.tau_max):
			unlagged = data_discrete[:-tau]
			lagged = np.roll(data_discrete, -tau)[:-tau]
			nmi.append(metrics.normalized_mutual_info_score(unlagged, lagged))

			if res is None and len(nmi) > 1 and nmi[-2] < nmi[-1]:
				res = tau - 1

		if res is None:
			res = 50

		return res, nmi


f:id:kijnmb-1110:20170718154249p:plain




求めた遅延時間から状態空間を再構築

遅延時間\Delta tがもとまると、

\vec{v}(t) = (y(t), y(t + \Delta t), ..., y(t + (k − 1) \Delta t))

によって、状態空間を再構築できる。

ただし、ターケンスの埋め込みの定理はあくまで十分条件なので、必ずしもこの定理を満たしていないと状態空間の再構築ができないと言うわけではない。
実際に、今回の場合も(よく行われている手法)元の3次元の状態空間を、
別の3次元空間に埋め込んでいるので、いわゆる k > 2m この部分を満たしていないが、
視覚的に、元の位相的構造を保持していると分かるので、時間遅れ座標系を使うのが有効と分かる。


f:id:kijnmb-1110:20170718171130p:plain



コードではこの部分。

# 状態空間を再構築
# 求まったτ(tau) × n (n = 0, 1, 2) 分だけ、np.roll()でシフトしてる
def reconstruct(self):
		_data1 = self.data[:-2]
		_data2 = np.roll(self.data, -1 * self.tau)[:-2]
		_data3 = np.roll(self.data, -2 * self.tau)[:-2]
		return _data1, _data2, _data3




次回、バンドパスフィルタについてです。
haretoke.hatenablog.jp


では、今回はターケンスの埋め込み定理でした。


f:id:kijnmb-1110:20170718164100p:plain