#........................................................
source("D:/VALERIE/ENSEIGNEMENT/INTRO_DM/COURS/em_gauss_mixture.R")
source("D:/VALERIE/ENSEIGNEMENT/INTRO_DM/COURS/kmeans_vm.R")

d <- 2
M <- 3
mu <- matrix(0,d,M)
sigma <- matrix(0,d^2,M)
p1 <- 1/2
p2 <- 1/4
p3 <- 1-(p1+p2)
mu[,1] <- c(5,4)
mu[,2] <- c(-3,2)
mu[,3] <- c(0,-1)

sigma[,1] <- c(2,1,1,1)
sigma[,2] <- c(4,0,0,0.2)
sigma[,3] <- c(2,0.2,.2,.5)

n <- 1000
u <- runif(n)
indic_1 <- u<p1
indic_2 <- u>(p1+p2)
indic_3 <- as.logical(1-(indic_1+indic_2))
n1 <- sum(indic_1)
n2 <- sum(indic_2)
n3 <- sum(indic_3)

x <- matrix(0,n,2)
x[indic_1,] <- rmvnorm(n1,mean=mu[,1],sigma=matrix(sigma[,1],2,2))
x[indic_2,] <- rmvnorm(n2,mean=mu[,2],sigma=matrix(sigma[,2],2,2))
x[indic_3,] <- rmvnorm(n3,mean=mu[,3],sigma=matrix(sigma[,3],2,2))
plot(x[indic_1,1],x[indic_1,2],col="black",xlim=c(min(x[,1]),max(x[,1])),ylim = c(min(x[,2]),max(x[,2])))
points(x[indic_2,1],x[indic_2,2],col="blue")
points(x[indic_3,1],x[indic_3,2],col="green")
title("Points simulés, une couleur par classe")

%%% Pause -> voir graphiques

data <- x
theta0<- NULL
theta0$pi <- c(p1,1-(p1+p3),p3)
theta0$mu <- matrix(c(-5,2,5,2, 0,-1),d,M)
sigma0 <- matrix(0,d^2,M)
sigma0[,1] <- c(1,0,0,1)
sigma0[,2] <- c(1,0,0,1)
sigma0[,3] <- c(1,0,0,1)

theta0$sigma <- sigma0

theta <- theta0
MaxIter <- 100
epsilon <- 1e-3
em_res <- EM_gauss(data,M,theta0,MaxIter,epsilon,1,1) 
em_res$theta$mu

%%% Pause -> voir graphiques

kmeans_res <- kmeans_vm(data,M,theta0,MaxIter,epsilon,1,1)
kmeans_res$theta$mu

%%% Pause -> voir graphiques

t1 <- seq(min(x[,1]),max(x[,1]), length=100)
t2 <- seq(min(x[,2]),max(x[,2]), length=100)
d1 <- expand.grid(x=t1,y=t2)

theta <- em_res$theta
p_pos <- proba_pos(d1,M,theta)
L <- 128
col <- rainbow(L)
col <- col[length(col) :1]
par(mfrow=c(2,M))
k <- 2
ind <- round((p_pos[,k]-min(p_pos[,k]))/(max(p_pos[,k])-min(p_pos[,k]))*(L/2-1))+1 
plot(d1[,1],d1[,2],col=col[ind],pch=19)
points(x[indic_1,1],x[indic_1,2],col="black",pch=20)
title(c("EM, classe",k))
k <- 1
ind <- round((p_pos[,k]-min(p_pos[,k]))/(max(p_pos[,k])-min(p_pos[,k]))*(L/2-1))+1 
plot(d1[,1],d1[,2],col=col[ind],pch=19)
points(x[indic_2,1],x[indic_2,2],col="black",pch=20)
title(c("EM, classe",k))
k <- 3
ind <- round((p_pos[,k]-min(p_pos[,k]))/(max(p_pos[,k])-min(p_pos[,k]))*(L/2-1))+1 
plot(d1[,1],d1[,2],col=col[ind],pch=19)
points(x[indic_3,1],x[indic_3,2],col="black",pch=20)
title(c("EM, classe",k))


theta <- kmeans_res$theta
p_pos <- proba_pos(d1,M,theta)
k <- 2
ind <- round((p_pos[,k]-min(p_pos[,k]))/(max(p_pos[,k])-min(p_pos[,k]))*(L/2-1))+1 
plot(d1[,1],d1[,2],col=col[ind],pch=19)
points(x[indic_1,1],x[indic_1,2],col="black",pch=20)
title(c("Kmeans, classe",k))
k <- 1
ind <- round((p_pos[,k]-min(p_pos[,k]))/(max(p_pos[,k])-min(p_pos[,k]))*(L/2-1))+1 
plot(d1[,1],d1[,2],col=col[ind],pch=19)
points(x[indic_2,1],x[indic_2,2],col="black",pch=20)
title(c("Kmeans, classe",k))
k <- 3
ind <- round((p_pos[,k]-min(p_pos[,k]))/(max(p_pos[,k])-min(p_pos[,k]))*(L/2-1))+1 
plot(d1[,1],d1[,2],col=col[ind],pch=19)
points(x[indic_3,1],x[indic_3,2],col="black",pch=20)
title(c("Kmeans, classe",k))