マルチフラクタル次元を計算するため,多次元データの度数分布表を作成するプログラム.
マルチフラクタル次元は,一般化次元(北海道大学の井上 純一先生作成の資料「2009年度 カオス・フラクタル 講義ノート」の第8回や第12回に示される)とも言われる.
以下のプログラムは効率が悪い.
データの特徴量次元を,データ点数を,ビンの分割数をとすると,
メモリ空間はに比例し,計算量(ループ回数)はに比例する.
ソース中のループの実行時間は,約1600秒.ループ一回当たりの実行時間は約0.03秒.
特徴量次元D=2,データ点数N=104847,ビンの分割数V=2 * 3 * 5 * 7とした.
import numpy as np import wave from scipy import fromstring, int16, int64 import itertools import seaborn as sns import matplotlib.pyplot as plt import time def probability_distribution(X, div_num=2 * 3 * 5 * 7): """ 多次元データの度数分布表を作成する. 多次元ヒストグラムあるいは多次元ヒートマップのような,数値の分布を計算する. :param X:入力データ.第0方向が時系列.第1方向が特徴量次元 :param div_num:ヒストグラムで言うところのビンの数 :return: 合計(div_num^D)個の要素を持つ,D次元配列を返す.ここで,Dは特徴量次元数. """ X_max = int64(np.max(X, axis=0)) # 入力データが小数点数ならint64以外にする X_min = int64(np.min(X, axis=0)) # ビン分割の最大値と最小値に余裕を持たせる bin_margin = 0.2 X_max = X_max * (1 + np.sign(X_max) * bin_margin) X_min = X_min * (1 - np.sign(X_min) * bin_margin) bin_resolution = np.array([(X_max - X_min) / div_num]) D = [] # 特徴量次元数 N = [] # データ点数 len_X_shape = X.shape.__len__() if len_X_shape == 2: D = X.shape[1] N = X.shape[0] elif len_X_shape == 1: D = 1 N = X.shape[0] else: return None too_many_D = 30 # 次元が多すぎる場合 if D > too_many_D: import sys print( f"High-dimensional data(data dimension={D}) was input to the function \"{sys._getframe().f_code.co_name}\".\n" f"Do you want to ontinue? [Y]/[N]") if "Y" != input(): print(f"Interrupted was processing of the function \"{sys._getframe().f_code.co_name}.") return None ind_tmp = np.array(list(range(div_num))) ind = np.array(list(itertools.product(ind_tmp, repeat=D))) N_i = np.zeros(np.repeat(div_num, D)) # 度数分布表の配列 start = time.time() # 下のループにかかる時間を測定する. for i in ind: range_min = bin_resolution * i + X_min range_max = bin_resolution * (i + 1) + X_min temp_bin_cond = sum(map(all, np.logical_and(range_min <= X, X < range_max))) N_i[tuple(i)] = temp_bin_cond pass end = time.time() print(f"ループ処理全体に掛かった時間:{(end - start)}") print(f"ループ一周あたりの平均時間:{(end - start) / len(ind)}") sns.heatmap(np.log10(N_i + 1e-10)) # ヒートマップを描画.値の大きさは対数にする. print(f"sum N_i:{np.sum(N_i)} (入力データのデータ数(次元数ではなく)に一致するはず)") return N_i # 度数分布表を返す pass if __name__ == "__main__": file_name = "RIFF wavファイルのパス" with wave.open(file_name, "r") as wf: fs = wf.getframerate() n = wf.getnframes() data = wf.readframes(wf.getnframes()) X = fromstring(data, dtype=int16) # t = np.linspace(0, 10, 1000) # 自分でデータを作る場合 # X = np.sin(2 * np.pi * np.sqrt(101) * t) * 100 # 自分でデータを作る場合 # 特徴量が2次元の時系列データを作成(X[:,0]は時系列データ,X[:,1]は時系列データの時間差分(時間微分)) data = np.concatenate([np.array([X[0:-1]]), np.array([np.diff(X)])], axis=0).T # 上で作成したデータから多次元度数分布表を計算する.ビン分割数は2*3*5*7個としている. probability_distribution(X=data, div_num=2 * 3 * 5 * 7) plt.figure() plt.plot(data[:, 0], data[:, 1]) plt.show() pass