### function to generate dataset
simSNP<-function(n.obs=500, n.snp=50, maf=0.25, beta0=0,beta=NULL,
	rand=NA, mean=0, sd=10){
	X<-data.frame(matrix(NA, nrow=n.obs, ncol=n.snp))
	colnames(X)<-paste('SNP', 1:n.snp, sep='')
	if(!is.na(rand)) set.seed(rand)
	if(length(maf)==1){
		X <- data.frame(sapply(1:n.snp, function(z) 
			       as.factor(sample(c(0,1,2), n.obs, replace=TRUE, 
                            prob=c({1-maf}**2, 2*{1-maf}*maf, maf**2)))))
	}
	else{
		#if(length(maf)!=n.snp) stop('length of "maf" and "n.snp" have to be equal') 
		X <- data.frame(sapply(1:n.snp, function(z)
			         sample(c(0,1,2), n.obs, replace=TRUE,  
                        prob=c({1-maf[z]}**2, 2*{1-maf[z]}*maf[z], maf[z]**2))))
	}
	colnames(X)<-paste('SNP', 1:n.snp, sep='')
        if (is.null(beta)) 
          {
          beta<-rep(0,n.snp)
          }
        linpred<-as.matrix(X)%*%matrix(beta,ncol=1)
	Y <- as.factor(rbinom(n.obs,size=1,prob=plogis(linpred)))
  for (j in 1:n.snp)
   {
   X[,j]<-as.factor(X[,j])
   }
  return(data.frame(cbind(y=Y, X)))

}

