独立成分分析の実装

独立成分分析とは、ブラインド信号源分離に用いられる。これは、2つの異なる音が混合した音から、２種類の音を判別するものである。
今回は、その手法の一つの射影追跡を用いて実装する。
射影追跡とは、2つの要素からより非ガウスの部分を抽出することである。これは、中心極限定理より、最も独立した音を抽出することと、非ガウスを見つけることは、等しい。

仮定として、まず音を聞き取る装置は2つあるものとする。音の大きさやその順序は問わない。

#モチベーション
課題

#アルゴリズム
射影追跡
ニュートン近似法

#乱数生成

import numpy as np
from scipy.linalg import sqrtm
import math
import matplotlib.pyplot as plt


def generate_data(n=1000):
    global M, b
    x = np.concatenate([np.random.rand(n, 1), np.random.randn(n, 1)], axis=1)
    x[0, 1] = 6   # outlier
    x = (x - np.mean(x, axis=0)) / np.std(x, axis=0)  # Standardization
    M = np.array([[1, 3], [5, 3]])
    x = x.dot(M.T)
    x = np.linalg.inv(sqrtm(np.cov(x, rowvar=False))).dot(x.T).T
    b  = np.random.uniform(0.08, -0.08, size= (2, 1))
    return x

X = generate_data()

for i in range(X.shape[0]):
  plt.scatter(generate_data()[i][0], generate_data()[i][1], color= 'red')
plt.show()

すでに上で示されている行列Mで音声が混合されて居る。

#導入する関数の定義

def g_s_3_dif(s):
  return 3 * (s**2)

def g_tan_dif (s):
  return 1 - (math.tanh(s))**2

def g_s_3 (s):
  return s**3

def g_tan (s):
  return math.tanh(s)

#球状化

#xの球状か
def kyujouka ():
  X_= []
  for i in range(X.shape[0]):
    x0, x1 = 0, 0
    n =X.shape[0]
    for j in range(X.shape[0]):
      x0 += X[j][1]
      x1 += X[j][1]
    x_ = (X[i] - [x0/n, x1/n])
    X_.append(x_.tolist())
  X_ = np.array(X_)
#   X_.append(x_)
  global X__
  X__ = []
  for l in range(X_.shape[0]):
    a = np.matmul(X_[l], X_[l].T)/n
    aa = 1/math.sqrt(a)
    x__ = aa*X_[l]
    X__.append(x__)
  return X__
X__ = kyujouka()
print((X__))

#尖度を用いたとき

difference = 1
b  = np.random.uniform(0.08, -0.08, size= (1, 2))
print(b)
X__ = np.array(kyujouka())
while difference > 0.1:
  n = X.shape[0]
  sum_dif, sum_ = 0, 0
  for i in range(X.shape[0]):
    a = np.matmul(b, X__[i].T)
    g = g_s_3(a[0])
    gdif = g_s_3_dif(a[0])
    sum_dif += gdif
    sum_ += X__[i]*g
  b_before = b
  b = (sum_dif/n)*b - (sum_/n)
  b = b/math.sqrt(np.matmul(b, b.T)[0][0])
  difference = abs((b[0][0]) - abs(b_before[0][0]))
y =[]
for i in range(X.shape[0]):
  a = np.matmul(b, X[i].T)
  y.append(a[0])

aa = np.histogram(y, bins = 50)
a_bins = aa[1]
a_hist = aa[0]
X1 = []
for i in range(1, len(a_bins)):
    X1.append((a_bins[i-1]+a_bins[i])/2)
plt.bar(X1,a_hist, width=0.08)

ガウスに近く、射影追跡がうまくいっていない。

#tanhを用いたとき

difference = 1
b  = np.random.uniform(0.08, -0.08, size= (1, 2))
print(b)
X__ = np.array(kyujouka())
while difference > 0.1:
  n = X.shape[0]
  sum_dif, sum_ = 0, 0
  for i in range(X.shape[0]):
    a = np.matmul(b, X__[i].T)
    g = g_tan(a[0])
    gdif = g_tan_dif(a[0])
    sum_dif += gdif
    sum_ += X__[i]*g
  b_before = b
  b = (sum_dif/n)*b - (sum_/n)
  print(np.matmul(b, b.T)[0][0], b)
  b = b/math.sqrt(np.matmul(b, b.T)[0][0])
  print(b)
  difference = abs((b[0][0]) - abs(b_before[0][0]))

  print(difference)

y =[]
for i in range(X.shape[0]):
  a = np.matmul(b, X[i].T)
  y.append(a[0])
aa = np.histogram(y, bins =50)
print(aa)
a_bins = aa[1]
a_hist = aa[0]
X1 = []
for i in range(1, len(a_bins)):
    X1.append((a_bins[i-1]+a_bins[i])/2)
plt.bar(X1,a_hist, width=0.05)

うまく分離できて居る