## Artificial data set 
 	n <- 100						# number of samples
 	p <- 30						# dimension of data
 	t.k <- 4						# number of factor variables
 	dgs <- 0.70						# degree of sparseness, which lies in 0 to 1
 	t.sigma <- seq(0.8, 1.5, length=t.k)	# variances of factor variables
 	t.lambda <- 0.05					# noise variance


  ## Construction of orthogonal sparse loading matrix
	## Sparsity configuration  
 	t.load <- matrix( sample(c(0,1), p*t.k, replace=TRUE, prob=c(dgs, 1-dgs)), nrow=p, ncol=t.k)
 	id <- which( t.load!=0, arr.ind = TRUE )	

 	## Sequential generation of loading vector orthogonal to the preceding ones 
	for(i in 2:t.k){
		prid <- c(1:i) 
		set1 <- which(is.element(id[ ,2], prid)==TRUE, arr.ind = FALSE)	
		set2 <- which(is.element(id[ ,2], i)==TRUE, arr.ind = FALSE)	
		set1 <- setdiff(set1, set2)						
		cset <- intersect(id[set1, 1], id[set2, 1])
		if(length(cset)==1){ t.load[cset, i] <- 0 }
		if(length(cset)>1){
			rmat <- matrix( rnorm(length(cset)*n), nrow=n, ncol=length(cset) )
			rmat <- t(rmat)%*%rmat
			idtmp <- prid[1:i-1]
			q <- solve( t(t.load[cset, idtmp])%*%t.load[cset, idtmp] )
			prj.mat <- (diag(rep(1, length=length(cset)))- t.load[cset, idtmp]%*%q%*%t(t.load[cset, idtmp]))
			res <- eigen( prj.mat%*%rmat%*%t(prj.mat) )
			t.load[cset, i] <- t(prj.mat)%*%res$vectors[ ,1]		
		}
 	} 
 	t.load <- t.load%*%diag(1/sqrt(colSums(t.load*t.load)))
 
 ## generation of data points

   ## Gaussian noise
   noise <- matrix(rnorm(n*p, mean=0, sd=sqrt(t.lambda)), ncol=n, nrow=p)
   ## Factor variables
   factor <- diag(sqrt(t.sigma))%*%matrix(rnorm(n*t.k, mean=0, sd=1), ncol=n, nrow=t.k)
   ## Data 
   data <- t.load%*%factor + noise