Step0. パッケージを用意

import numpy as np
import matplotlib.pyplot as plt
import random

Step1. 写像を定義

写像を定義

def f(x, alpha):
	y = 1 - alpha*x*x
	return y

def coupled_f(x_list, alpha, e):
	total = 0
	y_list = []
	num = len(x_list)

	for n in range(num):
		total += f(x_list[n], alpha)
	for n in range(num):
		y_list.append((1 - e) * f(x_list[n], alpha) + (e / num) * total)
	return y_list

Step2. 時系列データを生成しそれぞれの差に変換

時系列データを生成

def series(alpha, e, initial, length, trans_period):
	num_coupled = len(initial)  # 結合数
	x_list = [n + trans_period for n in range(length)]
	y_list = []

	for n in range(num_coupled):
		y_list.append([initial[n]])
	length -= 1
	for n1 in range(length + trans_period):
		pre_y_list = []
		for n2 in range(len(initial)):
			pre_y_list.append(y_list[n2][n1])

		after_y_list = coupled_f(pre_y_list, alpha, e)
		for n2 in range(len(initial)):
			y_list[n2].append(after_y_list[n2])

	for n1 in range(num_coupled):
		y_list[n1] = y_list[n1][trans_period:]
	return x_list, y_list

それぞれの差に変換

def product_dif_list(y_list):
	num_coupled = len(y_list[0])
	dif_list = []
	for n1 in range(len(y_list)):
		sub_list = y_list[n1]
		sub_dif_list = []
		for n2 in range(num_coupled):
			for n3 in range(n2 + 1, num_coupled):
				sub_dif_list.append(abs(sub_list[n2] - sub_list[n3]))
		dif_list.append(sub_dif_list)
	return dif_list

Step3. 確率分布のラベルを作成

ラベルを作成

def product_label(y_min, y_max, num_label):
	y_length = (y_max - y_min) / num_label
	x_list = []
	for n1 in range(num_label):
		y_min += y_length
		x_list.append(round(y_min, 4))
	return x_list

システム全体のラベルを作成

# 2進展開
def Dec_to_N(num, base):
	if num >= base:
			yield from Dec_to_N(num // base, base)
	yield num % base

# 状態空間を判別
def state_replace(state_matrix, convert_list):
	for i in range(len(state_matrix)):
		for j in range(len(state_matrix[i])):
			state_matrix[i][j] = convert_list[state_matrix[i][j]]
	return state_matrix

# 全体の状態ラベルの行列を作成
def product_whole_label_list(label_list, num_coupled):
	num_label = len(label_list) # number of state
	whole = list(range(num_label ** num_coupled))
	state_matrix = [list(Dec_to_N(n, num_label)) for n in whole]
	for n1 in range(len(state_matrix)):
			while len(state_matrix[n1]) < num_coupled:
					state_matrix[n1].insert(0, 0)
	whole_label_list = state_replace(state_matrix, label_list)
	return whole_label_list

Step4. 同時確率行列を生成

全体の確率行列を作成

def product_whole_prob_matrix(dif_list, dif_label):
	couple    = len(dif_list)
	whole_porb = [[0 for n in range(len(dif_label) ** couple)] for n in range(len(dif_label) ** couple)]
	for n1 in range(int(len(dif_list[0]) - 1)):
		present = [dif_list[n2][n1] for n2 in range(len(dif_list))]
		future  = [dif_list[n2][n1 + 1] for n2 in range(len(dif_list))]

		flag = False
		for n2 in range(len(whole_porb)):
			test_index = 0
			for n3 in range(couple):
				if present[n3] <= dif_label[n2]:
					row = n2
					test_index += 1
					if test_index == couple:
						flag = True
						break
			if flag:
				break

		flag = False
		for n2 in range(len(whole_porb)):
			test_index = 0
			for n3 in range(couple):
				if future[n3] <= dif_label[n2]:
					column = n2
					test_index += 1
					if test_index == couple:
						flag = True
						break
			if flag:
				break
		whole_porb[row][column] += 1
	total_sum = sum(sum(row) for row in whole_porb)
	joint_prob_matrix = [[element / total_sum for element in row] for row in whole_porb]
	return joint_prob_matrix

Appendix. 個別の要素の同時確率行列を作成

個別の要素の同時確率行列を作成

def product_prob_matrix(y_list, label_list):
	num_label = len(label_list)
	length = int(len(y_list) - 1)
	prob_matrix = []
	for n1 in range(num_label):
		prob_matrix.append([0 for n2 in range(num_label)])
	for n1 in range(length):
		present = y_list[n1]
		future = y_list[n1 + 1]
		for n2 in range(num_label):
			if present <= label_list[n2]:
				row = n2
				break
		for n2 in range(num_label):
			if future <= label_list[n2]:
				column = n2
				break
		prob_matrix[row][column] += 1
	total_sum = sum(sum(row) for row in prob_matrix)
	joint_prob_matrix = [[element / total_sum for element in row] for row in prob_matrix]
	return joint_prob_matrix

同時確率行列から遷移確率行列に変換

def product_tpm(whole_prob):
	tpm = []
	present = np.sum(np.array(whole_prob), axis=1)
	for row_idx in range(len(whole_prob)):
		tpm.append([])
		for col_idx in range(len(whole_prob)):
			if present[row_idx] > 0:
				tpm[row_idx].append(whole_prob[row_idx][col_idx] / present[row_idx])
			else:
				tpm[row_idx].append(0)
	return tpm