negentropyICA: Fixed-point algorithm based on negentropy

In hanspeter6/icaPlay: What the Package Does (One Line, Title Case)

Description

A deflationary fixed-point algorithm to estimate an unmixing matrix based on negentropy.

Usage

negentropyICA(x, g = g2, max.iters = 10, init.w = diag(m), tol = 1e-04)

Arguments

x	A mixture set with variables in rows and observations in columns.
g	One of g1, g2 or g3. Default = g2. See details for more information.
max.iters	The maximum number of iterations. Default = 10.
init.w	The initial unmixing matrix. Default is identity.
tol	The tolerance for theta (the angle between subsequent estimates of the relevant W vector). Default is 1e-4.

Details

This fixed-point algorithm using approximations of negentropy is an implementation of the deflationary algorithm described in Hyvarinen, Karunen and Oja.2001. Independent Component Analysis

More specifically this is the implementation of iteration equation 8.43 on page 189 and the algorithm indicated in section 8.4.2 on page 194.

The three options for the non-quadratic function g are g1 = tanh(a1y); g2 = yexp(-y^2/2) and g3 = y^3. The last equates to the use of kurtosis.

The unmixing matrix is normalised throughout and the data is whitened in a pre-processing step.

Value

ws	A list per estimated independent component of unmixing matrices by iteration
W	The final converged estimate of the unmixing matrix
S	The estimated components by row with observations by column
thetas	A list per component of angles in radians between subsequent estimates of W
iters	The number of iterations to convergence for each component

Author(s)

Hans-Peter Bakker

References

See details

Examples

library(rmutil)
s <- rbind(runif(1000), rlaplace(1000), rnorm(100))
a <- matrix(runif(9), nrow = 3)
x <- a 

out <- negentropyICA(x)

par(mfrow = c(3,3))
plot(density(s[1,]), main = "signal 1")
plot(density(s[2,]), main = "signal 2")
plot(density(s[3,]), main = "signal 3")
plot(density(x[1,]), main = "mixture 1")
plot(density(x[2,]), main = "mixture 2")
plot(density(x[3,]), main = "mixture 3")
plot(density(out$S[1,]), main = "Indep. Comp. 1")
plot(density(out$S[2,]), main = "Indep. Comp. 2")
plot(density(out$S[3,]), main = "Indep. Comp. 3")

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function (x, g = g2, max.iters = 10, init.w = diag(m), tol = 1e-04) 
{
    source("my_Whiten.R")
    g1 <- expression(tanh(a1 * y))
    g2 <- expression(y * exp(-y^2/2))
    g3 <- expression(y^3)
    norm_vec <- function(x) {
        sqrt(sum(x^2))
    }
    negent <- function(wi, z, g, a1 = 1) {
        m <- nrow(z)
        y <- wi %*% z
        f <- function(y, g, a1) {
            eval(g[[1]])
        }
        gy <- f(y, g, a1)
        dg <- D(g, "y")
        as.vector((rowMeans(z * matrix(rep(gy, m), byrow = TRUE, 
            nrow = m)) - rowMeans(matrix(rep(eval(dg), m), byrow = TRUE, 
            nrow = m)) * t(wi)))
    }
    m <- nrow(x)
    n <- ncol(x)
    z <- my_Whiten(t(x))$z
    w <- init.w/norm_vec(init.w)
    iters <- max.iters
    ws <- vector("list", m)
    thetas <- vector("list", m)
    is <- vector()
    for (p in 1:m) {
        wp <- w[p, ]
        theta_vector <- vector()
        i <- 1
        repeat {
            y <- wp %*% z
            wp <- negent(wp, z, g2, a1 = 1)
            if (p > 1) {
                terms <- matrix(rep(0, (p - 1) * m), ncol = m)
                for (j in 1:(p - 1)) {
                  terms[j, ] <- t(wp %*% w[j, ]) %*% w[j, ]
                }
                wp <- wp - colSums(terms)
            }
            wp <- wp/norm_vec(wp)
            val <- wp %*% w[p, ]/(length(wp) * length(w[p, ]))
            theta <- acos(round(val, 6))
            theta_vector[i] <- theta
            w[p, ] <- wp
            is[p] <- i
            if (i == iters | abs(theta) < tol | abs(theta - pi) < 
                tol) {
                break
            }
            i <- i + 1
        }
        ws[[p]] <- w
        thetas[[p]] <- theta_vector
    }
    out <- list(ws = ws, W = ws[[m]], S = ws[[m]] %*% z, thetas = thetas, 
        iters = is)
  }

hanspeter6/icaPlay documentation built on May 2, 2020, 2:34 p.m.