獨立成分分析(ICA)的模擬實驗(R語言)

本文永久連接: https://esl.hohoweiya.xyz/notes/ICA/index.html

本筆記是ESL14.7節圖14.42的模擬過程。第一部分將以ProDenICA法爲例試圖介紹ICA的整個計算過程；第二部分將比較ProDenICA、FastICA以及KernelICA這種方法，試圖重現圖14.42。

ICA的模擬過程

生成數據

首先咱們得有一組獨立（ICA的前提條件）分佈的數據$S$（未知），而後通過矩陣$A_0$混合以後獲得實際的觀測值$X$，即

$$ X= SA_0 $$

也能夠寫成 $$ S=XA_0^{-1} $$

用雞尾酒酒會的例子來講就是，來自不一樣個體的說話聲通過麥克風混合以後獲得咱們實際接收到的信號。假設有兩組獨立同分布的數據，分佈都爲n（對應圖14.42中的編號），每組數據個數均爲$N=1024$，混合矩陣爲A0，用R代碼描述這一過程以下

library(ProDenICA)
p = 2
dist = "n" 
N = 1024
A0 = mixmat(p)
s = scale(cbind(rjordan(dist,N),rjordan(dist,N)))
x = s %*% A0

最終咱們獲得觀測值x。

白化

在進行ICA時，也就是恢復$X=S\mathbf A$中的混合矩陣$\mathbf A$，都會假設$X$已經白化獲得$\mathrm{Cov}(X)=\mathbf I$，而這個處理過程能夠用SVD實現。對於中心化的$X$，根據

$$ X=\mathbf{UDV}^T= \sqrt{N}\mathbf U\frac{1}{\sqrt{N}}\mathbf{DV}^T=X^*\frac{1}{\sqrt{N}}\mathbf{DV}^T $$

獲得知足$Cov(X^)=\mathbf I$的$X^$，則

$$ S=XA_0^{-1}=X^*DV^TA_0^{-1}/\sqrt{N} $$

因而通過這個變換以後，混合矩陣變爲

$$ A = DV^TA_0^{-1}/\sqrt{N} $$

則 $$ X^*=SA^T $$

用R語言表示以下

x <- scale(x, TRUE, FALSE) # central
sx <- svd(x)	
x <- sqrt(N) * sx$u # satisfy cov(x) = I
target <- solve(A0)
target <- diag(sx$d) %*% t(sx$v) %*% target/sqrt(N) # new mixing maxtrix

ProDenICA法

細節再也不展開，直接利用ProDenICA中的包進行計算

W0 <- matrix(rnorm(2*2), 2, 2)
W0 <- ICAorthW(W0)
W1 <- ProDenICA(x, W0=W0,trace=TRUE,Gfunc=GPois)$W

獲得$A$的估計值W1

計算Amari距離

amari(W1, target)

比較兩種算法

這一部分試圖重現Fig. 14.42。

N = 1024

genData <- function(dist, N = 1024, p = 2){
  # original sources
  s = scale(cbind(rjordan(dist, N), rjordan(dist, N)))
  # mixing matrix 
  mix.mat = mixmat(2)
  # original observation
  x = s %*% mix.mat
  
  # central x
  x = scale(x, TRUE, FALSE)
  # whiten x
  xs = svd(x)
  x = sqrt(N) * xs$u # new observations
  mix.mat2 = diag(xs$d) %*% t(xs$v) %*% solve(mix.mat) / sqrt(N) # new mixing matrix
  return(list(x = x, A = mix.mat2))
}

res = array(NA, c(2, 18, 30))
for (i in c(1:18)){
  for (j in c(1:30)){
    data = genData(letters[i])
    x = data$x
    A = data$A
    W0 <- matrix(rnorm(2*2), 2, 2)
    W0 <- ICAorthW(W0)
    # ProDenICA
    W1 <- ProDenICA(x, W0=W0,trace=FALSE,Gfunc=GPois, restarts = 5)$W
    # FastICA
    W2 <- ProDenICA(x, W0=W0,trace=FALSE,Gfunc=G1, restarts = 5)$W
    res[1, i, j] = amari(W1, A)
    res[2, i, j] = amari(W2, A)  
  }
}

res.mean = apply(res, c(1,2), mean)
#offset = apply(res, c(1,2), sd)
#offset = 0
#res.max = res.mean + offset/4
#res.min = res.mean - offset/4
res.max = apply(res, c(1,2), max)
res.min = apply(res, c(1,2), min)

# plot
plot(1:18, res.mean[1, ], xlab = "Distribution", ylab = "Amari Distance from True A", xaxt = 'n', type = "o", col = "orange", pch = 19, lwd = 2, ylim = c(0, 0.5))
axis(1, at = 1:18, labels = letters[1:18])
lines(1:18, res.mean[2, ], type = "o", col = 'blue', pch = 19, lwd=2)
legend("topright", c("ProDenICA", "FastICA"), lwd = 2, pch = 19, col = c("orange", "blue"))
#for(i in 1:18)
#{
#  for (j in 1:2)
#  {
#    color = c("orange", "blue")
#    lines(c(i, i), c(res.min[j, i], res.max[j, i]), col = color[j], pch = 3)
#    lines(c(i-0.2, i+0.2), c(res.min[j, i], res.min[j, i]), col = color[j], pch = 3)
#    lines(c(i-0.2, i+0.2), c(res.max[j, i], res.max[j, i]), col = color[j], pch = 3)
#  }
#}

獲得下圖

與圖14.42的右圖中的FastICA和ProDenICA的圖象一致。

試圖繪製出圖中的變化範圍，但因爲書中並未指出變換範圍是什麼，嘗試了標準差及最大最小值，但效果不是很好，這是能夠繼續優化的一個方面。下圖是用四分之一的標準差做爲其波動範圍獲得的

待完善

加入kernelICA