Title: | Independent Component Analysis |
---|---|
Description: | Independent Component Analysis (ICA) using various algorithms: FastICA, Information-Maximization (Infomax), and Joint Approximate Diagonalization of Eigenmatrices (JADE). |
Authors: | Nathaniel E. Helwig <[email protected]> |
Maintainer: | Nathaniel E. Helwig <[email protected]> |
License: | GPL (>= 2) |
Version: | 1.0-3 |
Built: | 2024-12-19 06:43:40 UTC |
Source: | CRAN |
The Amari-Cichocki-Yang (ACY) error is an asymmetric measure of dissimilarity between two nonsingular matrices X
and Y
. The ACY error: (a) is invariant to permutation and rescaling of the columns of X
and Y
, (b) ranges between 0 and n-1
, and (c) equals 0 if and only if X
and Y
are identical up to column permutations and rescalings.
acy(X,Y)
acy(X,Y)
X |
Nonsingular matrix of dimension |
Y |
Nonsingular matrix of dimension |
The ACY error is defined as
where .
Returns a scalar (the ACY error).
If Y
is singular, function will produce an error.
Nathaniel E. Helwig <[email protected]>
Amari, S., Cichocki, A., & Yang, H.H. (1996). A new learning algorithm for blind signal separation. In D. S. Touretzky, M. C. Mozer, and M. E. Hasselmo (Eds.), Advances in Neural Information Processing Systems, 8. Cambridge, MA: MIT Press.
########## EXAMPLE ########## set.seed(1) X <- matrix(runif(16),4,4) Y <- matrix(runif(16),4,4) Z <- X[,c(3,1,2,4)]%*%diag(1:4) acy(X,Y) acy(X,Z)
########## EXAMPLE ########## set.seed(1) X <- matrix(runif(16),4,4) Y <- matrix(runif(16),4,4) Z <- X[,c(3,1,2,4)]%*%diag(1:4) acy(X,Y) acy(X,Z)
Computes ICA decomposition using Hyvarinen's (1999) FastICA algorithm, Bell and Sejnowski's (1995) Information-Maximization (Infomax) algorithm, or Cardoso and Souloumiac's (1993, 1996) Joint Approximate Diagonalization of Eigenmatrices (JADE) algorithm.
ica(X, nc, method = c("fast", "imax", "jade"), ...)
ica(X, nc, method = c("fast", "imax", "jade"), ...)
X |
Data matrix with |
nc |
Number of components to extract. |
method |
Method for decomposition. |
... |
Additional arguments to be passed to other ICA functions (see Details). |
ICA Model
The ICA model can be written as X = tcrossprod(S, M) + E
, where S
contains the source signals, M
is the mixing matrix, and E
contains the noise signals. Columns of X
are assumed to have zero mean. The goal is to find the unmixing matrix W
such that columns of S = tcrossprod(X, W)
are independent as possible.
Whitening
Without loss of generality, we can write M = P %*% R
where P
is a tall matrix and R
is an orthogonal rotation matrix. Letting Q
denote the pseudoinverse of P
, we can whiten the data using Y = tcrossprod(X, Q)
. The goal is to find the orthongal rotation matrix R
such that the source signal estimates S = Y %*% R
are as independent as possible. Note that W = crossprod(R, Q)
.
Method
This is a wrapper function for the functions icafast
, icaimax
, or icajade
. See the corresponding function for details on the method, as well as the available arguments (handled by the ...
argument).
S |
Matrix of source signal estimates ( |
M |
Estimated mixing matrix. |
W |
Estimated unmixing matrix ( |
Y |
Whitened data matrix. |
Q |
Whitening matrix. |
R |
Orthogonal rotation matrix. |
vafs |
Variance-accounted-for by each component. |
iter |
Number of algorithm iterations. |
converged |
Logical indicating if algorithm converged. |
... |
Other arguments (if |
Nathaniel E. Helwig <[email protected]>
Bell, A.J. & Sejnowski, T.J. (1995). An information-maximization approach to blind separation and blind deconvolution. Neural Computation, 7(6), 1129-1159. doi:10.1162/neco.1995.7.6.1129
Cardoso, J.F., & Souloumiac, A. (1993). Blind beamforming for non-Gaussian signals. IEE Proceedings-F, 140(6), 362-370. doi:10.1049/ip-f-2.1993.0054
Cardoso, J.F., & Souloumiac, A. (1996). Jacobi angles for simultaneous diagonalization. SIAM Journal on Matrix Analysis and Applications, 17(1), 161-164. doi:10.1137/S0895479893259546
Helwig, N.E. & Hong, S. (2013). A critique of Tensor Probabilistic Independent Component Analysis: Implications and recommendations for multi-subject fMRI data analysis. Journal of Neuroscience Methods, 213(2), 263-273. doi:10.1016/j.jneumeth.2012.12.009
Hyvarinen, A. (1999). Fast and robust fixed-point algorithms for independent component analysis. IEEE Transactions on Neural Networks, 10(3), 626-634. doi:10.1109/72.761722
icafast
for ICA via FastICA
icaimax
for ICA via Infomax
icajade
for ICA via JADE
########## EXAMPLE 1 ########## # generate noiseless data (p == r) set.seed(123) nobs <- 1000 Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs)) Bmat <- matrix(2 * runif(4), nrow = 2, ncol = 2) Xmat <- tcrossprod(Amat, Bmat) # ICA via different algorithms imod.fast <- ica(Xmat, nc = 2) imod.imax <- ica(Xmat, nc = 2, method = "imax") imod.jade <- ica(Xmat, nc = 2, method = "jade") # compare mixing matrix recovery acy(Bmat, imod.fast$M) acy(Bmat, imod.imax$M) acy(Bmat, imod.jade$M) # compare source signal recovery cor(Amat, imod.fast$S) cor(Amat, imod.imax$S) cor(Amat, imod.jade$S) ########## EXAMPLE 2 ########## # generate noiseless data (p != r) set.seed(123) nobs <- 1000 Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs)) Bmat <- matrix(2 * runif(200), nrow = 100, ncol = 2) Xmat <- tcrossprod(Amat, Bmat) # ICA via different algorithms imod.fast <- ica(Xmat, nc = 2) imod.imax <- ica(Xmat, nc = 2, method = "imax") imod.jade <- ica(Xmat, nc = 2, method = "jade") # compare source signal recovery cor(Amat, imod.fast$S) cor(Amat, imod.imax$S) cor(Amat, imod.jade$S) ########## EXAMPLE 3 ########## # generate noisy data (p != r) set.seed(123) nobs <- 1000 Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs)) Bmat <- matrix(2 * runif(200), 100, 2) Emat <- matrix(rnorm(10^5), nrow = 1000, ncol = 100) Xmat <- tcrossprod(Amat,Bmat) + Emat # ICA via different algorithms imod.fast <- ica(Xmat, nc = 2) imod.imax <- ica(Xmat, nc = 2, method = "imax") imod.jade <- ica(Xmat, nc = 2, method = "jade") # compare source signal recovery cor(Amat, imod.fast$S) cor(Amat, imod.imax$S) cor(Amat, imod.jade$S)
########## EXAMPLE 1 ########## # generate noiseless data (p == r) set.seed(123) nobs <- 1000 Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs)) Bmat <- matrix(2 * runif(4), nrow = 2, ncol = 2) Xmat <- tcrossprod(Amat, Bmat) # ICA via different algorithms imod.fast <- ica(Xmat, nc = 2) imod.imax <- ica(Xmat, nc = 2, method = "imax") imod.jade <- ica(Xmat, nc = 2, method = "jade") # compare mixing matrix recovery acy(Bmat, imod.fast$M) acy(Bmat, imod.imax$M) acy(Bmat, imod.jade$M) # compare source signal recovery cor(Amat, imod.fast$S) cor(Amat, imod.imax$S) cor(Amat, imod.jade$S) ########## EXAMPLE 2 ########## # generate noiseless data (p != r) set.seed(123) nobs <- 1000 Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs)) Bmat <- matrix(2 * runif(200), nrow = 100, ncol = 2) Xmat <- tcrossprod(Amat, Bmat) # ICA via different algorithms imod.fast <- ica(Xmat, nc = 2) imod.imax <- ica(Xmat, nc = 2, method = "imax") imod.jade <- ica(Xmat, nc = 2, method = "jade") # compare source signal recovery cor(Amat, imod.fast$S) cor(Amat, imod.imax$S) cor(Amat, imod.jade$S) ########## EXAMPLE 3 ########## # generate noisy data (p != r) set.seed(123) nobs <- 1000 Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs)) Bmat <- matrix(2 * runif(200), 100, 2) Emat <- matrix(rnorm(10^5), nrow = 1000, ncol = 100) Xmat <- tcrossprod(Amat,Bmat) + Emat # ICA via different algorithms imod.fast <- ica(Xmat, nc = 2) imod.imax <- ica(Xmat, nc = 2, method = "imax") imod.jade <- ica(Xmat, nc = 2, method = "jade") # compare source signal recovery cor(Amat, imod.fast$S) cor(Amat, imod.imax$S) cor(Amat, imod.jade$S)
Computes ICA decomposition using Hyvarinen's (1999) FastICA algorithm with various options.
icafast(X, nc, center = TRUE, maxit = 100, tol = 1e-6, Rmat = diag(nc), alg = "par", fun = "logcosh", alpha = 1)
icafast(X, nc, center = TRUE, maxit = 100, tol = 1e-6, Rmat = diag(nc), alg = "par", fun = "logcosh", alpha = 1)
X |
Data matrix with |
nc |
Number of components to extract. |
center |
If |
maxit |
Maximum number of algorithm iterations to allow. |
tol |
Convergence tolerance. |
Rmat |
Initial estimate of the |
alg |
Algorithm to use: |
fun |
Contrast function to use for negentropy approximation: |
alpha |
Tuning parameter for "logcosh" contrast function (1 <= |
ICA Model
The ICA model can be written as X = tcrossprod(S, M) + E
, where S
contains the source signals, M
is the mixing matrix, and E
contains the noise signals. Columns of X
are assumed to have zero mean. The goal is to find the unmixing matrix W
such that columns of S = tcrossprod(X, W)
are independent as possible.
Whitening
Without loss of generality, we can write M = P %*% R
where P
is a tall matrix and R
is an orthogonal rotation matrix. Letting Q
denote the pseudoinverse of P
, we can whiten the data using Y = tcrossprod(X, Q)
. The goal is to find the orthongal rotation matrix R
such that the source signal estimates S = Y %*% R
are as independent as possible. Note that W = crossprod(R, Q)
.
FastICA
The FastICA algorithm finds the orthogonal rotation matrix R
that (approximately) maximizes the negentropy of the estimated source signals. Negentropy is approximated using
where E denotes the expectation, G is the contrast function, and z is a standard normal variable. See Hyvarinen (1999) or Helwig and Hong (2013) for specifics of fixed-point algorithm.
S |
Matrix of source signal estimates ( |
M |
Estimated mixing matrix. |
W |
Estimated unmixing matrix ( |
Y |
Whitened data matrix. |
Q |
Whitening matrix. |
R |
Orthogonal rotation matrix. |
vafs |
Variance-accounted-for by each component. |
iter |
Number of algorithm iterations. |
alg |
Algorithm used (same as input). |
fun |
Contrast function (same as input). |
alpha |
Tuning parameter (same as input). |
converged |
Logical indicating if algorithm converged. |
Nathaniel E. Helwig <[email protected]>
Helwig, N.E. & Hong, S. (2013). A critique of Tensor Probabilistic Independent Component Analysis: Implications and recommendations for multi-subject fMRI data analysis. Journal of Neuroscience Methods, 213(2), 263-273. doi:10.1016/j.jneumeth.2012.12.009
Hyvarinen, A. (1999). Fast and robust fixed-point algorithms for independent component analysis. IEEE Transactions on Neural Networks, 10(3), 626-634. doi:10.1109/72.761722
icaimax
for ICA via Infomax
icajade
for ICA via JADE
########## EXAMPLE 1 ########## # generate noiseless data (p == r) set.seed(123) nobs <- 1000 Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs)) Bmat <- matrix(2 * runif(4), nrow = 2, ncol = 2) Xmat <- tcrossprod(Amat, Bmat) # ICA via FastICA with 2 components imod <- icafast(Xmat, nc = 2) acy(Bmat, imod$M) cor(Amat, imod$S) ########## EXAMPLE 2 ########## # generate noiseless data (p != r) set.seed(123) nobs <- 1000 Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs)) Bmat <- matrix(2 * runif(200), nrow = 100, ncol = 2) Xmat <- tcrossprod(Amat, Bmat) # ICA via FastICA with 2 components imod <- icafast(Xmat, nc = 2) cor(Amat, imod$S) ########## EXAMPLE 3 ########## # generate noisy data (p != r) set.seed(123) nobs <- 1000 Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs)) Bmat <- matrix(2 * runif(200), 100, 2) Emat <- matrix(rnorm(10^5), nrow = 1000, ncol = 100) Xmat <- tcrossprod(Amat,Bmat) + Emat # ICA via FastICA with 2 components imod <- icafast(Xmat, nc = 2) cor(Amat, imod$S)
########## EXAMPLE 1 ########## # generate noiseless data (p == r) set.seed(123) nobs <- 1000 Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs)) Bmat <- matrix(2 * runif(4), nrow = 2, ncol = 2) Xmat <- tcrossprod(Amat, Bmat) # ICA via FastICA with 2 components imod <- icafast(Xmat, nc = 2) acy(Bmat, imod$M) cor(Amat, imod$S) ########## EXAMPLE 2 ########## # generate noiseless data (p != r) set.seed(123) nobs <- 1000 Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs)) Bmat <- matrix(2 * runif(200), nrow = 100, ncol = 2) Xmat <- tcrossprod(Amat, Bmat) # ICA via FastICA with 2 components imod <- icafast(Xmat, nc = 2) cor(Amat, imod$S) ########## EXAMPLE 3 ########## # generate noisy data (p != r) set.seed(123) nobs <- 1000 Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs)) Bmat <- matrix(2 * runif(200), 100, 2) Emat <- matrix(rnorm(10^5), nrow = 1000, ncol = 100) Xmat <- tcrossprod(Amat,Bmat) + Emat # ICA via FastICA with 2 components imod <- icafast(Xmat, nc = 2) cor(Amat, imod$S)
Computes ICA decomposition using Bell and Sejnowski's (1995) Information-Maximization (Infomax) approach with various options.
icaimax(X, nc, center = TRUE, maxit = 100, tol = 1e-6, Rmat = diag(nc), alg = "newton", fun = "tanh", signs = rep(1, nc), signswitch = TRUE, rate = 1, rateanneal = NULL)
icaimax(X, nc, center = TRUE, maxit = 100, tol = 1e-6, Rmat = diag(nc), alg = "newton", fun = "tanh", signs = rep(1, nc), signswitch = TRUE, rate = 1, rateanneal = NULL)
X |
Data matrix with |
nc |
Number of components to extract. |
center |
If |
maxit |
Maximum number of algorithm iterations to allow. |
tol |
Convergence tolerance. |
Rmat |
Initial estimate of the |
alg |
Algorithm to use: |
fun |
Nonlinear (squashing) function to use for algorithm: |
signs |
Vector of length |
signswitch |
If |
rate |
Learing rate for gradient descent algorithm. Ignored if |
rateanneal |
Annealing angle and proportion for gradient descent learing rate (see Details). Ignored if |
ICA Model
The ICA model can be written as X = tcrossprod(S, M) + E
, where S
contains the source signals, M
is the mixing matrix, and E
contains the noise signals. Columns of X
are assumed to have zero mean. The goal is to find the unmixing matrix W
such that columns of S = tcrossprod(X, W)
are independent as possible.
Whitening
Without loss of generality, we can write M = P %*% R
where P
is a tall matrix and R
is an orthogonal rotation matrix. Letting Q
denote the pseudoinverse of P
, we can whiten the data using Y = tcrossprod(X, Q)
. The goal is to find the orthongal rotation matrix R
such that the source signal estimates S = Y %*% R
are as independent as possible. Note that W = crossprod(R, Q)
.
Infomax
The Infomax approach finds the orthogonal rotation matrix R
that (approximately) maximizes the joint entropy of a nonlinear function of the estimated source signals. See Bell and Sejnowski (1995) and Helwig and Hong (2013) for specifics of algorithms.
S |
Matrix of source signal estimates ( |
M |
Estimated mixing matrix. |
W |
Estimated unmixing matrix ( |
Y |
Whitened data matrix. |
Q |
Whitening matrix. |
R |
Orthogonal rotation matrix. |
vafs |
Variance-accounted-for by each component. |
iter |
Number of algorithm iterations. |
alg |
Algorithm used (same as input). |
fun |
Contrast function (same as input). |
signs |
Component signs (same as input). |
rate |
Learning rate (same as input). |
converged |
Logical indicating if algorithm converged. |
Nathaniel E. Helwig <[email protected]>
Bell, A.J. & Sejnowski, T.J. (1995). An information-maximization approach to blind separation and blind deconvolution. Neural Computation, 7(6), 1129-1159. doi:10.1162/neco.1995.7.6.1129
Helwig, N.E. & Hong, S. (2013). A critique of Tensor Probabilistic Independent Component Analysis: Implications and recommendations for multi-subject fMRI data analysis. Journal of Neuroscience Methods, 213(2), 263-273. doi:10.1016/j.jneumeth.2012.12.009
icafast
for FastICA
icajade
for ICA via JADE
########## EXAMPLE 1 ########## # generate noiseless data (p == r) set.seed(123) nobs <- 1000 Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs)) Bmat <- matrix(2 * runif(4), nrow = 2, ncol = 2) Xmat <- tcrossprod(Amat, Bmat) # ICA via Infomax with 2 components imod <- icaimax(Xmat, nc = 2) acy(Bmat, imod$M) cor(Amat, imod$S) ########## EXAMPLE 2 ########## # generate noiseless data (p != r) set.seed(123) nobs <- 1000 Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs)) Bmat <- matrix(2 * runif(200), nrow = 100, ncol = 2) Xmat <- tcrossprod(Amat, Bmat) # ICA via Infomax with 2 components imod <- icaimax(Xmat, nc = 2) cor(Amat, imod$S) ########## EXAMPLE 3 ########## # generate noisy data (p != r) set.seed(123) nobs <- 1000 Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs)) Bmat <- matrix(2 * runif(200), 100, 2) Emat <- matrix(rnorm(10^5), nrow = 1000, ncol = 100) Xmat <- tcrossprod(Amat,Bmat) + Emat # ICA via Infomax with 2 components imod <- icaimax(Xmat, nc = 2) cor(Amat, imod$S)
########## EXAMPLE 1 ########## # generate noiseless data (p == r) set.seed(123) nobs <- 1000 Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs)) Bmat <- matrix(2 * runif(4), nrow = 2, ncol = 2) Xmat <- tcrossprod(Amat, Bmat) # ICA via Infomax with 2 components imod <- icaimax(Xmat, nc = 2) acy(Bmat, imod$M) cor(Amat, imod$S) ########## EXAMPLE 2 ########## # generate noiseless data (p != r) set.seed(123) nobs <- 1000 Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs)) Bmat <- matrix(2 * runif(200), nrow = 100, ncol = 2) Xmat <- tcrossprod(Amat, Bmat) # ICA via Infomax with 2 components imod <- icaimax(Xmat, nc = 2) cor(Amat, imod$S) ########## EXAMPLE 3 ########## # generate noisy data (p != r) set.seed(123) nobs <- 1000 Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs)) Bmat <- matrix(2 * runif(200), 100, 2) Emat <- matrix(rnorm(10^5), nrow = 1000, ncol = 100) Xmat <- tcrossprod(Amat,Bmat) + Emat # ICA via Infomax with 2 components imod <- icaimax(Xmat, nc = 2) cor(Amat, imod$S)
Computes ICA decomposition using Cardoso and Souloumiac's (1993, 1996) Joint Approximate Diagonalization of Eigenmatrices (JADE) approach.
icajade(X, nc, center = TRUE, maxit = 100, tol = 1e-6, Rmat = diag(nc))
icajade(X, nc, center = TRUE, maxit = 100, tol = 1e-6, Rmat = diag(nc))
X |
Data matrix with |
nc |
Number of components to extract. |
center |
If |
maxit |
Maximum number of algorithm iterations to allow. |
tol |
Convergence tolerance. |
Rmat |
Initial estimate of the |
ICA Model
The ICA model can be written as X = tcrossprod(S, M) + E
, where S
contains the source signals, M
is the mixing matrix, and E
contains the noise signals. Columns of X
are assumed to have zero mean. The goal is to find the unmixing matrix W
such that columns of S = tcrossprod(X, W)
are independent as possible.
Whitening
Without loss of generality, we can write M = P %*% R
where P
is a tall matrix and R
is an orthogonal rotation matrix. Letting Q
denote the pseudoinverse of P
, we can whiten the data using Y = tcrossprod(X, Q)
. The goal is to find the orthongal rotation matrix R
such that the source signal estimates S = Y %*% R
are as independent as possible. Note that W = crossprod(R, Q)
.
JADE
The JADE approach finds the orthogonal rotation matrix R
that (approximately) diagonalizes the cumulant array of the source signals. See Cardoso and Souloumiac (1993,1996) and Helwig and Hong (2013) for specifics of the JADE algorithm.
S |
Matrix of source signal estimates ( |
M |
Estimated mixing matrix. |
W |
Estimated unmixing matrix ( |
Y |
Whitened data matrix. |
Q |
Whitening matrix. |
R |
Orthogonal rotation matrix. |
vafs |
Variance-accounted-for by each component. |
iter |
Number of algorithm iterations. |
converged |
Logical indicating if algorithm converged. |
Nathaniel E. Helwig <[email protected]>
Cardoso, J.F., & Souloumiac, A. (1993). Blind beamforming for non-Gaussian signals. IEE Proceedings-F, 140(6), 362-370. doi:10.1049/ip-f-2.1993.0054
Cardoso, J.F., & Souloumiac, A. (1996). Jacobi angles for simultaneous diagonalization. SIAM Journal on Matrix Analysis and Applications, 17(1), 161-164. doi:10.1137/S0895479893259546
Helwig, N.E. & Hong, S. (2013). A critique of Tensor Probabilistic Independent Component Analysis: Implications and recommendations for multi-subject fMRI data analysis. Journal of Neuroscience Methods, 213(2), 263-273. doi:10.1016/j.jneumeth.2012.12.009
icafast
for FastICA
icaimax
for ICA via Infomax
########## EXAMPLE 1 ########## # generate noiseless data (p == r) set.seed(123) nobs <- 1000 Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs)) Bmat <- matrix(2 * runif(4), nrow = 2, ncol = 2) Xmat <- tcrossprod(Amat, Bmat) # ICA via JADE with 2 components imod <- icajade(Xmat, nc = 2) acy(Bmat, imod$M) cor(Amat, imod$S) ########## EXAMPLE 2 ########## # generate noiseless data (p != r) set.seed(123) nobs <- 1000 Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs)) Bmat <- matrix(2 * runif(200), nrow = 100, ncol = 2) Xmat <- tcrossprod(Amat, Bmat) # ICA via JADE with 2 components imod <- icajade(Xmat, nc = 2) cor(Amat, imod$S) ########## EXAMPLE 3 ########## # generate noisy data (p != r) set.seed(123) nobs <- 1000 Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs)) Bmat <- matrix(2 * runif(200), 100, 2) Emat <- matrix(rnorm(10^5), nrow = 1000, ncol = 100) Xmat <- tcrossprod(Amat,Bmat) + Emat # ICA via JADE with 2 components imod <- icajade(Xmat, nc = 2) cor(Amat, imod$S)
########## EXAMPLE 1 ########## # generate noiseless data (p == r) set.seed(123) nobs <- 1000 Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs)) Bmat <- matrix(2 * runif(4), nrow = 2, ncol = 2) Xmat <- tcrossprod(Amat, Bmat) # ICA via JADE with 2 components imod <- icajade(Xmat, nc = 2) acy(Bmat, imod$M) cor(Amat, imod$S) ########## EXAMPLE 2 ########## # generate noiseless data (p != r) set.seed(123) nobs <- 1000 Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs)) Bmat <- matrix(2 * runif(200), nrow = 100, ncol = 2) Xmat <- tcrossprod(Amat, Bmat) # ICA via JADE with 2 components imod <- icajade(Xmat, nc = 2) cor(Amat, imod$S) ########## EXAMPLE 3 ########## # generate noisy data (p != r) set.seed(123) nobs <- 1000 Amat <- cbind(icasamp("a", "rnd", nobs), icasamp("b", "rnd", nobs)) Bmat <- matrix(2 * runif(200), 100, 2) Emat <- matrix(rnorm(10^5), nrow = 1000, ncol = 100) Xmat <- tcrossprod(Amat,Bmat) + Emat # ICA via JADE with 2 components imod <- icajade(Xmat, nc = 2) cor(Amat, imod$S)
Plot density (pdf) and kurtosis for the 18 source signal distributions used in Bach and Jordan (2002); see icasamp
for more information.
icaplot(xseq = seq(-2,2,length.out=500), xlab = "", ylab = "", lty = 1, lwd = 1, col = "black", ...)
icaplot(xseq = seq(-2,2,length.out=500), xlab = "", ylab = "", lty = 1, lwd = 1, col = "black", ...)
xseq |
Sequence of ordered data values for plotting density. |
xlab |
X-axis label for plot (default is no label). |
ylab |
Y-axis label for plot (default is no label). |
lty |
Line type for each density (scalar or vector of length 18). |
lwd |
Line width for each density (scalar or vector of length 18). |
col |
Line color for each density (scalar or vector of length 18). |
... |
Optional inputs for |
Produces a plot with NULL
return value.
Nathaniel E. Helwig <[email protected]>
Bach, F.R. (2002). kernel-ica. MATLAB toolbox (ver 1.2) http://www.di.ens.fr/~fbach/kernel-ica/.
Bach, F.R. & Jordan, M.I. (2002). Kernel independent component analysis. Journal of Machine Learning Research, 3, 1-48.
## Not run: ########## EXAMPLE ########## quartz(height=9,width=7) par(mar=c(3,3,3,3)) icaplot() ## End(Not run)
## Not run: ########## EXAMPLE ########## quartz(height=9,width=7) par(mar=c(3,3,3,3)) icaplot() ## End(Not run)
Sample observations from the 18 source signal distributions used in Bach and Jordan (2002). Can also return density values and kurtosis for each distribution. Use icaplot
to plot distributions.
icasamp(dname, query = c("rnd","pdf","kur"), nsamp = NULL, data = NULL)
icasamp(dname, query = c("rnd","pdf","kur"), nsamp = NULL, data = NULL)
dname |
Distribution name: letter "a" through "r" (see Bach & Jordan, 2002). |
query |
What to return: |
nsamp |
Number of observations to sample. Only used if |
data |
Data values for density evaluation. Only used if |
Inspired by usr_distrib.m
from Bach's (2002) kernel-ica
MATLAB toolbox.
If query="rnd"
, returns random sample of size nsamp
.
If query="pdf"
, returns density for input data
.
If query="kur"
, returns kurtosis of distribution.
Nathaniel E. Helwig <[email protected]>
Bach, F.R. (2002). kernel-ica. MATLAB toolbox (ver 1.2) http://www.di.ens.fr/~fbach/kernel-ica/.
Bach, F.R. & Jordan, M.I. (2002). Kernel independent component analysis. Journal of Machine Learning Research, 3, 1-48.
########## EXAMPLE ########## # sample 1000 observations from distribution "f" set.seed(123) mysamp <- icasamp("f","rnd",nsamp=1000) xr <- range(mysamp) hist(mysamp,freq=FALSE,ylim=c(0,.8),breaks=sqrt(1000)) # evaluate density of distribution "f" xseq <- seq(-5,5,length.out=1000) mypdf <- icasamp("f","pdf",data=xseq) lines(xseq,mypdf) # evaluate kurtosis of distribution "f" icasamp("f","kur")
########## EXAMPLE ########## # sample 1000 observations from distribution "f" set.seed(123) mysamp <- icasamp("f","rnd",nsamp=1000) xr <- range(mysamp) hist(mysamp,freq=FALSE,ylim=c(0,.8),breaks=sqrt(1000)) # evaluate density of distribution "f" xseq <- seq(-5,5,length.out=1000) mypdf <- icasamp("f","pdf",data=xseq) lines(xseq,mypdf) # evaluate kurtosis of distribution "f" icasamp("f","kur")