Type:	Package
Title:	Blind Source Separation Methods Based on Joint Diagonalization and Some BSS Performance Criteria
Version:	2.0-4
Date:	2023-09-17
Author:	Klaus Nordhausen [aut, cre], Jean-Francois Cardoso [aut], Jari Miettinen [aut], Hannu Oja [aut], Esa Ollila [aut], Sara Taskinen [aut]
Maintainer:	Klaus Nordhausen <klausnordhausenR@gmail.com>
Imports:	clue, graphics
Suggests:	ICS, ICSNP
Description:	Cardoso's JADE algorithm as well as his functions for joint diagonalization are ported to R. Also several other blind source separation (BSS) methods, like AMUSE and SOBI, and some criteria for performance evaluation of BSS algorithms, are given. The package is described in Miettinen, Nordhausen and Taskinen (2017) <doi:10.18637/jss.v076.i02>.
License:	GPL-2 | GPL-3 [expanded from: GPL (≥ 2)]
LazyData:	true
NeedsCompilation:	yes
Packaged:	2023-09-17 14:20:22 UTC; admin
Repository:	CRAN
Date/Publication:	2023-09-17 23:02:33 UTC

Blind Source Separation Methods Based on Joint Diagonalization and Some BSS Performance Criteria

Description

Cardoso's JADE algorithm as well as his functions for joint diagonalization are ported to R. Also several other blind source separation (BSS) methods, like AMUSE and SOBI, and some criteria for performance evaluation of BSS algorithms, are given. The package is described in Miettinen, Nordhausen and Taskinen (2017) <doi:10.18637/jss.v076.i02>.

Details

Index of help topics:

AMUSE                   AMUSE Method for Blind Source Separation
CPPdata                 Cocktail Party Problem Data
ComonGAP                Comon's Gap
FG                      Joint Diagonalization of Real Positive-definite
                        Matrices
FOBI                    Function to perform FOBI for ICA
JADE                    JADE Algorithm for ICA
JADE-package            Blind Source Separation Methods Based on Joint
                        Diagonalization and Some BSS Performance
                        Criteria
MD                      Minimum Distance index MD
NSS.JD                  NSS.JD Method for Nonstationary Blind Source
                        Separation
NSS.SD                  NSS.SD Method for Nonstationary Blind Source
                        Separation
NSS.TD.JD               NSS.TD.JD Method for Nonstationary Blind Source
                        Separation
SIR                     Signal to Interference Ratio
SOBI                    SOBI Method for Blind Source Separation
amari.error             Amari Error
bss.components          Function to Extract Estimated Sources from an
                        Object of Class bss
cjd                     Joint Diagonalization of Complex Matrices
coef.bss                Coefficients of a bss Object
djd                     Function for Joint Diagonalization of k Square
                        Matrices in a Deflation Based Manner
k_JADE                  Fast Equivariant k-JADE Algorithm for ICA
multscatter             Function to Compute Several Scatter Matrices
                        for the Same Data
plot.bss                Plotting an Object of Class bss
print.bss               Printing an Object of Class bss
rjd                     Joint Diagonalization of Real Matrices

Author(s)

Klaus Nordhausen [aut, cre] (<https://orcid.org/0000-0002-3758-8501>), Jean-Francois Cardoso [aut], Jari Miettinen [aut] (<https://orcid.org/0000-0002-3270-7014>), Hannu Oja [aut] (<https://orcid.org/0000-0002-4945-5976>), Esa Ollila [aut], Sara Taskinen [aut] (<https://orcid.org/0000-0001-9470-7258>)

Maintainer: Klaus Nordhausen <klausnordhausenR@gmail.com>

References

Miettinen, J., Nordhausen, K. and Taskinen, S. (2017), Blind Source Separation Based on Joint Diagonalization in R: The Packages JADE and BSSasymp, Journal of Statistical Software, 76, 1–31, <doi:10.18637/jss.v076.i02>.

AMUSE Method for Blind Source Separation

Description

AMUSE method for the second order blind source separation problem. The function estimates the unmixing matrix in a second order stationary source separation model by jointly diagonalizing the covariance matrix and an autocovariance matrix at lag k.

Usage

AMUSE(x, ...)

## Default S3 method:
AMUSE(x, k = 1, ...)
## S3 method for class 'ts'
AMUSE(x, ...)

Arguments

Details

The lag k has a huge effect on the performance and it should be chosen so that the eigenvalues of autocovariance matrix are distinct. The function assumes always as many sources as there are time series.

Value

A list with class 'bss' containing the following components:

Author(s)

Klaus Nordhausen

References

Tong, L., Soon, V.C., Huang, Y.F. and Liu, R. (1990), AMUSE: a new blind identification algorithm, in Proceedings of IEEE International Symposium on Circuits and Systems 1990, 1784–1787.

Miettinen, J., Nordhausen, K., Oja, H. and Taskinen, S. (2012), Statistical properties of a blind source separation estimator for stationary time series, Statistics & Probability Letters, 82, 1865–1873.

Miettinen, J., Nordhausen, K. and Taskinen, S. (2017), Blind Source Separation Based on Joint Diagonalization in R: The Packages JADE and BSSasymp, Journal of Statistical Software, 76, 1–31, <doi:10.18637/jss.v076.i02>.

See Also

ts

Examples

# creating some toy data
A<- matrix(rnorm(9),3,3)
s1 <- arima.sim(list(ar=c(0.3,0.6)),1000)
s2 <- arima.sim(list(ma=c(-0.3,0.3)),1000)
s3 <- arima.sim(list(ar=c(-0.8,0.1)),1000)

S <- cbind(s1,s2,s3)
X <- S %*% t(A)

res1<-AMUSE(X)
res1
coef(res1)
plot(res1) # compare to plot.ts(S)
MD(coef(res1),A)

# input of a time series
X2<- ts(X, start=c(1961, 1), frequency=12)
plot(X2)
res2<-AMUSE(X2, k=2)
plot(res2)

Cocktail Party Problem Data

Description

This data set is a toy example for the so called cocktail party problem. In this case three sounds are mixed together with one noise source using four microphones.

Usage

data("CPPdata")

Format

A data frame with 50000 observations on the following 4 variables.

Mic1

the mixture recorded by the first microphone.

Mic2

the mixture recorded by the second microphone.

Mic3

the mixture recorded by the third microphone.

Mic4

the mixture recorded by the fourth microphone.

Details

The three original source files were kindly provided by Ella Bingham and are also available online at the following locations: https://research.ics.aalto.fi/ica/cocktail/source5.wav, https://research.ics.aalto.fi/ica/cocktail/source7.wav and https://research.ics.aalto.fi/ica/cocktail/source9.wav.

Note that the original sound files are included in the package's subfolder datafiles. In the example section we illustrate how the CPPdata was created. An example analysis of the data is given in Miettinen et al. (2017).

Source

Ella Bingham

References

Miettinen, J., Nordhausen, K. and Taskinen, S. (2017), Blind Source Separation Based on Joint Diagonalization in R: The Packages JADE and BSSasymp, Journal of Statistical Software, 76, 1–31, <doi:10.18637/jss.v076.i02>.

Examples

## Not run: 
# the data was created as follows:
library("tuneR")
S1 <- readWave(system.file("datafiles/source5.wav", package = "JADE"))
S2 <- readWave(system.file("datafiles/source7.wav", package = "JADE"))
S3 <- readWave(system.file("datafiles/source9.wav", package = "JADE"))

set.seed(321)
NOISE <- noise("white", duration = 50000)
S <- cbind(S1@left, S2@left, S3@left, NOISE@left)
S <- scale(S, center = FALSE, scale = apply(S, 2, sd))
St <- ts(S, start = 0, frequency = 8000)
p <- 4
A <- matrix(runif(p^2, 0, 1), p, p)
A

X <- tcrossprod(St, A)
Xt <- as.ts(X)

colnames(X) <- c("Mic1", "Mic2", "Mic3", "Mic4")
CPPdata <- as.data.frame(X)

## End(Not run)

Comon's Gap

Description

Comon's GAP criterion to evaluate the performance of an ICA algorithm.

Usage

ComonGAP(A, A.hat)

Arguments

Details

Comon's GAP criterion is permutation and scale invariant. It can take every positive value and 0 corresponds to an optimal separation. If A is however nearly singular the values of the criterion can be huge.

Note that this function assumes the ICA model is X = S A', as is assumed by JADE and ics. However fastICA and PearsonICA assume X = S A. Therefore matrices from those functions have to be transposed first.

Value

The value of the Comon's GAP.

Author(s)

Klaus Nordhausen

References

Comon, P., (1994), Independent Component Analysis, A new concept?, Signal Processing, 36, 287–314.

See Also

amari.error, SIR

Examples

S <- cbind(rt(1000, 4), rnorm(1000), runif(1000))
A <- matrix(rnorm(9), ncol = 3)
X <- S %*% t(A)

A.hat <- JADE(X, 3)$A
ComonGAP(A, A.hat)

Joint Diagonalization of Real Positive-definite Matrices

Description

This is a slightly modified version of Flury's FG algorithm for the joint diagonalization of k positive-definite matrices. The underlying function is written in C.

Usage

FG(X, weight = NULL, init = NULL, maxiter = 100, eps = 1e-06, na.action = na.fail)

Arguments

Value

A list with the components

Author(s)

Jari Miettinen

References

Flury, B. D. (1998), Common principal components and related models, Wiley, New York.

See Also

rjd, rjd.fortran

Function to perform FOBI for ICA

Description

The FOBI method for independent component analysis (ICA). We assume that all components have different kurtosis values.

Usage

FOBI(X, na.action = na.fail)

Arguments

Value

A list with class 'bss' containing the following components:

Note

More general is the function ics in the ICS package.

Author(s)

Klaus Nordhausen

References

Cardoso, J.-F. (1989), Source separation using higher order moments, in Proceedings of IEEE International Conference on Accoustics, Speech and Signal Processing, 2109–2112.

Miettinen, J., Taskinen S., Nordhausen, K. and Oja, H. (2015), Fourth Moments and Independent Component Analysis, Statistical Science, 30, 372–390.

Miettinen, J., Nordhausen, K. and Taskinen, S. (2017), Blind Source Separation Based on Joint Diagonalization in R: The Packages JADE and BSSasymp, Journal of Statistical Software, 76, 1–31, <doi:10.18637/jss.v076.i02>.

See Also

ics

Examples

# 3 source and 3 signals

S <- cbind(rt(1000, 4), rnorm(1000), runif(1000))
A <- matrix(rnorm(9), ncol = 3)
X <- S %*% t(A)
res<-FOBI(X)
MD(coef(res),A)

JADE Algorithm for ICA

Description

This is an R version of Cardoso's JADE ICA algorithm (for real data) ported from matlab. The ported version is 1.5. Some minor changes compared to the matlab function are explained in the details section. The matlab code can be found for example on the ICA central homepage.

The function uses frjd for the joint diagonalization.

Usage

JADE(X, n.comp = NULL, eps = 1e-06, maxiter = 100, na.action = na.fail)

Arguments

Details

Some minor modifications were done when porting the function to R, and they are:

1

the model assumed here is X=S A' +\mu. Therefore S and X have one row per observation. Note that this still differs from the model definition in R of FastICA and PearsonICA but agrees with ics.

2

The whitening covariance matrix is divided by n-1 and not n (n = number of observations).

3

The initial value for the joint diagonalisation is always I.

4

The original eps would be \frac{1}{100\sqrt{n}}.

It is also worth mentioning that the estimated independent components S are scaled to unit variance and are ordered in such a way, that their fourth moments are in the decreasing order. The signs of the unmixing matrix W are fixed so that the sum of the elements on each row is positive.

The code is based on the original matlab code ("MatlabjadeR.m").

Value

A list with class 'bss' containing the following components:

Author(s)

Jean-Francois Cardoso. Ported to R by Klaus Nordhausen

References

Cardoso, J.-F. and Souloumiac, A., (1993), Blind beamforming for non Gaussian signals, IEE Proceedings-F, 140, 362–370.

Miettinen, J., Taskinen S., Nordhausen, K. and Oja, H. (2015), Fourth Moments and Independent Component Analysis, Statistical Science, 30, 372–390.

Miettinen, J., Nordhausen, K. and Taskinen, S. (2017), Blind Source Separation Based on Joint Diagonalization in R: The Packages JADE and BSSasymp, Journal of Statistical Software, 76, 1–31, <doi:10.18637/jss.v076.i02>.

Examples

# 3 source and 3 signals

S <- cbind(rt(1000, 4), rnorm(1000), runif(1000))
A <- matrix(rnorm(9), ncol = 3)
X <- S %*% t(A)
res<-JADE(X,3)
res$A
res$W
res$S[1:10,]
(sweep(X,2,res$Xmu) %*% t(res$W))[1:10,]
round(res$W %*% A,4) 

# 2 sources and 3 signals

S2 <- cbind(rt(1000, 4), rnorm(1000)) 
A2 <- matrix(rnorm(6), ncol = 2)
X2 <- S2 %*% t(A2)
res2 <-JADE(X2,2)
res2$A
res2$W
res2$S[1:10,]
(sweep(X2,2,res2$Xmu) %*% t(res2$W))[1:10,]
SIR(S2,res2$S)

Minimum Distance index MD

Description

Computes the Minimum Distance index MD to evaluate the performance of an ICA algorithm.

Usage

MD(W.hat, A)

Arguments

Details

MD(\hat{W},A)=\frac{1}{\sqrt{p-1}} \inf_{P D}{||PD \hat{W} A-I||,}

where P is a permutation matrix and D a diagonal matrix with nonzero diagonal entries.

The step that minimizes the index of the set over all permutation matrix can be expressed as a linear sum assignment problem (LSAP) for which we use as solver the Hungarian method implemented as solve_LASP in the clue package.

Note that this function assumes the ICA model is X = S A', as is assumed by JADE and ics. However fastICA and PearsonICA assume X = S A. Therefore matrices from those functions have to be transposed first.

The MD index is scaled in such a way, that it takes a value between 0 and 1. And 0 corresponds to an optimal separation.

Value

The value of the MD index.

Author(s)

Klaus Nordhausen

References

Ilmonen, P., Nordhausen, K., Oja, H. and Ollila, E. (2010), A New Performance Index for ICA: Properties, Computation and Asymptotic Analysis. In Vigneron, V., Zarzoso, V., Moreau, E., Gribonval, R. and Vincent, E. (editors) Latent Variable Analysis and Signal Separation, 229–236, Springer.

Nordhausen, K., Ollila, E. and Oja, H. (2011), On the Performance Indices of ICA and Blind Source Separation. In the Proceedings of 2011 IEEE 12th International Workshop on Signal Processing Advances in Wireless Communications (SPAWC 2011), 486–490.

Miettinen, J., Nordhausen, K. and Taskinen, S. (2017), Blind Source Separation Based on Joint Diagonalization in R: The Packages JADE and BSSasymp, Journal of Statistical Software, 76, 1–31, <doi:10.18637/jss.v076.i02>.

See Also

ComonGAP, SIR, amari.error, solve_LSAP

Examples

S <- cbind(rt(1000, 4), rnorm(1000), runif(1000))
A <- matrix(rnorm(9), ncol = 3)
X <- S %*% t(A)

W.hat <- JADE(X, 3)$W
MD(W.hat, A)

NSS.JD Method for Nonstationary Blind Source Separation

Description

The NSS.JD method for nonstationary blind source separation. The method first whitens the complete data and then divides it into K time intervals. Then frjd is used to jointly diagonalize the covariance matrices computed for the individual time intervals to find the sources.

Usage

NSS.JD(X, ...)

## Default S3 method:
NSS.JD(X,  K=12, Tau=0, n.cuts=NULL, eps = 1e-06, maxiter = 100, ...)
## S3 method for class 'ts'
NSS.JD(X, ...)

Arguments

Details

The model assumes that the mean of the p-variate time series is constant but the variances change over time.

Value

A list with class 'bss' containing the following components:

Author(s)

Klaus Nordhausen

References

Choi S. and Cichocki A. (2000), Blind separation of nonstationary sources in noisy mixtures, Electronics Letters, 36, 848–849.

Choi S. and Cichocki A. (2000), Blind separation of nonstationary and temporally correlated sources from noisy mixtures, Proceedings of the 2000 IEEE Signal Processing Society Workshop Neural Networks for Signal Processing X, 1, 405–414.

Nordhausen K. (2014), On robustifying some second order blind source separation methods for nonstationary time series, Statistical Papers, 55, 141–156.

Miettinen, J., Nordhausen, K. and Taskinen, S. (2017), Blind Source Separation Based on Joint Diagonalization in R: The Packages JADE and BSSasymp, Journal of Statistical Software, 76, 1–31, <doi:10.18637/jss.v076.i02>.

See Also

ts, NSS.SD, NSS.TD.JD

Examples

n <- 1000
s1 <- rnorm(n)
s2 <- 2*sin(pi/200*1:n)* rnorm(n)
s3 <- c(rnorm(n/2), rnorm(100,0,2), rnorm(n/2-100,0,1.5))
S <- cbind(s1,s2,s3)
plot.ts(S)
A<-matrix(rnorm(9),3,3)
X<- S%*%t(A)

NSS2 <- NSS.JD(X)
NSS2
MD(coef(NSS2),A)
plot(NSS2)
cor(NSS2$S,S)

NSS2b <- NSS.JD(X, Tau=1)
MD(coef(NSS2b),A)

NSS2c <- NSS.JD(X, n.cuts=c(1,300,500,600,1000))
MD(coef(NSS2c),A)

NSS.SD Method for Nonstationary Blind Source Separation

Description

The NSS.SD method for nonstationary blind source separation. The function estimates the unmixing matrix in a nonstationary source separation model by simultaneously diagonalizing two covariance matrices computed for different time intervals.

Usage

NSS.SD(X, ...)

## Default S3 method:
NSS.SD(X, n.cut=NULL, ...)
## S3 method for class 'ts'
NSS.SD(X, ...)

Arguments

Details

The model assumes that the mean of the p-variate time series is constant but the variances change over time.

Value

A list with class 'bss' containing the following components:

Author(s)

Klaus Nordhausen

References

Choi S. and Cichocki A. (2000), Blind separation of nonstationary sources in noisy mixtures, Electronics Letters, 36, 848–849.

Choi S. and Cichocki A. (2000), Blind separation of nonstationary and temporally correlated sources from noisy mixtures, Proceedings of the 2000 IEEE Signal Processing Society Workshop Neural Networks for Signal Processing X, 1, 405–414.

Nordhausen K. (2014), On robustifying some second order blind source separation methods for nonstationary time series, Statistical Papers, 55, 141–156.

Miettinen, J., Nordhausen, K. and Taskinen, S. (2017), Blind Source Separation Based on Joint Diagonalization in R: The Packages JADE and BSSasymp, Journal of Statistical Software, 76, 1–31, <doi:10.18637/jss.v076.i02>.

See Also

ts, NSS.JD, NSS.TD.JD

Examples

n <- 1000
s1 <- rnorm(n)
s2 <- 2*sin(pi/200*1:n)* rnorm(n)
s3 <- c(rnorm(n/2), rnorm(100,0,2), rnorm(n/2-100,0,1.5))
S <- cbind(s1,s2,s3)
plot.ts(S)
A<-matrix(rnorm(9),3,3)
X<- S%*%t(A)

NSS1 <- NSS.SD(X)
NSS1
MD(coef(NSS1),A)
plot(NSS1)
cor(NSS1$S,S)

NSS1b <- NSS.SD(X, n.cut=400)
MD(coef(NSS1b),A)

NSS1c <- NSS.SD(X, n.cut=c(1,600,1000))
MD(coef(NSS1c),A)

NSS.TD.JD Method for Nonstationary Blind Source Separation

Description

The NSS.TD.JD method for nonstationary blind source separation. The method first whitens the complete data and then divides it into K time intervals. It is then assumed that within each interval the time series is approximately second order stationary and within each interval L autocovariance are computed. The underlying sources are then found by jointly diagonalizing the K*L autocovariance matrices using frjd.

Usage

NSS.TD.JD(X, ...)

## Default S3 method:
NSS.TD.JD(X, K=12, Tau=0:11, n.cuts=NULL, eps = 1e-06, maxiter = 100, ...)
## S3 method for class 'ts'
NSS.TD.JD(X, ...)

Arguments

Details

The model assumes that the mean of the p-variate time series is constant but the variances change over time.

Value

A list with class 'bss' containing the following components:

Author(s)

Klaus Nordhausen

References

Choi S. and Cichocki A. (2000), Blind separation of nonstationary sources in noisy mixtures, Electronics Letters, 36, 848–849.

Choi S. and Cichocki A. (2000), Blind separation of nonstationary and temporally correlated sources from noisy mixtures, Proceedings of the 2000 IEEE Signal Processing Society Workshop Neural Networks for Signal Processing X, 1, 405–414.

Nordhausen K. (2014), On robustifying some second order blind source separation methods for nonstationary time series, Statistical Papers, 55, 141–156.

Miettinen, J., Nordhausen, K. and Taskinen, S. (2017), Blind Source Separation Based on Joint Diagonalization in R: The Packages JADE and BSSasymp, Journal of Statistical Software, 76, 1–31, <doi:10.18637/jss.v076.i02>.

See Also

ts, NSS.JD, NSS.JD, SOBI

Examples

n <- 1000
s1 <- rnorm(n)
s2 <- 2*sin(pi/200*1:n)* rnorm(n)
s3 <- c(rnorm(n/2), rnorm(100,0,2), rnorm(n/2-100,0,1.5))
S <- cbind(s1,s2,s3)
plot.ts(S)
A<-matrix(rnorm(9),3,3)
X<- S%*%t(A)

NSS3 <- NSS.TD.JD(X)
NSS3
MD(coef(NSS3),A)
plot(NSS3)
cor(NSS3$S,S)

NSS3b <- NSS.TD.JD(X, Tau=c(0,3,7,12), K=6)
MD(coef(NSS3b),A)

NSS3c <- NSS.TD.JD(X, n.cuts=c(1,300,500,600,1000))
MD(coef(NSS3c),A)

Signal to Interference Ratio

Description

Computes the signal to interference ratio between true and estimated signals

Usage

SIR(S, S.hat)

Arguments

Details

The signal to interference ratio is measured in dB and values over 20 are thought to be good. It is scale and permutation invariant and can be seen as measuring the correlation between the matched true and estimated signals.

Value

The value of the signal to interference ratio.

Author(s)

Klaus Nordhausen

References

Eriksson, J., Karvanen, J. and Koivunen, V. (2000), Source distribution adaptive maximum likelihood estimation in ICA model, Proceedings of the second international workshop on independent component analysis and blind source separation (ICA 2000), 227–232.

See Also

amari.error, ComonGAP

Examples

S <- cbind(rt(1000, 4), rnorm(1000), runif(1000))
A <- matrix(rnorm(9), ncol = 3)
X <- S %*% t(A)

S.hat <- JADE(X, 3)$S
SIR(S, S.hat)

SOBI Method for Blind Source Separation

Description

The SOBI method for the second order blind source separation problem. The function estimates the unmixing matrix in a second order stationary source separation model by jointly diagonalizing the covariance matrix and several autocovariance matrices at different lags.

Usage

SOBI(X, ...)

## Default S3 method:
SOBI(X, k=12, method="frjd", eps = 1e-06, maxiter = 100, ...)
## S3 method for class 'ts'
SOBI(X, ...)

Arguments

.

Details

The order of the estimated components is fixed so that the sums of squared autocovariances are in the decreasing order.

Value

A list with class 'bss' containing the following components:

Author(s)

Klaus Nordhausen

References

Belouchrani, A., Abed-Meriam, K., Cardoso, J.F. and Moulines, R. (1997), A blind source separation technique using second-order statistics, IEEE Transactions on Signal Processing, 434–444.

Miettinen, J., Nordhausen, K., Oja, H. and Taskinen, S. (2014), Deflation-based Separation of Uncorrelated Stationary Time Series, Journal of Multivariate Analysis, 123, 214–227.

Miettinen, J., Illner, K., Nordhausen, K., Oja, H., Taskinen, S. and Theis, F.J. (2016), Separation of Uncorrelated Stationary Time Series Using Autocovariance Matrices, Journal of Time Series Analysis, 37, 337–354.

Miettinen, J., Nordhausen, K. and Taskinen, S. (2017), Blind Source Separation Based on Joint Diagonalization in R: The Packages JADE and BSSasymp, Journal of Statistical Software, 76, 1–31, <doi:10.18637/jss.v076.i02>.

See Also

ts

Examples

# creating some toy data
A<- matrix(rnorm(9),3,3)
s1 <- arima.sim(list(ar=c(0.3,0.6)),1000)
s2 <- arima.sim(list(ma=c(-0.3,0.3)),1000)
s3 <- arima.sim(list(ar=c(-0.8,0.1)),1000)

S <- cbind(s1,s2,s3)
X <- S %*% t(A)

res1<-SOBI(X)
res1
coef(res1)
plot(res1) # compare to plot.ts(S)
MD(coef(res1),A)

# input of a time series
X2<- ts(X, start=c(1961, 1), frequency=12)
plot(X2)
res2<-SOBI(X2, k=c(5,10,1,4,2,9,10))
plot(res2)

Amari Error

Description

Computes the Amari Error to evaluate the performance of an ICA algorithm.

Usage

amari.error(W.hat, A, standardize = F)

Arguments

Details

The Amari Error can be used in simulation studies to evaluate the performance of an ICA algorithm. The Amari error is permutation invariant but not scale invariant. Therefore if different algorithms should be compared the matrices should be scaled in the same way. If standardize is TRUE, this will be done by the function by standardizing 'W.hat' and the inverse of 'A' in such a way, that every row has length 1, the largest absolute value of the row has a positive sign and the rows are ordered decreasingly according to their largest values.

Note that this function assumes the ICA model is X = S A', as is assumed by JADE and ics. However fastICA and PearsonICA assume X = S A. Therefore matrices from those functions have to be transposed first.

The Amari Error is scaled in such a way, that it takes a value between 0 and 1. And 0 corresponds to an optimal separation.

Value

The value of the Amari Error.

Author(s)

Klaus Nordhausen

References

Amari, S., Cichocki, A. and Yang, H.H. (1996), A new learning algorithm for blind signal separation, Advances in Neural Information Processing Systems, 8, 757–763.

Nordhausen, K., Ollila, E. and Oja, H. (2011), On the Performance Indices of ICA and Blind Source Separation. In the Proceedings of 2011 IEEE 12th International Workshop on Signal Processing Advances in Wireless Communications (SPAWC 2011), 486–490.

See Also

ComonGAP, SIR

Examples

S <- cbind(rt(1000, 4), rnorm(1000), runif(1000))
A <- matrix(rnorm(9), ncol = 3)
X <- S %*% t(A)

W.hat <- JADE(X, 3)$W
amari.error(W.hat, A)
amari.error(W.hat, A, TRUE)

Function to Extract Estimated Sources from an Object of Class bss

Description

Extracts the sources estimated by an bss method.

Usage

bss.components(object)

Arguments

Author(s)

Klaus Nordhausen

Examples

A<- matrix(rnorm(9),3,3)
s1 <- arima.sim(list(ar=c(0.3,0.6)),1000)
s2 <- arima.sim(list(ma=c(-0.3,0.3)),1000)
s3 <- arima.sim(list(ar=c(-0.8,0.1)),1000)

S <- cbind(s1,s2,s3)
X <- S %*% t(A)

res1<-AMUSE(X)
head(bss.components(res1))
colMeans(bss.components(res1))
cov(bss.components(res1))

Joint Diagonalization of Complex Matrices

Description

This is an R version of Cardoso's joint_diag matlab function for joint diagonalization of k complex-valued square matrices.

Usage

cjd(X, eps = 1e-06, maxiter = 100)

Arguments

Value

Author(s)

Jean-Francois Cardoso. Ported to R by Klaus Nordhausen.

References

Cardoso, J.-F. and Souloumiac, A., (1996), Jacobi angles for simultaneous diagonalization, SIAM J. Mat. Anal. Appl., 17, 161–164.

See Also

rjd, rjd.fortran

Examples

D1 <- diag(complex(real=runif(3,0,2), imaginary=runif(3)))
D2 <- diag(complex(real=runif(3,0,2), imaginary=runif(3)))
D3 <- diag(complex(real=runif(3,0,2), imaginary=runif(3)))
D4 <- diag(complex(real=runif(3,0,2), imaginary=runif(3)))

Z <- matrix(runif(9), ncol = 3)
V <- eigen(Z %*% t(Z))$vectors

M1 <- t(V)%*%D1%*%V
M2 <- t(V)%*%D2%*%V
M3 <- t(V)%*%D3%*%V
M4 <- t(V)%*%D4%*%V
MS <- rbind(M1,M2,M3,M4)
Ms <- array(0,dim=c(3,3,4))
Ms[,,1]<-M1
Ms[,,3]<-M3
Ms[,,2]<-M2
Ms[,,4]<-M4
res.array <- cjd(Ms)
res.mat <- cjd(MS)
Re(res.array$V)
V
round(V%*%Re(res.array$V),2)
round(V%*%Re(res.mat$V),2)

Coefficients of a bss Object

Description

Extracts the estimated unmixing matrix from an object of class bss.

Usage

## S3 method for class 'bss'
coef(object, ...)

Arguments

Author(s)

Klaus Nordhausen

Examples

A<- matrix(rnorm(9),3,3)
s1 <- arima.sim(list(ar=c(0.3,0.6)),1000)
s2 <- arima.sim(list(ma=c(-0.3,0.3)),1000)
s3 <- arima.sim(list(ar=c(-0.8,0.1)),1000)

S <- cbind(s1,s2,s3)
X <- S %*% t(A)

res1<-AMUSE(X)
coef(res1)
coef(res1) %*% A # should be a matrix with one dominant element in each row and column

Function for Joint Diagonalization of k Square Matrices in a Deflation Based Manner

Description

This function jointly diagonalizes k real-valued square matrices by searching an orthogonal matrix in a deflation based manner.

Usage

djd(X, G = "max", r = 2, eps = 1e-06, maxiter = 500)

Arguments

Details

Denote the square matrices as A_i, i=1,\ldots,k. This algorithm searches then an orthogonal matrix W so that D_i=W'A_iW is diagonal for all i. If the A_i commute then there is an exact solution. If not, the function will perform an approximate joint diagonalization by maximizing \sum G(w_j' A_i w_j) where w_j are the orthogonal vectors in W.

The function G can be choosen to be of the form G(x) = |x|^r or G(x) = log(x). If G="max" is chosen, the function G is of the form G(x) = |x|^r, and the diagonalization criterion will be maximized globally at each stage by choosing an appropriate initial value from a set of random vectors. If G="pow" or G="log" are chosen, the initial values are the eigenvectors of A_1 which plays hence a special role.

Value

The matrix W

Author(s)

Klaus Nordhausen, Jari Miettinen

References

Nordhausen, K., Gutch, H. W., Oja, H. and Theis, F.J. (2012): Joint Diagonalization of Several Scatter Matrices for ICA, in LVA/ICA 2012, LNCS 7191, pp. 172–179.

Miettinen, J., Nordhausen, K., Oja, H. and Taskinen, S. (2014), Deflation-based Separation of Uncorrelated Stationary Time Series, Journal of Multivariate Analysis, 123, 214–227.

Miettinen, J., Nordhausen, K. and Taskinen, S. (2017), Blind Source Separation Based on Joint Diagonalization in R: The Packages JADE and BSSasymp, Journal of Statistical Software, 76, 1–31, <doi:10.18637/jss.v076.i02>.

Examples

Z <- matrix(runif(9), ncol = 3)
U <- eigen(Z %*% t(Z))$vectors
D1 <- diag(runif(3))
D2 <- diag(runif(3))
D3 <- diag(runif(3))
D4 <- diag(runif(3))

X.matrix <- array(0, dim=c(3, 3, 4))
X.matrix[,,1] <- t(U) %*% D1 %*% U
X.matrix[,,2] <- t(U) %*% D2 %*% U
X.matrix[,,3] <- t(U) %*% D3 %*% U
X.matrix[,,4] <- t(U) %*% D4 %*% U

W1 <- djd(X.matrix)
round(U %*% W1, 4) # should be a signed permutation 
                     # matrix if W1 is correct.

W2 <- djd(X.matrix, r=1)
round(U %*% W2, 4) # should be a signed permutation 
                     # matrix if W2 is correct.

W3 <- djd(X.matrix, G="l")
round(U %*% W3, 4) # should be a signed permutation 
                     # matrix if W3 is correct.

Fast Equivariant k-JADE Algorithm for ICA

Description

This algorithm generalizes the JADE algorithm, as it provides JADE when k is set to the number of dimensions. Otherwise k can be considered as a way to reduce the number of cumulant matrices to be jointly diagonalized. Hence small values of k speed up the method considerably in high-dimensional cases. In general, k can be considered as maximum number of underlying identical sources.

The function uses FOBI as an initial estimate and frjd for the joint diagonalization.

Usage

k_JADE(X, k = 1, eps = 1e-06, maxiter = 100, na.action = na.fail)

Arguments

Details

The order of the estimated components is fixed so that their fourth moments are in the decreasing order.

Value

A list with class 'bss' containing the following components:

Note

The function uses FOBI as initial estimate and frjd for the joint diagonalization.

Author(s)

Jari Miettinen

References

Miettinen, J., Nordhausen, K., Oja, H. and Taskinen, S. (2013), Fast Equivariant JADE, In the Proceedings of 38th IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP 2013), 6153–6157.

Miettinen, J., Nordhausen, K. and Taskinen, S. (2017), Blind Source Separation Based on Joint Diagonalization in R: The Packages JADE and BSSasymp, Journal of Statistical Software, 76, 1–31, <doi:10.18637/jss.v076.i02>.

See Also

JADE, FOBI, frjd

Examples

# 3 source and 3 signals

S <- cbind(rt(1000, 4), rnorm(1000), runif(1000))
A <- matrix(rnorm(9), ncol = 3)
X <- S %*% t(A)
res_k1<-k_JADE(X,1)
res_k1$A
res_k1$W
res_k1$S[1:10,]

MD(coef(res_k1),A) 

Function to Compute Several Scatter Matrices for the Same Data

Description

The function can be used to compute several scatter matrices for the same data.

Usage

multscatter(scatterlist, X, toshape = TRUE)

Arguments

Details

It is important that the functions do not need any additional imput and that they return only the p times p scatter matrix. Hence it might be sometimes necessary to write wrappers for some of the functions. See examples.

Value

An array of dimension c(p,p,k) where k is the number of scatter matrices.

Author(s)

Klaus Nordhausen

Examples

# example requires the packages ICS and ICSNP
library(ICSNP)
X <- cbind(rexp(1000), rt(1000,6), runif(1000)) 

my.tM1 <- function(X,df=1) tM(X,)$V
my.tM2 <- function(X,df=2) tM(X,)$V

multscatter(c("cov","cov4","HP1.shape","my.tM1", "my.tM2"), X)
multscatter(c("cov","cov4","HP1.shape","my.tM1", "my.tM2"), X, toshape=FALSE)

Plotting an Object of Class bss

Description

Plots the estimated sources resulting from an bss method. If the bss method is based on second order assumptions and returned the sources as a time series object it will plot the sources using plot.ts, otherwise it will plot a scatter plot matrix using pairs or plot if there are only two sources.

Usage

## S3 method for class 'bss'
plot(x, ...)

Arguments

Author(s)

Klaus Nordhausen

See Also

plot.ts, pairs, plot

Examples

A<- matrix(rnorm(9),3,3)
s1 <- arima.sim(list(ar=c(0.3,0.6)),1000)
s2 <- arima.sim(list(ma=c(-0.3,0.3)),1000)
s3 <- arima.sim(list(ar=c(-0.8,0.1)),1000)

S <- cbind(s1,s2,s3)
X <- S %*% t(A)

res1 <- AMUSE(X)
plot(res1)
# not so useful:
plot(res1, plot.type = "single", col=1:3)

# not meaningful for this data
res2 <- JADE(X)
plot(res2)

Printing an Object of Class bss

Description

Prints an object of class bss. It prints all elements of the list of class bss except the component S which is the source matrix.

Usage

## S3 method for class 'bss'
print(x, ...)

Arguments

Author(s)

Klaus Nordhausen

Joint Diagonalization of Real Matrices

Description

This is an R version of Cardoso's rjd matlab function for joint diagonalization of k real-valued square matrices. A version written in C is also available and preferrable.

Usage

rjd(X, eps = 1e-06, maxiter = 100, na.action = na.fail)
frjd(X, weight = NULL, maxiter = 100, eps = 1e-06, na.action = na.fail)
frjd.int(X, maxiter = 100, eps = 1e-06)
rjd.fortran(X, weight = NULL, maxiter = 100, eps = 1e-06, na.action = na.fail)

Arguments

Details

Denote the square matrices as A_i, i=1,\ldots,k. The algorithm searches an orthogonal matrix V so that D_i=V'A_iV is diagonal for all i. If the A_i commute then there is an exact solution. Otherwise, the function will perform an approximate joint diagonalization by trying to make the D_i as diagonal as possible.

Cardoso points out that notion of approximate joint diagonalization is ad hoc and very small values of eps make in that case not much sense since the diagonality criterion is ad hoc itself.

rjd, frjd and rjd.fortran terminate with an error in case maxiter is reach without convergence whereas frjd_int returns the current state at when maxiter is reached and does not warn about convergence problems.

Value

A list with the components

Author(s)

Jean-Francois Cardoso. Ported to R by Klaus Nordhausen. C code by Jari Miettinen

References

Cardoso, J.-F. and Souloumiac, A., (1996), Jacobi angles for simultaneous diagonalization, SIAM J. Mat. Anal. Appl., 17, 161–164.

Miettinen, J., Nordhausen, K. and Taskinen, S. (2017), Blind Source Separation Based on Joint Diagonalization in R: The Packages JADE and BSSasymp, Journal of Statistical Software, 76, 1–31, <doi:10.18637/jss.v076.i02>.

Examples

Z <- matrix(runif(9), ncol = 3)
U <- eigen(Z %*% t(Z))$vectors
D1 <- diag(runif(3))
D2 <- diag(runif(3))
D3 <- diag(runif(3))
D4 <- diag(runif(3))

X.matrix <- rbind(t(U) %*% D1 %*% U, t(U) %*% D2 %*% U,
                  t(U) %*% D3 %*% U, t(U) %*% D4 %*% U)
res.matrix <- rjd(X.matrix)
res.matrix$V
round(U %*% res.matrix$V, 4) # should be a signed permutation 
                             # matrix if V is correct.

round(res.matrix$D, 4)

# compare to C version

#res.matrix.C <- frjd(X.matrix)
#res.matrix.C$V
#round(U %*% res.matrix.C$V, 4)
#round(res.matrix.C$D, 4)

X.array <- aperm(array(t(X.matrix), dim = c(3,3,4)), c(2,1,3))

res.array <- rjd(X.array)
round(res.array$D, 4)

res.array.C <- frjd(X.array)
round(res.array.C$D, 4)

res.array.C2 <- frjd.int(X.array)
round(res.array.C2$D, 4)