Zenn
🟠

拡張クロスデータ行列法 with Python

2023/11/03に公開
Python
math
統計
統計学
データサイエンス
tech

はじめに

拡張クロスデータ行列法(extended-cross-data-matrix methodology (ECDM))は、CDMのデータ分割をいろんなパターンでするようにした方法です。

CDMについてはこちらをチェック。
https://zenn.dev/totopironote/articles/b63807d5fb8333

今回は分割方法を紹介し、その分割するpythonコードを載せてあります。
ECDMを用いて推定量を求めるところまではしてないので詳しく知りたい方は参考文献をご確認ください。

拡張クロスデータ行列法

n_{(1)}= \lceil n/2 \rceil とし、n_{(2)} = n-n_{(1)}と定義します。

また、以下のように定義します：

\begin{align*} V_{n(1)(k)} &= \begin{cases} \{ \lfloor k/2 \rfloor -n_{(1)}+1, \ldots, \lfloor k/2 \rfloor \} ,& \lfloor k/2 \rfloor \ge n_{(1)}\\ \{ 1, \ldots, \lfloor k/2 \rfloor \} \cup \{ \lfloor k/2 \rfloor +n_{(2)}+1, \ldots, n \} ,& それ以外 \end{cases}\\ V_{n(2)(k)} &= \begin{cases} V_{n(2)(k)} = \{ \lfloor k/2 \rfloor +1, \ldots, \lfloor k/2 \rfloor +n_{(2)} \}, & \lfloor k/2 \rfloor \le n_{(1)}\\ \{1, \ldots, \lfloor k/2 \rfloor -n_{(1)} \} \cup \{ \lfloor k/2 \rfloor +1, \ldots, n \}, & それ以外 \end{cases} \end{align*}

注意: \#V_{n(l)(k)}=n_{(l)}（l=1,2）、V_{n(1)(k)} \cap V_{n(2)(k)}=\emptyset および

V_{n(1)(k)} \cup V_{n(2)(k)}=\{1, \ldots, n\}（k=3, \ldots, 2n-1）です。

また、i<j\ (\le n) の場合、i\in V_{n(1)(i+j)} および j\in V_{n(2)(i+j)} です。

さらに、以下のように定義します：k=3, \ldots, 2n-1 に対して

\overline{x}_{n{(1)}(k)}=n_{(1)}^{-1} \sum_{j \in V_{n{(1)}(k)}} x_j \hspace{5pt}, \hspace{5pt}\overline{x}_{n(2)(k)}=n_{(2)}^{-1} \sum_{j\in V_{n(2)(k)}} x_j

ECDMを用いてtr(\Sigma^2)を推定したりしている。

V_{n(1)(i+j)} と V_{n(2)(i+j)} に対応してデータ集合を2分割し、それらに基づいて不偏推定量を1つ計算する。これを， i < j (≤ n)のすべての組合せで繰り返し，得られる不偏推定量の平均をとる。

Pythonコード

def ECDM(X):
  d, n = X.shape
  n1 = np.ceil(n / 2).astype(int)
  n2 = n - n1

  data_n1 = []
  data_n2 = []

  for k in range(3, 2 * n):
          if np.floor(k / 2) >= n1:
              v_n1_k = set(range(int(np.floor(k / 2) - n1 + 1), int(np.floor(k / 2) + 1)))
          else:
              v_n1_k = set(range(1, int(np.floor(k / 2) + 1))).union(set(range(int(np.floor(k / 2) + n2 + 1), n + 1)))
          
          if np.floor(k / 2) <= n1:
              v_n2_k = set(range(int(np.floor(k / 2) + 1), int(np.floor(k / 2) + n2 + 1)))
          else:
              v_n2_k = set(range(1, int(np.floor(k / 2) - n1 + 1))).union(set(range(int(np.floor(k / 2) + 1), n + 1)))

          a = np.array(list(v_n1_k))-1
          b = np.array(list(v_n2_k))-1

          data_n1.append(X[:, a])
          data_n2.append(X[:, b])

  return np.array(data_n1), np.array(data_n2)

参考文献

Discussion