??!

q1.r

library("testit")
library("MASS")
library("mvtnorm")

em.ll <- function(x, means, covariances, mix.props) {
	probabilities <- mapply(dmvnorm, means, covariances, MoreArgs=list(x=x), SIMPLIFY=F)
	mixed_probs <- Map("*", probabilities, mix.props)
	mixed_probs <- Reduce("+", mixed_probs)
	Reduce("+", Map(log, mixed_probs))
}



em.sigmas <- function(x, mus, gs, n.k) {
	normed <- lapply(mus, function(mu, x) {x - mu}, x)
	sigmas <- Map(function(norm, g){

		sigma <- Map(function(x){
			t(x) %*% t(t(x))
		}, split(norm, rownames(norm)))
		
		sigma <- Map("*", sigma, g)
		sigma <- Reduce("+", sigma)

	}, normed, gs)

	sigmas <- Map("/",sigmas, n.k)
	
	assert('need to be symmetric', isSymmetric(sigmas[[1]]))
	return(sigmas)
}

em.mus <- function(x, gs, n.k){
	Map("/", lapply(lapply(gs, "*", x), colSums), n.k)
}

em.responsibilities <- function(mixed_probs){
	ret <- lapply(mixed_probs, "/", Reduce("+", mixed_probs))
	print(ret)
	ret
}

em.mixedNormal <- function(x, means, covariances, pi_i) {
	probabilities <- mapply(dmvnorm, means, covariances, MoreArgs=list(x=x), SIMPLIFY=F)
	Map("*", probabilities, pi_i)
}	

em.norm <- function(x,means,covariances,mix.prop, max_iterations = 150){
	old_LL <- -Inf
	i <- 0
	new_LL <- em.ll(x, means, covariances, mix.prop)
	ll <- list()
	while(i < max_iterations){
		
		mixed_probs <- em.mixedNormal(x, means, covariances, mix.prop)
		resp <- em.responsibilities(mixed_probs)

		n.k <- lapply(resp, sum)
		props <- lapply(n.k, "/", nrow(x))
		mu <- em.mus(x, resp, n.k)
		cov <- em.sigmas(x, means, resp, n.k)
		
		i = i + 1
		new_LL <- em.ll(x, mu, cov, props)
		print(i)
		print(new_LL)

		ll[[i]] <- new_LL

		if(old_LL == new_LL){
			break
		} else {
			old_LL <- new_LL
			means <- mu
			covariances <- cov
			mix.prop <- props
		}
	}

	return(list(means=means, covariances=covariances, mix.props=mix.prop, LL=old_LL, p=ll))
}

#prepare data

get.means <- function(x,K){ 
	split(kmeans(x, K)$centers, seq_len(K))
}

get.covariances <- function(x, K) {
	replicate(K, cov(x), simplify = F)
}



x <- synth.te[1:10,-3]
sigma <- cov(x)
prop <- list(0.01, 0.79, 0.2)
mus <- get.means(3)
covs <- get.covariances(3)


run.K <- function(x, k){
	mix.props <- lapply(rep(1,k),"/", k)
	means <- get.means(x, k)
	covariances <- get.covariances(x,k)
	print(mix.props)
	em.norm(x, means, covariances, mix.props)
}