R/Rutils.R

Defines functions generateSparsePoints dat3sample dat2sample dat1sample generate2Dtransform appendToList getVarbayesResp getBic getDataLikelihood pca setDomain spiralgen pixmapToVector reBuild binnedEntropy normalizeVariable semispheregen gramschmidt eigenMppca normMppca appendToGmm newGmm gmmToMppca newMppca subGmm ZtoLabels incremMerge subMppca appendToMppca subVarbayes isNonVoid getQforComp gaussianKL mppcaToGmm l2norm circlegen covgen datagen plotGmm randomGmm

Documented in appendToGmm appendToList appendToMppca binnedEntropy circlegen covgen dat1sample dat2sample dat3sample datagen eigenMppca gaussianKL generate2Dtransform generateSparsePoints getBic getDataLikelihood getQforComp getVarbayesResp gmmToMppca gramschmidt incremMerge isNonVoid l2norm mppcaToGmm newGmm newMppca normalizeVariable normMppca pca pixmapToVector plotGmm randomGmm reBuild semispheregen setDomain spiralgen subGmm subMppca subVarbayes ZtoLabels

# Copyright (C) 2011 Pierrick Bruneau, see README for full notice



randomGmm <- function(domain=10) {
	# generates random 2D GMMs
	mod <- newGmm()
	K <- rpois(1, 5)
	if(K == 0) K <- rpois(1, 5)
	# generate weights from dirichlet
	mod$w <- rDirichlet(K)
	
	# generate means & covs from runif and covgen over some specified domain
	for(i in 1:K) {
		mod$mean[[i]] <- numeric()
		for(j in 1:2) {
			mod$mean[[i]][j] <- runif(1, min=-domain, max=domain)
		}
		
		mod$cov[[i]] <- covgen(bounds=c(2, 10))
	}
	return(mod)
}

plotGmm <- function(mod, steps=200) {
	# get domain from means and covariances
	K <- length(mod$w)
	if(length(mod$mean[[1]]) != 2) stop("only 2D GMMs can be used")
	
	low_left <- c(+Inf, +Inf)
	up_right <- c(-Inf, -Inf)
	
	for(i in 1:K) {
		cur_mean <- mod$mean[[i]]
		cur_vars <- diag(mod$cov[[i]])
		for(d in 1:2) {
			if((cur_mean[d] - 2 * cur_vars[d]) < low_left[d]) low_left[d] <- cur_mean[d] - 2 * cur_vars[d]
			if((cur_mean[d] + 2 * cur_vars[d]) > up_right[d]) up_right[d] <- cur_mean[d] + 2 * cur_vars[d]
		}
	}
	
	stepx <- (up_right[1] - low_left[1]) / steps
	stepy <- (up_right[2] - low_left[2]) / steps
	
	
	pdf <- matrix(nrow=steps+1, ncol=steps+1)
	
	for(i in 0:(steps+1)) {
		for(j in 0:(steps+1)) {
			pdf[i,j] <- gmmdensity(mod, t(as.matrix(c(low_left[1] + i * stepx, low_left[2] + j * stepy))))
		}
	}
	
	lattice::wireframe(pdf, shade = TRUE, 
		aspect = c(1, 1), 
		xlab="", ylab="", zlab="", light.source = c(10,0,10))


}


datagen <- function(dreal=2, deff=6, npts=200, noise=0.1, genmean=rep(0,dreal), genspan=6, iso=FALSE) {
	data <- matrix(0, nrow=npts, ncol=deff)
	
	combs <- rnorm(dreal * (deff - dreal), sd=1)
	combs <- matrix(combs, nrow=deff-dreal, ncol=dreal)
	
	gencov <- matrix(0, nrow=dreal, ncol=dreal)
	# random choice of diagonal elements
	for(i in 1:dreal) {
		gencov[i,i] <- runif(1, min=1, max=genspan)
	}

	noise <- genspan / 60

	# generate random corelations if isotropic is set to F
	if(!iso) {
		for(i in 1:dreal) {
			for(j in i:dreal) {
				if(i != j) {
					# through property found in wikipedia:Positive definiteness
					thres <- (gencov[i,i] + gencov[j,j]) / 2
					gencov[i,j] = runif(1, min=-thres, max=thres)
					gencov[j,i] = gencov[i,j]
				}
			}
		}
	}
	
	print(gencov)
	
	data[,1:dreal] <- mnormt::rmnorm(npts, mean=genmean, varcov=gencov)
	
	if(dreal < deff) {
		for(i in (dreal+1):deff) {
			data[,i] <- rnorm(npts, sd=noise)
			for(j in 1:dreal) {
				data[,i] <- data[,i] + combs[i-dreal, j] * data[,j]
			}
		}
	}
	
	return(list(data=data, genmean=genmean, gencov=gencov))
}





covgen <- function(d=2, bounds=c(1, 5)) {
	gencov <- matrix(0,nrow=d, ncol=d)
	done <- FALSE
	while(!done) {
		for(i in 1:d) {
			gencov[i,i] <- runif(1, min=bounds[1], max=bounds[2])
		}
	
		for(i in 1:d) {
			for(j in i:d) {
				if(i != j) {
					thres <- (gencov[i,i] + gencov[j,j]) / 2
					gencov[i,j] <- gencov[j,i] <- runif(1, min=-thres, max=thres)
				}
			}
		}
		# test for positive definiteness
		vals <- eigen(gencov, only.values=TRUE)$values
		cond <- vals[length(vals)] / vals[1]
		if(cond > 10^(-10)) done <- TRUE
	}
	return(gencov)
}




circlegen <- function(npts=200, radius=10, noise=1) {
	dat <- matrix(0, nrow=npts, ncol=2)
	
	for(i in 1:npts) {
		cur <- runif(1, min=0, max=2*pi)
		dat[i,] <- c(radius*cos(cur), radius*sin(cur)) + mnormt::rmnorm(1, varcov=noise*diag(2))
	}
	
	return(dat)
}

l2norm <- function(vec) {
	val <- 0
	for(i in 1:length(vec)) {
		val <- val + vec[i]^2
	}
	return(sqrt(val))
}

mppcaToGmm <- function(model, notau=FALSE) {
	# convert W nu... format to w mean cov
	output <- list()
	output$w <- numeric()
	output$mean <- list()
	output$cov <- list()
	
	d <- length(as.numeric(model$mumean[[1]]))
	
	output$w <- model$alpha
	output$w <- output$w / sum(output$w)
	for(i in 1:length(model$alpha)) {
		output$mean[[i]] <- model$mumean[[i]]
		if(!notau) {
			output$cov[[i]] <- model$wmean[[i]] %*% t(model$wmean[[i]]) + (1/model$taumoment[i]) * diag(d)
		} else {
			output$cov[[i]] <- model$wmean[[i]] %*% t(model$wmean[[i]]) + (10^(-1)) * diag(d)
		}
	}

	return(output)
}






gaussianKL <- function(N0, N1) {
	d <- length(N0[1,])
	return(log(det(N1)) - log(det(N0)) + sum(diag(solve(N1) %*% N0)) - d)
}

getQforComp <- function(loadings, tau=1.0, verbose=FALSE, quick=FALSE) {
	d <- length(loadings[,1])
	q <- length(loadings[1,])
	
	if(!quick) {
		cur <- loadings[,1] %*% t(loadings[,1]) + (1/tau) * diag(d)
		curind <- 1
		delta <- +Inf
		while((delta > 0.1) && (curind < q)) {
			curind <- curind + 1
			new <- loadings[,1:curind] %*% t(loadings[,1:curind]) + (1/tau) * diag(d)
			delta <- gaussianKL(new, cur)
			cur <- new

		}
		if(delta > 0.1) {
			# last dimension added still has an interest : curind is q.
			curind <- q
			if(verbose) print(paste("all dimensions are used, q is", curind))
		} else {
			if(verbose) print(paste("adding dim.", curind, "makes delta=", delta, ", q set to", curind-1))
			curind <- curind - 1
		}
	} else {
		# quicker solution : analyze column norms
		norms <- numeric()
		for(i in 1:q) {
			norms[i] <- sqrt(sum(loadings[,i] * loadings[,i]))
		}
		print(paste("min:", min(norms), "max:", max(norms)))
		tot <- sum(norms)
		cur <- 0
		curind <- 0
		while((cur/tot)<0.95) {
			curind <- curind + 1
			cur <- cur + norms[curind]
		}
	}
	return(curind)
	
}

isNonVoid <- function(loadings) {
	d <- length(loadings[,1])
	q <- length(loadings[1,])
	result <- FALSE
	for(i in 1:q) {
		for(j in 1:d) {
			if(loadings[j,i] > 0.1) result <- TRUE
		}
	}
	return(result)	
}


subVarbayes <- function(model, thres=2.001) {
	# "sub" a model returned by varbayes
	d <- length(model$data[1,])
	
	filter1 <- list()
	filter1$model <- list()
	filter1$model$alpha <- numeric()
	filter1$model$beta <- numeric()
	filter1$model$nu <- numeric()
	filter1$model$mean <- list()
	filter1$model$wish <- list()
	#filter1$model$resp <- numeric()
	filter1$data <- model$data
	
	selection <- numeric()
	for(i in 1:length(model$model$alpha)) {
		if(model$model$alpha[i] > thres) {
			selection <- c(selection, i)
			filter1$model$alpha <- c(filter1$model$alpha, model$model$alpha[i])
			filter1$model$nu <- c(filter1$model$nu, model$model$nu[i])
			filter1$model$beta <- c(filter1$model$beta, model$model$beta[i])
			filter1$model$mean <- appendToList(filter1$model$mean, model$model$mean[[i]])
			filter1$model$wish <- appendToList(filter1$model$wish, model$model$wish[[i]])
		}
	}
	filter1$model$resp <- model$model$resp[,selection]	
	return(filter1)
}


appendToMppca <- function(mppca1, mppca2) {
	# append mppca2 to mppca1
	# uses appendToList
	mppca1$alpha <- c(mppca1$alpha, mppca2$alpha)
	mppca1$numoment <- appendToList(mppca1$numoment, mppca2$numoment, appendList=TRUE)
	mppca1$nua <- c(mppca1$nua, mppca2$nua)
	mppca1$nub <- appendToList(mppca1$nub, mppca2$nub, appendList=TRUE)
	mppca1$taumoment <- c(mppca1$taumoment, mppca2$taumoment)
	mppca1$taua <- c(mppca1$taua, mppca2$taua)
	mppca1$taub <- c(mppca1$taub, mppca2$taub)
	mppca1$wmean <- appendToList(mppca1$wmean, mppca2$wmean, appendList=TRUE)
	mppca1$wsigma <- appendToList(mppca1$wsigma, mppca2$wsigma, appendList=TRUE)
	mppca1$xsigma <- appendToList(mppca1$xsigma, mppca2$xsigma, appendList=TRUE)
	mppca1$mumean <- appendToList(mppca1$mumean, mppca2$mumean, appendList=TRUE)
	mppca1$musigma <- appendToList(mppca1$musigma, mppca2$musigma, appendList=TRUE)
	mppca1$mustar <- appendToList(mppca1$mustar, mppca2$mustar, appendList=TRUE)
	
	return(mppca1)
}




subMppca <- function(model, prune=FALSE, thres=2.001, quick=FALSE, noxmean=TRUE) {
	# ncol for selecting default number of column factors
	
	d <- length(as.numeric(model$mumean[[1]]))
	q <- length(model$numoment[[1]])

	# select xmean, x1mean or x2mean, depending on cases.

	filter1 <- list()
	filter1$alpha <- numeric()
	filter1$numoment <- list()
	filter1$nua <- numeric()
	filter1$nub <- list()
	filter1$taumoment <- numeric()
	filter1$taua <- numeric()
	filter1$taub <- numeric()
	filter1$wmean <- list()
	filter1$wsigma <- list()
	filter1$xsigma <- list()
	filter1$mumean <- list()
	filter1$musigma <- list()
	filter1$mustar <- list()
	if(!is.null(model$xmean)) {
		if(!noxmean) filter1$xmean <- list()
	} else if(!is.null(model$x1mean)) {
		if(!noxmean) {
			filter1$x1mean <- list()
			filter1$x2mean <- list()
		}
	}



	for(i in 1:length(model$alpha)) {
		if(model$alpha[i] > thres) {
			filter1$alpha <- c(filter1$alpha, model$alpha[i])
			filter1$numoment <- appendToList(filter1$numoment, model$numoment[[i]])
			filter1$nua <- c(filter1$nua, model$nua[i])
			filter1$nub <- appendToList(filter1$nub, model$nub[[i]])
			filter1$taumoment <- c(filter1$taumoment, model$taumoment[i])
			filter1$taua <- c(filter1$taua, model$taua[i])
			filter1$taub <- c(filter1$taub, model$taub[i])
			filter1$wmean <- appendToList(filter1$wmean, model$wmean[[i]])
			filter1$wsigma <- appendToList(filter1$wsigma, model$wsigma[[i]])
			filter1$xsigma <- appendToList(filter1$xsigma, model$xsigma[[i]])
			filter1$mumean <- appendToList(filter1$mumean, model$mumean[[i]])
			filter1$musigma <- appendToList(filter1$musigma, model$musigma[[i]])
			filter1$mustar <- appendToList(filter1$mustar, model$mustar[[i]])
			if(!is.null(filter1$xmean)) {
				filter1$xmean <- appendToList(filter1$xmean, model$xmean[[i]])
			} else if(!is.null(filter1$x1mean)) {
				filter1$x1mean <- appendToList(filter1$x1mean, model$x1mean[[i]])
				filter1$x2mean <- appendToList(filter1$x2mean, model$x2mean[[i]])
			}

		}
	}
	


	qmax <- 1
	

	
	if(prune) {
		for(i in 1:length(filter1$alpha)) {
			if(isNonVoid(filter1$wmean[[i]])) {
				curq <- getQforComp(filter1$wmean[[i]], filter1$taumoment[i], quick=quick)
				if(curq > qmax) qmax <- curq 
			}
		}
	} else {
		qmax <- q
	}
	
	print(paste(qmax, "factors selected"))

	

	filter2 <- list()

	filter2$alpha <- filter1$alpha
	filter2$numoment <- list()
	filter2$nua <- filter2$nua
	filter2$nub <- list()
	filter2$taumoment <- filter1$taumoment
	filter2$taua <- filter1$taua
	filter2$taub <- filter1$taub
	filter2$wmean <- list()
	filter2$wsigma <- list()
	filter2$xsigma <- list()
	filter2$mumean <- filter1$mumean
	filter2$musigma <- filter1$musigma
	filter2$mustar <- filter1$mustar
	if(!is.null(filter1$xmean)) {
		if(!noxmean) filter2$xmean <- list()
	} else if(!is.null(filter1$x1mean)) {
		if(!noxmean) {
			filter2$x1mean <- list()
			filter2$x2mean <- list()
		}
	}


	for(i in 1:length(filter1$alpha)) {
		filter2$numoment[[i]] <- filter1$numoment[[i]][1:qmax]
		filter2$nub[[i]] <- filter1$nub[[i]][1:qmax]
		if(qmax > 1) {
			filter2$wmean[[i]] <- filter1$wmean[[i]][,1:qmax]
			filter2$wsigma[[i]] <- filter1$wsigma[[i]][1:qmax, 1:qmax]
			filter2$xsigma[[i]] <- filter1$xsigma[[i]][1:qmax, 1:qmax]
			if(!is.null(filter2$xmean)) {
				filter2$xmean[[i]] <- filter1$xmean[[i]][,1:qmax]
			} else if(!is.null(filter2$x1mean)) {
				filter2$x1mean[[i]] <- filter1$x1mean[[i]][,1:qmax]
				filter2$x2mean[[i]] <- list()
				for(j in 1:length(filter1$x2mean[[i]])) {
					filter2$x2mean[[i]][[j]] <- filter1$x2mean[[i]][[j]][1:qmax, 1:qmax]
				}
			}
		} else {
			filter2$wmean[[i]] <- as.matrix(filter1$wmean[[i]][,1])
			filter2$wsigma[[i]] <- filter1$wsigma[[i]][1,1]
			filter2$xsigma[[i]] <- filter1$xsigma[[i]][1,1]
			if(!is.null(filter2$xmean)) {
				filter2$xmean[[i]] <- as.matrix(filter1$xmean[[i]][,1])
			} else if(!is.null(filter2$x1mean)) {
				filter2$x1mean[[i]] <- as.matrix(filter1$x1mean[[i]][,1])
				filter2$x2mean[[i]] <- list()
				for(j in 1:length(filter1$x2mean[[i]])) {
					filter2$x2mean[[i]][[j]] <- filter1$x2mean[[i]][[j]][1,1]
				}

			}	

		}
	}

	return(filter2)
	
}



incremMerge <- function(modref, newmod, k=200, nit=100, quick=FALSE) {
	modref$alpha <- modref$alpha * 100000 / sum(modref$alpha)
	newmod$alpha <- newmod$alpha * 100000 / sum(newmod$alpha)
	# inputs are assumed to be correctly filtered
	
	modref <- appendToMppca(modref, newmod)
	modref <- normMppca(modref)

	modref <- mmppca(modref, k, maxit=nit)
	modref <- subMppca(modref)

	return(modref)
}



ZtoLabels <- function(resp) {
	labs <- numeric()
	for(i in 1:length(resp[,1])) {
		labs[i] <- which.max(resp[i,])
	}
	return(labs)
}





subGmm <- function(model, dims=c(1,2), inds=NULL) {
	output <- list()
	output$mean <- list()
	output$cov <- list()

	
	if(is.null(inds)) inds <- 1:length(model$w)
	output$w <- model$w[inds]
	output$w <- output$w / sum(output$w)
	
	#for(i in 1:length(model$w)) {
	for(i in 1:length(inds)) {
		cur <- model$mean[[inds[i]]]
		cur <- as.numeric(cur)
		cur <- cur[dims]
		cur <- as.matrix(cur)
		output$mean[[i]] <- cur
		output$cov[[i]] <- model$cov[[inds[i]]][dims, dims]
	}
	return(output)
}


newMppca <- function() {
	out <- list()
	out$alpha <- numeric()
	out$wmean <- list()
	out$mumean <- list()
	out$numoment <- list()
	out$taumoment <- numeric()
	return(out)
}


gmmToMppca <- function(model, alpha=500) {
	# alpha is the overall population assicated with output MPPCA
	output <- newMppca()
	k <- length(model$w)
	d <- length(as.numeric(model$mean[[1]]))
	
	q <- 1
	eigens <- list()
	# get the dimensionality to retain
	for(i in 1:k) {
		eigens[[i]] <- eigen(model$cov[[i]])
		ratio <- 0
		curq <- 0
		while(ratio < 0.8) {
			curq <- curq + 1
			ratio <- sum(eigens[[i]]$values[1:curq]) / sum(eigens[[i]]$values)
		}
		if(curq > q) q <- curq		
	}
	
	print(paste("q selected :", q))
	
	for(i in 1:k) {
		output$alpha[i] <- model$w[i] * alpha
		output$mumean[[i]] <- model$mean[[i]]
		if(q > 1) {
			output$wmean[[i]] <- eigens[[i]]$vectors[,1:q] %*% diag(eigens[[i]]$values[1:q])
		} else {
			output$wmean[[i]] <- as.matrix(eigens[[i]]$vectors[,1]) * eigens[[i]]$values[1]
		}
		# set columns to same sign
		for(j in 1:q) {
			if(output$wmean[[i]][1,j] < 0) output$wmean[[i]][,j] <- - output$wmean[[i]][,j]
		}
		
		output$numoment[[i]] <- 1 / eigens[[i]]$values[1:q]
		output$taumoment[i] <- 1 / mean(eigens[[i]]$values[(q+1):d])
	}
	
	return(output)
}



newGmm <- function() {
	newmod <- list()
	newmod$w <- numeric()
	newmod$mean <- list()
	newmod$cov <- list()
	newmod$a <- numeric()
	return(newmod)
}

appendToGmm <- function(mod1, mod2) {
	if(length(mod1$a)==0) {
		mod2$a <- rep(0, length(mod2$w))
	} else {
		mod2$a <- rep(mod1$a[length(mod1$a)] + 1, length(mod2$w))
	}
	mod1$a <- c(mod1$a, mod2$a)
	mod1$w <- c(mod1$w, mod2$w)
	mod1$mean <- appendToList(mod1$mean, mod2$mean, appendList=TRUE)
	mod1$cov <- appendToList(mod1$cov, mod2$cov, appendList=TRUE)
	return(mod1)
}



normMppca <- function(mppca1) {
	maxq <- 0
	len <- length(mppca1$alpha)
	for(i in 1:len) {
		curq <- length(mppca1$numoment[[i]])
		if(curq > maxq) maxq <- curq
	}
	for(i in 1:len) {
		curq <- length(mppca1$numoment[[i]])
		d <- length(mppca1$wmean[[i]][,1])
		if(curq < maxq) {
			newwmean <- matrix(0, nrow=d, ncol=maxq)
			newwmean[,1:curq] <- mppca1$wmean[[i]]
			mppca1$wmean[[i]] <- newwmean
			mppca1$numoment[[i]][(curq+1):maxq] <- 100
		}
	}
	print(paste("normalized at", maxq, "factors"))
	return(mppca1)
}



eigenMppca <- function(mod) {
	# diagonalize t(W) %*% W => eigenvectors (cols) are t(R)
	# input is assumed to having been reduced (subMppca)
	
	d <- length(mod$mumean[[1]])
	
	for(i in 1:length(mod$alpha)) {
		# first build the matrix to diagonalize
		mat <- t(mod$wmean[[i]]) %*% mod$wmean[[i]]
		
		# perform the diag
		op <- eigen(mat)
		Rt <- op$vectors
		vals <- op$values
		
		# see tipping PPCA
		mod$wmean[[i]] <- mod$wmean[[i]] %*% Rt
		
	}
	return(mod)
	
}

gramschmidt <- function(mat) {
	# perform gramschmidt normalisation on the set of column vectors in mat
	
	d <- length(mat[1,])
	for(i in 1:d) {
		for(j in 1:i) {
			if(j < i) {
				# previous vectors are normalized : no norm term.
				proj <- mat[,j] * as.numeric(t(mat[,i]) %*% mat[,j])
				mat[,i] <- mat[,i] - proj
			}
		}
		# when finished, normalize current vector
		norm <- as.numeric(t(mat[,i]) %*% mat[,i])
		mat[,i] <- mat[,i] / sqrt(norm)
	}
	return(mat)
}


semispheregen <- function(npts=200, radius=10, noise=1) {
	dat <- matrix(nrow=npts, ncol=3)
	
	for(i in 1:npts) {
		theta <- runif(1, min=0, max=pi/2)
		phi <- runif(1, min=-pi, max=pi)
		dat[i,] <- c(radius*cos(theta)*cos(phi), radius*cos(theta)*sin(phi), radius*sin(theta)) + 
			mnormt::rmnorm(1, varcov=noise*diag(3))
	}
	return(dat)
}







normalizeVariable <- function(v) {
	themin <- min(v)
	themax <- max(v)
	v <- v - themin
	v <- v / (themax - themin)
	return(v)
}


binnedEntropy <- function(v, nbins=100) {
	v <- normalizeVariable(v)
	step <- 1/nbins
	counts <- numeric()
	cumul <- 0
	
	for(i in 1:nbins) {
		cur <- sum(as.numeric(v <= step*i)) - cumul
		counts <- c(counts, cur)
		cumul <- cumul + cur
	}
	counts <- counts / sum(counts)
	
	entropy <- 0
	for(i in 1:nbins) {
		if(counts[i] != 0) {
			entropy <- entropy - counts[i] * log(counts[i]) 
		}
	}
	return(entropy)
}


reBuild <- function(v, voids, nonvoids, domains, placeholder=1) {
	# rescale v to correct domain
	v <- t(as.matrix(as.numeric(v)))
	
	v <- as.numeric(setDomain(v, domains, 10))
	v2 <- numeric()
	for(i in 1:length(v)) {
		v2[nonvoids[i]] <- v[i]
	}
	
	v2[voids] <- placeholder
	v2 <- pixmap::pixmapGrey(matrix(v2, nrow=sqrt(length(v2)), ncol=sqrt(length(v2))))

	return(v2)
	
}



pixmapToVector <- function(p) {
	return(as.numeric(pixmap::getChannels(p, "grey")))
}



spiralgen <- function(radius=10, n=1000, laps=2, noise=1) {
	nptsperlap <- n / laps
	res <- matrix(nrow=n, ncol=2)
	for(i in 1:n) {
		curangle <- ((i %% nptsperlap) / nptsperlap) * (2 * pi)
		curradius <- (1+ (i / nptsperlap)) * (radius)
		res[i,1] <- curradius * cos(curangle) + rnorm(1, sd=noise)
		res[i,2] <- curradius * sin(curangle) + rnorm(1, sd=noise)
	}
	return(res)
	
}








setDomain <- function(dat, span=10, oldspan=NULL) {
	# possible lengths for span :
	# - class = numeric et 1 = all axes on [-span, span]
	# - class = numeric and 2 = all axes on [span[1], span[2]]
	# - class = list et length 2*d = full specification of axes.
	# see "set domain demonstration"
	
	n <- length(dat[,1])
	d <- length(dat[1,])
	
	maxs2 <- numeric()
	mins2 <- numeric()
	maxs1 <- numeric()
	mins1 <- numeric()
	
	# symmetric or asymetric domains ; depends on length of "span" argument
	if((class(span) == "numeric") && (length(span) == 1)) {
		if(span < 0) span <- -span
		maxs2 <- rep(span, d)
		mins2 <- rep(-span, d)
	} else if(class(span) == "numeric" && (length(span) == 2)) {
		if(span[1] > span[2]) {
			temp <- span[1]
			span[1] <- span[2]
			span[2] <- temp
		}
		maxs2 <- rep(span[2], d)
		mins2 <- rep(span[1], d)
	} else if((class(span) == "list") && (length(span) == 2)) {
		for(i in 1:d) {
			if(span[[1]][i] > span[[2]][i]) {
				temp <- span[[1]][i]
				span[[1]][i] <- span[[2]][i]
				span[[2]][i] <- temp
			}
			maxs2[i] <- span[[2]][i]
			mins2[i] <- span[[1]][i]
		}
	} else {
		stop("inadequate 'span' argument")
	}
	
	

	if(is.null(oldspan)) {
		mins1 <- apply(dat, 2, min)
		maxs1 <- apply(dat, 2, max)		
	} else if((class(oldspan) == "numeric") && (length(oldspan) == 1)) {
		if(oldspan < 0) oldspan <- -oldspan
		maxs1 <- rep(oldspan, d)
		mins1 <- rep(-oldspan, d)
	} else if(class(oldspan) == "numeric" && (length(oldspan) == 2)) {
		if(oldspan[1] > oldspan[2]) {
			temp <- oldspan[1]
			oldspan[1] <- oldspan[2]
			oldspan[2] <- temp
		}
		maxs1 <- rep(oldspan[2], d)
		mins1 <- rep(oldspan[1], d)
	} else if((class(oldspan) == "list") && (length(oldspan) == 2)) {
		for(i in 1:d) {
			if(oldspan[[1]][i] > oldspan[[2]][i]) {
				temp <- oldspan[[1]][i]
				oldspan[[1]][i] <- oldspan[[2]][i]
				oldspan[[2]][i] <- temp
			}
			maxs1[i] <- oldspan[[2]][i]
			mins1[i] <- oldspan[[1]][i]
		}
	} else {
		stop("inadequate 'oldspan' argument")
	}
	
	
	
	
	for(i in 1:n) {
		for(j in 1:d) {
			dat[i,j] <- (mins2[j] * (maxs1[j] - dat[i,j]) + maxs2[j] * (dat[i,j] - mins1[j])) / (maxs1[j] - mins1[j]) 
		}
	}
	
	return(dat)
}




pca <- function(dat, ncomp=NULL) {
	# force-cast to matrix
	dat <- as.matrix(dat)

	# if ncomp is NULL, all columns are selected
	if(is.null(ncomp)) ncomp <- dim(dat)[2]
	
	# nombre d'individus = longueur d'une colonne
	n <- length(dat[,1])
	
	# matrice des poids, vecteur de 1
	D <- diag(rep(1/n, n))
	one <- as.matrix(rep(1,n))
	
	# centred data
	Y <- (diag(n) - one %*% t(one) %*% D) %*% dat
	
	# variance
	V <- t(Y) %*% D %*% Y
	
	# D(1/s2) metric
	d1s2 <- diag(1/diag(V))
	
	# diagonalisation
	eig <- eigen(V %*% d1s2)

	eigvecs <- Re(eig$vectors)
	
	# fonction pour M-normer chaque vecteur de notre ensemble de vecteurs propres
	# en effet eigen les I-norme a 1, nous on veut utiliser D(1/s2)
	A <- apply(eigvecs, 2, function(X) { norm <- sqrt(as.numeric(t(X) %*% d1s2 %*% X))
	return(X / norm)})
	
	#info <- "eigenvalues :"
	#for(i in 1:length(eig$values)) {
	#	info <- paste(info, eig$values[i])
	#}
	#
	#print(info)
	eigvals <- Re(eig$values)
	message("ncomp=", ncomp, ", with ", format(100 * sum(eigvals[1:ncomp]) / sum(eigvals), digits=2),
		"% of retained variance")
	
	# U = transposee(A^{-1}), cf cours et TP ACP
	U <- t(solve(A))
	
	# transformed data (XU dans le TP)
	reduce <- Y %*% U
	
	# on retourne les ncomp premieres composantes (2D par defaut)
	return(reduce[,1:ncomp])
}



getDataLikelihood <- function(gmm, dat) {
	densities <- gmmdensity(gmm, dat)
	res <- sum(log(densities))
	return(res)
}


getBic <- function(gmm, dat) {
	lnL <- getDataLikelihood(gmm, dat)
	n <- length(dat[,1])
	d <- length(dat[1,])
	k <- length(gmm$w)
	# constraint on w
	paramsize <- (k-1) + k*d + k*(d*(d+1))/2
	res <- -2 * lnL + paramsize * log(n)
	return(res)
}


getVarbayesResp <- function(data, model) {
	# get responsibilities according to varbayes model
	n <- dim(data)[1]
	d <- dim(data)[2]
	k <- length(model$model$alpha)
	resp <- matrix(nrow=n, ncol=k)
	gammamat <- matrix(nrow=k, ncol=d)
	kappalog <- numeric()
	
	for(j in 1:d) {
		#print(model$model$nu)
		#print(digamma((model$model$nu - j + 1) / 2))
		#print(k)
		gammamat[,j] <- digamma((model$model$nu - j + 1) / 2)
	}
	
	for(j in 1:k) {
		kappalog[j] <- sum(gammamat[j,]) + 2 * d + log(det(model$model$wish[[j]]))
	}
	
	for(i in 1:n) {
		resp[i,] <- digamma(model$model$alpha) - digamma(sum(model$model$alpha))
		resp[i,] <- resp[i,] + 0.5 * kappalog
		for(j in 1:k) {
			resp[i,j] <- resp[i,j] - 0.5 * (d / model$model$beta[j] + model$model$nu[j] * t(as.matrix(data[i,] - model$model$mean[[j]])) %*% model$model$wish[[j]] %*% as.matrix(data[i,] - model$model$mean[[j]]))
		}
		themax <- max(resp[i,])
		resp[i,] <- resp[i,] - themax
		resp[i,] <- exp(resp[i,])
		resp[i,] <- resp[i,] / sum(resp[i,])
	}
	return(resp)
}


appendToList <- function(lst, obj, appendList=FALSE) {
	if(appendList) {
		for(i in 1:length(obj)) {
			lst[[length(lst) + 1]] <- obj[[i]]
		}
	} else {
		lst[[length(lst) + 1]] <- obj
	}
	return(lst)
}



generate2Dtransform <- function(dims=4) {
	transform <- matrix(rnorm(2*dims), nrow=dims, ncol=2)
	# check linear independency <=> X^T %*% X has full rank
	
	while(det(t(transform) %*% transform) < 0.001) {
		transform <- matrix(rnorm(2*dims), nrow=dims, ncol=2)
		message("degenerate 2D transform - trying again...")
	}
	
	transform <- gramschmidt(transform)
	# 2 first dims = original signal
	transform[1:2, 1:2] <- diag(2)
	return(transform)
}



dat1sample <- function(nelts, gmm, noise, transform=generate2Dtransform(2), oldbounds=NULL, newbounds=NULL) {
	# generate elements
	d <- dim(transform)[1]
	dat <- gmmgen(gmm, nelts)[[1]]
	dat <- dat %*% t(transform) + matrix(rnorm(nelts * d, sd=noise), nrow=nelts, ncol=d)
	
	# change domain if necessary
	if(!is.null(newbounds)) {
		if(is.null(oldbounds)) {
			oldbounds <- list()
			oldbounds[[1]] <- apply(dat, 2, min)
			oldbounds[[2]] <- apply(dat, 2, max)
		}
		dat <- setDomain(dat, newbounds, oldbounds)
	}
	return(dat)
}



dat2sample <- function(nelts, radius, noise, oldbounds=NULL, newbounds=NULL) {
	dat <- semispheregen(nelts, radius, noise)
	
	# change domain if necessary
	if(!is.null(newbounds)) {
		if(is.null(oldbounds)) {
			oldbounds <- list()
			oldbounds[[1]] <- apply(dat, 2, min)
			oldbounds[[2]] <- apply(dat, 2, max)
		}
		dat <- setDomain(dat, newbounds, oldbounds)
	}
	return(dat)
	
}

dat3sample <- function(nelts, radius, noise, transform=generate2Dtransform(2), oldbounds=NULL, newbounds=NULL) {
	# noise for transform = circle noise / 10
	d <- dim(transform)[1]
	dat <- circlegen(nelts, radius, noise)
	dat <- dat %*% t(transform) + matrix(rnorm(nelts * d, sd=noise/10), nrow=nelts, ncol=d)	
	
	# change domain if necessary
	if(!is.null(newbounds)) {
		if(is.null(oldbounds)) {
			oldbounds <- list()
			oldbounds[[1]] <- apply(dat, 2, min)
			oldbounds[[2]] <- apply(dat, 2, max)
		}
		dat <- setDomain(dat, newbounds, oldbounds)
	}
	return(dat)
	
}

generateSparsePoints <- function(npoints, dim=2, span=10, mindist=2, maxit=20) {
	converge <- FALSE
	thepoints <- matrix(nrow=npoints, ncol=dim)
	notfar <- 1:npoints
	if(span < 0) span <- -span
	nit <- 1
	while((!converge) && (nit < maxit)) {
		print(paste("it", nit, ",", length(notfar), "points to gen."))	
		thepoints[notfar,] <- runif(length(notfar) * dim, -span, span)
		notfar <- which(.Call("control", thepoints, mindist, PACKAGE="VBmix") == 0)
		converge <- (length(notfar) == 0)
		nit <- nit + 1
	}
	
	if(!converge) print("could not converge, current points are returned")
	return(thepoints)
}

Try the VBmix package in your browser

Any scripts or data that you put into this service are public.

VBmix documentation built on May 30, 2017, 2:34 a.m.