# NOT RUN {
## Example
set.seed(1)
# Generate data from a block-wise variance model
d <- 2
m <- 10
n <- 5000
group_index <- rep(c(1,2), each=n)
partition_index <- rep(rep(1:m, each=n/m), 2)
S <- matrix(NA, 2*n, d)
H <- matrix(NA, 2*n, d)
for(i in unique(group_index)){
varH <- abs(rnorm(d))/4
H[group_index==i, ] <- matrix(rnorm(d*n)*rep(varH, each=n), n, d)
for(j in unique(partition_index[group_index==i])){
varS <- abs(rnorm(d))
index <- partition_index==j & group_index==i
S[index,] <- matrix(rnorm(d*n/m)*rep(varS, each=n/m),
n/m, d)
}
}
A <- matrix(rnorm(d^2), d, d)
A <- A%*%t(A)
X <- t(A%*%t(S+H))
# Apply coroICA
res <- coroICA(X, group_index, partition_index, pairing="allpairs", rank_components=TRUE)
# Compare results
par(mfrow=c(2,2))
plot((S+H)[,1], type="l", main="true source 1", ylab="S+H")
plot(res$Shat[,1], type="l", main="estimated source 1", ylab="Shat")
plot((S+H)[,2], type="l", main="true source 2", ylab="S+H")
plot(res$Shat[,2], type="l", main="estimated source 2", ylab="Shat")
cor(res$Shat, S+H)
# }
Run the code above in your browser using DataLab