, Jean-Francois
Cardoso [aut], Jari Miettinen [aut]
, Hannu Oja [aut]
, Esa Ollila [aut],
Sara Taskinen [aut] | Title: | 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>. |
| Authors: | Klaus Nordhausen [aut, cre]
, Jean-Francois
Cardoso [aut], Jari Miettinen [aut]
, Hannu Oja [aut]
, Esa Ollila [aut],
Sara Taskinen [aut] |
| Maintainer: | Klaus Nordhausen <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 2.0-4 |
| Built: | 2024-08-14 06:07:47 UTC |
| Source: | CRAN |
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>.
| Package: | JADE |
| Type: | Package |
| Title: | Blind Source Separation Methods Based on Joint Diagonalization and Some BSS Performance Criteria |
| Version: | 2.0-4 |
| Date: | 2023-09-17 |
| Authors@R: | c(person(given = "Klaus", family = "Nordhausen", role = c("aut", "cre"), email = "[email protected]", comment = c(ORCID = "0000-0002-3758-8501")), person(given = "Jean-Francois", family = "Cardoso", role = "aut"), person(given = "Jari", family = "Miettinen", role = "aut", comment = c(ORCID = "0000-0002-3270-7014")), person(given = "Hannu", family = "Oja", role = "aut", comment = c(ORCID = "0000-0002-4945-5976")), person(given = "Esa", family = "Ollila", role = "aut"), person(given = "Sara", family = "Taskinen", role = "aut", comment = c(ORCID = "0000-0001-9470-7258"))) |
| Author: | 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 <[email protected]> |
| 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) |
| LazyData: | true |
| NeedsCompilation: | yes |
| Packaged: | 2023-09-17 14:20:22 UTC; admin |
| Repository: | CRAN |
| Date/Publication: | 2023-09-17 23:02:33 UTC |
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
Further information is available in the following vignettes:
JADE-BSSasymp |
Blind Source Separation Based on Joint Diagonalization in R: The Packages JADE and BSSasymp (source, pdf) |
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 <[email protected]>
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>.
Computes the Amari Error to evaluate the performance of an ICA algorithm.
amari.error(W.hat, A, standardize = F)amari.error(W.hat, A, standardize = F)
W.hat |
The estimated square unmixing matrix W. |
A |
The true square mixing matrix A. |
standardize |
Logical value if A and W.hat need to be standardized. Default is FALSE. |
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 , as is assumed by JADE and ics. However fastICA and
PearsonICA assume . 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.
The value of the Amari Error.
Klaus Nordhausen
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.
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)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)
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.
AMUSE(x, ...) ## Default S3 method: AMUSE(x, k = 1, ...) ## S3 method for class 'ts' AMUSE(x, ...)AMUSE(x, ...) ## Default S3 method: AMUSE(x, k = 1, ...) ## S3 method for class 'ts' AMUSE(x, ...)
x |
a numeric matrix or a multivariate time series object of class |
k |
integer lag for the autocovariance matrix, must be larger than 0. Default is 1. |
... |
further arguments to be passed to or from methods. |
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.
A list with class 'bss' containing the following components:
W |
estimated unmixing matrix. |
EV |
eigenvectors of autocovariance matrix. |
k |
lag of the autocovariance matrix used. |
S |
estimated sources as time series objected standardized to have mean 0 and unit variances. |
Klaus Nordhausen
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>.
# 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)# 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)
Extracts the sources estimated by an bss method.
bss.components(object)bss.components(object)
object |
object of class bss |
Klaus Nordhausen
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))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))
This is an R version of Cardoso's joint_diag matlab function for joint diagonalization of k complex-valued square matrices.
cjd(X, eps = 1e-06, maxiter = 100)cjd(X, eps = 1e-06, maxiter = 100)
X |
A matrix of k stacked pxp complex matrices with dimension c(kp,p) or an array with dimension c(p,p,k). |
eps |
Convergence tolerance. |
maxiter |
Maximum number of iterations. |
V |
An orthogonal matrix. |
D |
A stacked matrix with the diagonal matrices or an array with the diagonal matrices. The form of the output depends on the form of the input. |
Jean-Francois Cardoso. Ported to R by Klaus Nordhausen.
Cardoso, J.-F. and Souloumiac, A., (1996), Jacobi angles for simultaneous diagonalization, SIAM J. Mat. Anal. Appl., 17, 161–164.
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)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)
Extracts the estimated unmixing matrix from an object of class bss.
## S3 method for class 'bss' coef(object, ...)## S3 method for class 'bss' coef(object, ...)
object |
object of class bss. |
... |
further arguments to be passed to or from methods. |
Klaus Nordhausen
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 columnA<- 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
Comon's GAP criterion to evaluate the performance of an ICA algorithm.
ComonGAP(A, A.hat)ComonGAP(A, A.hat)
A |
The true square mixing matrix. |
A.hat |
The estimated square mixing matrix. |
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 , as is assumed by JADE and ics. However fastICA and
PearsonICA assume . Therefore matrices from those functions have to be transposed first.
The value of the Comon's GAP.
Klaus Nordhausen
Comon, P., (1994), Independent Component Analysis, A new concept?, Signal Processing, 36, 287–314.
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)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)
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.
data("CPPdata")data("CPPdata")
A data frame with 50000 observations on the following 4 variables.
Mic1the mixture recorded by the first microphone.
Mic2the mixture recorded by the second microphone.
Mic3the mixture recorded by the third microphone.
Mic4the mixture recorded by the fourth microphone.
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).
Ella Bingham
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>.
## 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)## 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)
This function jointly diagonalizes k real-valued square matrices by searching an orthogonal matrix in a deflation based manner.
djd(X, G = "max", r = 2, eps = 1e-06, maxiter = 500)djd(X, G = "max", r = 2, eps = 1e-06, maxiter = 500)
X |
an array containing the k p times p real valued matrices of dimension c(p, p, k). |
G |
criterion function used for the the algorithm. Options are |
r |
power value used if |
eps |
convergence tolerance. |
maxiter |
maximum number of iterations. |
Denote the square matrices as , . This algorithm searches then an orthogonal matrix W
so that is diagonal for all . If the commute then there is an exact solution. If not, the function
will perform an approximate joint diagonalization by maximizing where are the orthogonal vectors in W.
The function G can be choosen to be of the form or . If G="max" is chosen, the function G is of the form , 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 which plays hence a special role.
The matrix W
Klaus Nordhausen, Jari Miettinen
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>.
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.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.
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.
FG(X, weight = NULL, init = NULL, maxiter = 100, eps = 1e-06, na.action = na.fail)FG(X, weight = NULL, init = NULL, maxiter = 100, eps = 1e-06, na.action = na.fail)
X |
A matrix of k stacked pxp matrices with dimension c(kp,p) or an array with dimension c(p,p,k). |
weight |
A vector of length k to give weight to the different matrices, if NULL, all matrices have equal weight. |
init |
Initial value for the orthogonal matrix to be estimated, if NULL, the identity matrix is used. |
maxiter |
Maximum number of iterations. |
eps |
Convergence tolerance. |
na.action |
A function which indicates what should happen when the data contain 'NA's. Default is to fail. |
A list with the components
V |
An orthogonal matrix. |
D |
A stacked matrix with the diagonal matrices or an array with the diagonal matrices. The form of the output depends on the form of the input. |
iter |
The Fortran function returns also the number of iterations. |
Jari Miettinen
Flury, B. D. (1998), Common principal components and related models, Wiley, New York.
The FOBI method for independent component analysis (ICA). We assume that all components have different kurtosis values.
FOBI(X, na.action = na.fail)FOBI(X, na.action = na.fail)
X |
a numeric matrix. |
na.action |
A function which indicates what should happen when the data contain 'NA's. Default is to fail. |
A list with class 'bss' containing the following components:
W |
estimated unmixing matrix. |
EV |
eigenvectors of autocovariance matrix. |
Xmu |
the original mean of the data. |
S |
estimated sources as time series objected standardized to have mean 0 and unit variances. |
More general is the function ics in the ICS package.
Klaus Nordhausen
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>.
# 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)# 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)
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.
JADE(X, n.comp = NULL, eps = 1e-06, maxiter = 100, na.action = na.fail)JADE(X, n.comp = NULL, eps = 1e-06, maxiter = 100, na.action = na.fail)
X |
Numeric data matrix or dataframe. |
n.comp |
Number of components to extract. |
eps |
Convergence tolerance. |
maxiter |
Maximum number of iterations. |
na.action |
A function which indicates what should happen when the data contain 'NA's. Default is to fail. |
Some minor modifications were done when porting the function to R, and they are:
the model assumed here is . Therefore and have one row per observation. Note that this still differs from
the model definition in R of FastICA and PearsonICA but agrees with ics.
The whitening covariance matrix is divided by n-1 and not n (n = number of observations).
The initial value for the joint diagonalisation is always I.
The original eps would be .
It is also worth mentioning that the estimated independent components 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 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").
A list with class 'bss' containing the following components:
A |
The estimated mixing matrix. |
W |
The estimated unmixing matrix. |
S |
Dataframe with the estimated independent components. |
Xmu |
The location of the original data. |
Jean-Francois Cardoso. Ported to R by Klaus Nordhausen
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>.
# 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)# 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)
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.
k_JADE(X, k = 1, eps = 1e-06, maxiter = 100, na.action = na.fail)k_JADE(X, k = 1, eps = 1e-06, maxiter = 100, na.action = na.fail)
X |
Numeric data matrix or dataframe. |
k |
integer value between 1 and the number of columns of X. Default is 1. |
eps |
Convergence tolerance. |
maxiter |
Maximum number of iterations. |
na.action |
A function which indicates what should happen when the data contain 'NA's. Default is to fail. |
The order of the estimated components is fixed so that their fourth moments are in the decreasing order.
A list with class 'bss' containing the following components:
A |
The estimated mixing matrix. |
W |
The estimated unmixing matrix. |
S |
Matrix with the estimated independent components. |
Xmu |
The location of the original data. |
The function uses FOBI as initial estimate and frjd for the joint diagonalization.
Jari Miettinen
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>.
# 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)# 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)
Computes the Minimum Distance index MD to evaluate the performance of an ICA algorithm.
MD(W.hat, A)MD(W.hat, A)
W.hat |
The estimated square unmixing matrix W. |
A |
The true square mixing matrix A. |
where is a permutation matrix and 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 , as is assumed by JADE and ics. However fastICA and
PearsonICA assume . 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.
The value of the MD index.
Klaus Nordhausen
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>.
ComonGAP, SIR, amari.error, solve_LSAP
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)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)
The function can be used to compute several scatter matrices for the same data.
multscatter(scatterlist, X, toshape = TRUE)multscatter(scatterlist, X, toshape = TRUE)
scatterlist |
a vector with the names of the scatter matrices to be computed. Note that each of these functions should only return a matrix of size p times p. |
X |
the n times p data matrix for which the scatter should be computed. |
toshape |
logical, whether scatter matrices should be converted to shape matrices. If TRUE, all matrices will have determinant 1. |
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.
An array of dimension c(p,p,k) where k is the number of scatter matrices.
Klaus Nordhausen
# 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)# 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)
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.
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, ...)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, ...)
X |
a numeric matrix or a multivariate time series object of class |
K |
number of intervals to be used. |
Tau |
By default 0 which means covariance are computed of each time interval, if Tau is an integer > 0 then rather autocovariance matrices at lag Tau are used for the joint diagonaliation. |
n.cuts |
if NULL, then the time series is divided into K equally long intervals. To specify intervals n.cuts should be given in the form c(1,n.cut.1,...,n.cut.k, nrow(X)) to specify where to split the time series. |
eps |
maximum number of iterations for |
maxiter |
convergence tolerance for |
... |
further arguments to be passed to or from methods. |
The model assumes that the mean of the p-variate time series is constant but the variances change over time.
A list with class 'bss' containing the following components:
W |
estimated unmixing matrix. |
k |
the lag used for the autocovariance matrix. |
n.cut |
specifying the intervals where data is split |
K |
the number of intervals used |
S |
estimated sources as time series objected standardized to have mean 0 and that the variance of the sources are 1. |
Klaus Nordhausen
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>.
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)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)
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.
NSS.SD(X, ...) ## Default S3 method: NSS.SD(X, n.cut=NULL, ...) ## S3 method for class 'ts' NSS.SD(X, ...)NSS.SD(X, ...) ## Default S3 method: NSS.SD(X, n.cut=NULL, ...) ## S3 method for class 'ts' NSS.SD(X, ...)
X |
a numeric matrix or a multivariate time series object of class |
n.cut |
either an integer between 1 and nrow(X) or an vector of length 3 of the form c(1,n.cut,nrow(X)) to specify where to split the time series. If NULL, then c(1,floor(nrow(X)/2),nrow(X)) is used. |
... |
further arguments to be passed to or from methods. |
The model assumes that the mean of the p-variate time series is constant but the variances change over time.
A list with class 'bss' containing the following components:
W |
estimated unmixing matrix. |
EV |
eigenvalues from the eigenvalue-eigenvector decomposition. |
n.cut |
specifying the intervals where data is split |
S |
estimated sources as time series objected standardized to have mean 0 and that the sources in the first interval are 1. |
Klaus Nordhausen
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>.
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)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)
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.
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, ...)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, ...)
X |
a numeric matrix or a multivariate time series object of class |
K |
number of intervals to be used. |
Tau |
Lags for the autovariance matrices to be computed within each interval. |
n.cuts |
if NULL, then the time series is divided into K equally long intervals. To specify intervals n.cuts should be given in the form c(1,n.cut.1,...,n.cut.k, nrow(X)) to specify where to split the time series. |
eps |
maximum number of iterations for |
maxiter |
convergence tolerance for |
... |
further arguments to be passed to or from methods. |
The model assumes that the mean of the p-variate time series is constant but the variances change over time.
A list with class 'bss' containing the following components:
W |
estimated unmixing matrix. |
k |
the lags used for the autocovariance matrix used in each interval. |
n.cut |
specifying the intervals where data is split |
K |
the number of intervals used |
S |
estimated sources as time series objected standardized to have mean 0 and that the sources 1. |
Klaus Nordhausen
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>.
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)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)
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.
## S3 method for class 'bss' plot(x, ...)## S3 method for class 'bss' plot(x, ...)
x |
object of class bss. |
... |
further arguments to be passed to or from methods. |
Klaus Nordhausen
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)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)
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.
## S3 method for class 'bss' print(x, ...)## S3 method for class 'bss' print(x, ...)
x |
object of class bss. |
... |
further arguments to be passed to or from methods. |
Klaus Nordhausen
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.
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)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)
X |
A matrix of k stacked pxp matrices with dimension c(kp,p) or an array with dimension c(p,p,k). In case of |
weight |
A vector of length k to give weight to the different matrices, if NULL, all matrices have equal weight |
eps |
Convergence tolerance. |
maxiter |
Maximum number of iterations. |
na.action |
A function which indicates what should happen when the data contain 'NA's. Default is to fail. |
Denote the square matrices as , . The algorithm searches an orthogonal matrix so that is diagonal for all . If the commute then there is an exact solution. Otherwise, the function will perform an approximate joint diagonalization by trying to make the 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.
A list with the components
V |
An orthogonal matrix. |
D |
A stacked matrix with the diagonal matrices or an array with the diagonal matrices. The form of the output depends on the form of the input. |
iter |
The |
Jean-Francois Cardoso. Ported to R by Klaus Nordhausen. C code by Jari Miettinen
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>.
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)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)
Computes the signal to interference ratio between true and estimated signals
SIR(S, S.hat)SIR(S, S.hat)
S |
Matrix or dataframe with the true numeric signals. |
S.hat |
Matrix or dataframe with the estimated numeric signals. |
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.
The value of the signal to interference ratio.
Klaus Nordhausen
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.
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)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)
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.
SOBI(X, ...) ## Default S3 method: SOBI(X, k=12, method="frjd", eps = 1e-06, maxiter = 100, ...) ## S3 method for class 'ts' SOBI(X, ...)SOBI(X, ...) ## Default S3 method: SOBI(X, k=12, method="frjd", eps = 1e-06, maxiter = 100, ...) ## S3 method for class 'ts' SOBI(X, ...)
X |
a numeric matrix or a multivariate time series object of class |
k |
if a single integer, then the lags 1:k are used, if an integer vector, then these are used as the lags. |
method |
method to use for the joint diagonalization, options are |
.
eps |
convergence tolerance. |
maxiter |
maximum number of iterations. |
... |
further arguments to be passed to or from methods. |
The order of the estimated components is fixed so that the sums of squared autocovariances are in the decreasing order.
A list with class 'bss' containing the following components:
W |
estimated unmixing matrix. |
k |
lags used. |
method |
method used for the joint diagonalization. |
S |
estimated sources as time series objected standardized to have mean 0 and unit variances. |
Klaus Nordhausen
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>.
# 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)# 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)