From f227455a1604906b255ef366d64c10a93e796983 Mon Sep 17 00:00:00 2001
From: emilie <emilie@devijver.org>
Date: Thu, 16 Mar 2017 15:49:30 +0100
Subject: [PATCH] update EMGLLF.R and some few details

---
 .Rbuildignore                        |   2 +
 .Rhistory                            | 512 +++++++++++++++++++++++++++
 .gitignore                           |   1 +
 R/essai16mars.R                      |  88 +++++
 R/generateSampleInputs.R             |  19 +-
 R/gridLambda.R                       |   4 +-
 R/initSmallEM.R                      |   5 +-
 src/test/generate_test_data/EMGLLF.R |  28 +-
 8 files changed, 631 insertions(+), 28 deletions(-)
 create mode 100644 .Rbuildignore
 create mode 100644 .Rhistory
 create mode 100644 R/essai16mars.R

diff --git a/.Rbuildignore b/.Rbuildignore
new file mode 100644
index 0000000..91114bf
--- /dev/null
+++ b/.Rbuildignore
@@ -0,0 +1,2 @@
+^.*\.Rproj$
+^\.Rproj\.user$
diff --git a/.Rhistory b/.Rhistory
new file mode 100644
index 0000000..6616dbe
--- /dev/null
+++ b/.Rhistory
@@ -0,0 +1,512 @@
+lines(1:15, rep(0,15), col='red')
+N       = 300
+mugamma = c(1,2,3)
+Gamma   = matrix(c(5,-0.02,0.01,-0.02,1,0.1,0.01,0.1,1),3,3)
+s2      = 0.01
+Si      = getbasismatrix(seq(0,1,length=15),create.fourier.basis(nbasis = 3))
+gammai  = t(sapply(1:N,FUN = function(i){
+rmvnorm(1,mugamma,Gamma)
+}))
+Xdata   = t(sapply(1:N,FUN = function(i){
+as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
+}))
+matplot(t(Xdata),type="l")
+#beta = c(.1,.5,-.3)
+beta = c(1,-3.5,4.5,-2.5)
+pydata = sapply(1:N,FUN=function(i){
+a = sum(beta*c(1,gammai[i,]))
+return(exp(a)/(1+exp(a)))
+})
+U     = runif(length(pydata))
+ydata = as.numeric(U<pydata)
+plot(sort(pydata),col=ydata[order(pydata)]+2)
+wi = Si %*% beta[-1]
+plot(wi, type='l')
+lines(1:15,rep(0,15))
+plot( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),typ="l",lty=1,col="green")
+for (ite in 1:30){
+print(ite)
+Xdata   = t(sapply(1:N,FUN = function(i){
+as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
+}))
+#beta = c(.1,.5,-.3)
+beta = c(1,-3.5,4.5,-2.5)
+pydata = sapply(1:N,FUN=function(i){
+a = sum(beta*c(1,gammai[i,]))
+return(exp(a)/(1+exp(a)))
+})
+U     = runif(length(pydata))
+ydata = as.numeric(U<pydata)
+lines( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),col="green")
+}
+lines(1:15, rep(0,15), col='red')
+matplot(Xdata)
+matplot(Xdata,type='l')
+matplot(t(Xdata),type='l')
+matplot(apply(Xdata, 2, mean),type='l')
+mean(apply(Xdata, 2, mean))
+N       = 3000
+mugamma = c(1,2,3)
+Gamma   = matrix(c(5,-0.02,0.01,-0.02,1,0.1,0.01,0.1,1),3,3)
+s2      = 0.01
+Si      = getbasismatrix(seq(0,1,length=15),create.fourier.basis(nbasis = 3))
+gammai  = t(sapply(1:N,FUN = function(i){
+rmvnorm(1,mugamma,Gamma)
+}))
+Xdata   = t(sapply(1:N,FUN = function(i){
+as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
+}))
+matplot(t(Xdata),type="l")
+#beta = c(.1,.5,-.3)
+beta = c(1,-3.5,4.5,-2.5)
+pydata = sapply(1:N,FUN=function(i){
+a = sum(beta*c(1,gammai[i,]))
+return(exp(a)/(1+exp(a)))
+})
+U     = runif(length(pydata))
+ydata = as.numeric(U<pydata)
+plot(sort(pydata),col=ydata[order(pydata)]+2)
+wi = Si %*% beta[-1]
+plot(wi, type='l')
+lines(1:15,rep(0,15))
+plot( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),typ="l",lty=1,col="green")
+for (ite in 1:30){
+print(ite)
+Xdata   = t(sapply(1:N,FUN = function(i){
+as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
+}))
+#beta = c(.1,.5,-.3)
+beta = c(1,-3.5,4.5,-2.5)
+pydata = sapply(1:N,FUN=function(i){
+a = sum(beta*c(1,gammai[i,]))
+return(exp(a)/(1+exp(a)))
+})
+U     = runif(length(pydata))
+ydata = as.numeric(U<pydata)
+lines( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),col="green")
+}
+lines(1:15, rep(0,15), col='red')
+Gamma
+Gamma   = matrix(c(5,-0.02,0.01,-0.02,5,0.1,0.01,0.1,5),3,3)
+s2      = 0.01
+Si      = getbasismatrix(seq(0,1,length=15),create.fourier.basis(nbasis = 3))
+gammai  = t(sapply(1:N,FUN = function(i){
+rmvnorm(1,mugamma,Gamma)
+}))
+Xdata   = t(sapply(1:N,FUN = function(i){
+as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
+}))
+matplot(t(Xdata),type="l")
+#beta = c(.1,.5,-.3)
+plot(sort(sapply(1:N,FUN=function(i) 1/(1+exp((-1)*(resEM$betah[1,ITfinal]+sum(resEM$betah[-1,ITfinal]*resEM$gammahat[i,]))))) ))
+plot(wi, type='l')
+wi = Si %*% beta[-1]
+plot(wi, type='l')
+lines(1:15,rep(0,15))
+plot( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),typ="l",lty=1,col="green")
+for (ite in 1:30){
+print(ite)
+Xdata   = t(sapply(1:N,FUN = function(i){
+as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
+}))
+#beta = c(.1,.5,-.3)
+beta = c(1,-3.5,4.5,-2.5)
+pydata = sapply(1:N,FUN=function(i){
+a = sum(beta*c(1,gammai[i,]))
+return(exp(a)/(1+exp(a)))
+})
+U     = runif(length(pydata))
+ydata = as.numeric(U<pydata)
+lines( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),col="green")
+}
+lines(1:15, rep(0,15), col='red')
+N       = 200
+M       = 100
+mugamma = c(1,2,3,2,4,6)
+LGamma  = matrix(0,nrow=6,ncol=6)
+LGamma[lower.tri(LGamma,diag=FALSE)] = rnorm(6*5/2,0,0.05)
+LGamma  = LGamma + diag(runif(6,0.2,1))
+Gamma   = LGamma%*%t(LGamma)
+s2      = 0.2
+freq      = seq(0,1,length=M)
+intervals = list(12:25,67:89)
+Si        = as.matrix(bdiag(getbasismatrix(freq[12:25],create.fourier.basis(nbasis = 3)), getbasismatrix(freq[67:89],create.fourier.basis(nbasis = 3))))
+gammai  = t(sapply(1:N,FUN = function(i){
+rmvnorm(1,mugamma,Gamma)
+}))
+Xdata   = t(sapply(1:N,FUN = function(i){
+as.vector(Si%*%gammai[i,])+rnorm(nrow(Si),0,s2)
+}))
+xdata = matrix(0,nrow=N,ncol=M)
+xdata[,intervals[[1]] ] = Xdata[,1:14]
+xdata[,intervals[[2]] ] = Xdata[,15:37]
+beta = c(2,0.5,-2,1,0.2,1,-1)
+pydata = sapply(1:N,FUN=function(i){
+a = sum(beta*c(1,gammai[i,]))
+return(exp(a)/(1+exp(a)))
+})
+U     = runif(length(pydata))
+ydata = as.numeric(U<pydata)
+plot(sort(pydata),col=ydata[order(pydata)]+1)
+nBasis    = c(3,3)
+L   = 20
+eer = rep(0,L)
+for (ell in 1:L){
+kapp      = sample(1:200,0.8*200,replace=FALSE)
+xdata.app = xdata[kapp,]
+ydata.app = ydata[kapp]
+ktest      = setdiff(1:200,kapp)
+xdata.test = xdata[ktest,]
+ydata.test = ydata[ktest]
+resEM    = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=1000,25,eps=10^-8,keep=TRUE)
+Yres     = Y.predB(xdata.test,freq,intervals,resEM$muh,resEM$Gammah,resEM$s2h,resEM$betah,nBasis=c(3,3))
+eer[ell] = mean(Yres$Ypred==ydata.test)
+cat("\n", "ell =", ell)
+}
+source('~/Dropbox/projet spectres/algoEM/EMalgorithmBandes.R')
+for (ell in 1:L){
+kapp      = sample(1:200,0.8*200,replace=FALSE)
+xdata.app = xdata[kapp,]
+ydata.app = ydata[kapp]
+ktest      = setdiff(1:200,kapp)
+xdata.test = xdata[ktest,]
+ydata.test = ydata[ktest]
+resEM    = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=1000,25,eps=10^-8,keep=TRUE)
+Yres     = Y.predB(xdata.test,freq,intervals,resEM$muh,resEM$Gammah,resEM$s2h,resEM$betah,nBasis=c(3,3))
+eer[ell] = mean(Yres$Ypred==ydata.test)
+cat("\n", "ell =", ell)
+}
+N       = 3000
+mugamma = c(1,2,3)
+Gamma   = matrix(c(5,-0.02,0.01,-0.02,5,0.1,0.01,0.1,5),3,3)
+s2      = 0.01
+Si      = getbasismatrix(seq(0,1,length=15),create.fourier.basis(nbasis = 3))
+gammai  = t(sapply(1:N,FUN = function(i){
+rmvnorm(1,mugamma,Gamma)
+}))
+Xdata   = t(sapply(1:N,FUN = function(i){
+as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
+}))
+matplot(t(Xdata),type="l")
+#beta = c(.1,.5,-.3)
+beta = c(1,-3.5,4.5,-2.5)
+pydata = sapply(1:N,FUN=function(i){
+a = sum(beta*c(1,gammai[i,]))
+return(exp(a)/(1+exp(a)))
+})
+U     = runif(length(pydata))
+ydata = as.numeric(U<pydata)
+plot(sort(pydata),col=ydata[order(pydata)]+2)
+wi = Si %*% beta[-1]
+plot(wi, type='l')
+lines(1:15,rep(0,15))
+plot( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),typ="l",lty=1,col="green")
+for (ite in 1:30){
+print(ite)
+Xdata   = t(sapply(1:N,FUN = function(i){
+as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
+}))
+#beta = c(.1,.5,-.3)
+beta = c(1,-3.5,4.5,-2.5)
+pydata = sapply(1:N,FUN=function(i){
+a = sum(beta*c(1,gammai[i,]))
+return(exp(a)/(1+exp(a)))
+})
+U     = runif(length(pydata))
+ydata = as.numeric(U<pydata)
+lines( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),col="green")
+}
+lines(1:15, rep(0,15), col='red')
+N       = 3000
+mugamma = c(1,2,3)
+Gamma   = matrix(c(5,-0.02,0.01,-0.02,1,0.1,0.01,0.1,5),3,3)
+s2      = 0.01
+Si      = getbasismatrix(seq(0,1,length=15),create.fourier.basis(nbasis = 3))
+gammai  = t(sapply(1:N,FUN = function(i){
+rmvnorm(1,mugamma,Gamma)
+}))
+Xdata   = t(sapply(1:N,FUN = function(i){
+as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
+}))
+matplot(t(Xdata),type="l")
+#beta = c(.1,.5,-.3)
+beta = c(1,-3.5,4.5,-2.5)
+pydata = sapply(1:N,FUN=function(i){
+a = sum(beta*c(1,gammai[i,]))
+return(exp(a)/(1+exp(a)))
+})
+U     = runif(length(pydata))
+ydata = as.numeric(U<pydata)
+plot(sort(pydata),col=ydata[order(pydata)]+2)
+# ydata = ydata.noise
+wi = Si %*% beta[-1]
+plot(wi, type='l')
+lines(1:15,rep(0,15))
+plot( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),typ="l",lty=1,col="green")
+for (ite in 1:30){
+print(ite)
+Xdata   = t(sapply(1:N,FUN = function(i){
+as.vector(Si%*%gammai[i,])+rnorm(15,0,s2)
+}))
+#beta = c(.1,.5,-.3)
+beta = c(1,-3.5,4.5,-2.5)
+pydata = sapply(1:N,FUN=function(i){
+a = sum(beta*c(1,gammai[i,]))
+return(exp(a)/(1+exp(a)))
+})
+U     = runif(length(pydata))
+ydata = as.numeric(U<pydata)
+lines( apply(t(Xdata[ydata==1,]),1,mean)-apply(t(Xdata[ydata==0,]),1,mean),col="green")
+}
+lines(1:15, rep(0,15), col='red')
+matplot(t(Xdata[ydata==1,]),typ="l",lty=1,col="red")
+matplot(t(Xdata[ydata==0,]),typ="l",lty=1,col="black",add=TRUE)
+N       = 200
+M       = 100
+mugamma = c(1,2,3,2,4,6)
+LGamma  = matrix(0,nrow=6,ncol=6)
+LGamma[lower.tri(LGamma,diag=FALSE)] = rnorm(6*5/2,0,0.05)
+LGamma  = LGamma + diag(runif(6,0.2,1))
+Gamma   = LGamma%*%t(LGamma)
+s2      = 0.2
+freq      = seq(0,1,length=M)
+intervals = list(12:25,67:89)
+Si        = as.matrix(bdiag(getbasismatrix(freq[12:25],create.fourier.basis(nbasis = 3)), getbasismatrix(freq[67:89],create.fourier.basis(nbasis = 3))))
+gammai  = t(sapply(1:N,FUN = function(i){
+rmvnorm(1,mugamma,Gamma)
+}))
+Xdata   = t(sapply(1:N,FUN = function(i){
+as.vector(Si%*%gammai[i,])+rnorm(nrow(Si),0,s2)
+}))
+xdata = matrix(0,nrow=N,ncol=M)
+xdata[,intervals[[1]] ] = Xdata[,1:14]
+xdata[,intervals[[2]] ] = Xdata[,15:37]
+beta = c(2,0.5,-2,1,0.2,1,-1)
+pydata = sapply(1:N,FUN=function(i){
+a = sum(beta*c(1,gammai[i,]))
+return(exp(a)/(1+exp(a)))
+})
+U     = runif(length(pydata))
+ydata = as.numeric(U<pydata)
+plot(sort(pydata),col=ydata[order(pydata)]+1)
+nBasis    = c(3,3)
+resEM = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=500,100,eps=10^-8,keep=TRUE)
+source('~/Dropbox/projet spectres/algoEM/EMalgorithmBandes.R')
+resEM = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=500,100,eps=10^-8,keep=TRUE)
+L   = 20
+eer = rep(0,L)
+for (ell in 1:L){
+kapp      = sample(1:200,0.8*200,replace=FALSE)
+xdata.app = xdata[kapp,]
+ydata.app = ydata[kapp]
+ktest      = setdiff(1:200,kapp)
+xdata.test = xdata[ktest,]
+ydata.test = ydata[ktest]
+resEM    = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=1000,25,eps=10^-8,keep=FALSE)
+Yres     = Y.predB(xdata.test,freq,intervals,resEM$muh,resEM$Gammah,resEM$s2h,resEM$betah,nBasis=c(3,3))
+eer[ell] = mean(Yres$Ypred==ydata.test)
+cat("\n", "ell =", ell)
+}
+source('~/Dropbox/projet spectres/algoEM/EMalgorithmBandes.R')
+for (ell in 1:L){
+kapp      = sample(1:200,0.8*200,replace=FALSE)
+xdata.app = xdata[kapp,]
+ydata.app = ydata[kapp]
+ktest      = setdiff(1:200,kapp)
+xdata.test = xdata[ktest,]
+ydata.test = ydata[ktest]
+resEM    = EM.Bands.function(ydata,xdata,freq,intervals,nBasis,MCit=1000,25,eps=10^-8,keep=FALSE)
+Yres     = Y.predB(xdata.test,freq,intervals,resEM$muh,resEM$Gammah,resEM$s2h,resEM$betah,nBasis=c(3,3))
+eer[ell] = mean(Yres$Ypred==ydata.test)
+cat("\n", "ell =", ell)
+}
+eer
+mean(eer)
+source('~/Dropbox/projet spectres/cluster/EMLiqArtBandeHugues.R')
+library(fda)
+library(Matrix)
+library(mvtnorm)
+source('EMalgorithm.R')
+source('EMLiqArtBandeHugues.R')
+for (nVC in 1:100){
+EMLiqArtBandeHugues(nVC)
+}
+setwd("~/Dropbox/projet spectres/cluster")
+library(fda)
+library(Matrix)
+library(mvtnorm)
+source('EMalgorithm.R')
+source('EMLiqArtBandeHugues.R')
+for (nVC in 1:100){
+EMLiqArtBandeHugues(nVC)
+}
+sims= c(63:67,69,71:74,76:85,88:90,93:100)
+#simulation_settings=expand.grid(c(1,2,3,4,5),c(1,2),c(600),c(300))
+#for (i in 1:nrow(simulation_settings)){
+aa=cbind(sims,rep(9, each =length(sims)))
+colnames(aa)=c("isim","model")
+write.table(aa,file=paste("SimulationSettings9.txt",sep=""),row.names=FALSE,sep=",")
+setwd("~/Dropbox/WarpMixedModel/warpingMixedModel/SuperComputer")
+write.table(aa,file=paste("SimulationSettings9.txt",sep=""),row.names=FALSE,sep=",")
+sims= 1:100
+colnames(sims) = "isim"
+sims= 1:100
+#simulation_settings=expand.grid(c(1,2,3,4,5),c(1,2),c(600),c(300))
+#for (i in 1:nrow(simulation_settings)){
+aa=cbind(sims,rep(10, each =length(sims)))
+colnames(aa)=c("isim","model")
+write.table(aa,file=paste("SimulationSettings10.txt",sep=""),row.names=FALSE,sep=",")
+load("~/valse/data/data.RData")
+load("~/valse/data/data.RData")
+X
+Y
+library("devtools", lib.loc="~/R/x86_64-pc-linux-gnu-library/3.3")
+library(roxygen2)
+setwd("~/valse")
+document()
+document()
+setwd("~/")
+install("valse")
+c(9.75,
+20,
+11.75,
+11,
+8.25,
+12.5,
+11,
+18,
+7,
+12.75,
+13,
+14.75,
+4.75,
+11,
+20,
+13.5,
+8,
+13.25,
+6,
+17.5,
+13.25,
+8.5,
+9.5,
+16,
+8,
+9.5,
+10.25,
+13.75,
+9,
+14.75,
+12.5,
+19.5,
+17.5,
+7,
+11.5,
+4,
+7.5,
+13.25,
+10.5
+)
+notes = c(9.75,
+20,
+11.75,
+11,
+8.25,
+12.5,
+11,
+18,
+7,
+12.75,
+13,
+14.75,
+4.75,
+11,
+20,
+13.5,
+8,
+13.25,
+6,
+17.5,
+13.25,
+8.5,
+9.5,
+16,
+8,
+9.5,
+10.25,
+13.75,
+9,
+14.75,
+12.5,
+19.5,
+17.5,
+7,
+11.5,
+4,
+7.5,
+13.25,
+10.5
+)
+hist(notes)
+hist(notes, nclass = 20)
+notesBis = c(7.375,19.375,10.125,9,6,12.25,11,18.75,10.5,11.5,15.5,13,2.375,12.75,17.75,15.375,8.125,16,8.5,14.5,11.625,11.25,9.625,14,
+6.5,11.375,11.875,16.25,10.125,12,6.25,9.75,8.75,3.5,0,0,0,5.75,2,3.75,6.625,5.25)
+hist(notesBis)
+notes = c(7,
+19,
+10,
+9,
+6,
+12,
+11,
+19,
+11,
+12,
+16,
+13,
+2,
+13,
+18,
+15,
+8,
+16,
+9,
+15,
+12,
+11,
+10,
+14,
+7,
+11,
+12,
+16,
+13,
+15,
+14,
+19,
+14,
+4,  9,
+8,
+7,
+8,
+8
+)
+hist(notes)
+shiny::runApp('Dropbox/cranview-master')
+packages = 'shock'
+min_date = Sys.Date() - 1
+for (pkg in packages)
+{
+# api data for package. we want the initial release - the first element of the "timeline"
+pkg_data = httr::GET(paste0("http://crandb.r-pkg.org/", pkg, "/all"))
+pkg_data = httr::content(pkg_data)
+initial_release = pkg_data$timeline[[1]]
+min_date = min(min_date, as.Date(initial_release))
+}
+min_date
+load("~/valse/data/data.RData")
+k = 2
+init = initSmallEM(k, X, Y)
+setwd("~/valse")
+library("devtools", lib.loc="~/R/x86_64-pc-linux-gnu-library/3.3")
+library("roxygen2", lib.loc="~/R/x86_64-pc-linux-gnu-library/3.3")
diff --git a/.gitignore b/.gitignore
index 812257f..b8ee8f9 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,2 +1,3 @@
 # Ignore Rstudio project files
 *.Rproj*
+.Rproj.user
diff --git a/R/essai16mars.R b/R/essai16mars.R
new file mode 100644
index 0000000..d835599
--- /dev/null
+++ b/R/essai16mars.R
@@ -0,0 +1,88 @@
+meanX = rep(0,6)
+covX = 0.1*diag(6)
+
+covY = array(dim = c(5,5,2))
+covY[,,1] = 0.1*diag(5)
+covY[,,2] = 0.2*diag(5)
+
+beta = array(dim = c(6,5,2))
+beta[,,2] = matrix(c(rep(2,12),rep(0, 18)))
+beta[,,1] = matrix(c(rep(1,12),rep(0, 18)))
+
+n = 500
+
+pi = c(0.4,0.6)
+
+source('~/valse/R/generateSampleInputs.R')
+data = generateXY(meanX,covX,covY, pi, beta, n)
+
+X = data$X
+Y = data$Y
+
+k = 2
+
+n = nrow(Y)
+m = ncol(Y)
+p = ncol(X)
+
+Zinit1 = array(0, dim=c(n))
+betaInit1 = array(0, dim=c(p,m,k))
+sigmaInit1 = array(0, dim = c(m,m,k))
+phiInit1 = array(0, dim = c(p,m,k))
+rhoInit1 = array(0, dim = c(m,m,k))
+Gam = matrix(0, n, k)
+piInit1 = matrix(0,k)
+gamInit1 = array(0, dim=c(n,k))
+LLFinit1 = list()
+
+require(MASS) #Moore-Penrose generalized inverse of matrix
+
+  distance_clus = dist(X)
+  tree_hier = hclust(distance_clus)
+  Zinit1 = cutree(tree_hier, k)
+  sum(Zinit1==1)
+  
+  for(r in 1:k)
+  {
+    Z = Zinit1
+    Z_indice = seq_len(n)[Z == r] #renvoit les indices où Z==r
+    if (length(Z_indice) == 1) {
+      betaInit1[,,r] = ginv(crossprod(t(X[Z_indice,]))) %*%
+        crossprod(t(X[Z_indice,]), Y[Z_indice,])
+    } else {
+      betaInit1[,,r] = ginv(crossprod(X[Z_indice,])) %*%
+        crossprod(X[Z_indice,], Y[Z_indice,])
+    }
+    sigmaInit1[,,r] = diag(m)
+    phiInit1[,,r] = betaInit1[,,r] #/ sigmaInit1[,,r]
+    rhoInit1[,,r] = solve(sigmaInit1[,,r])
+    piInit1[r] = mean(Z == r)
+  }
+  
+  for(i in 1:n)
+  {
+    for(r in 1:k)
+    {
+      dotProduct = tcrossprod(Y[i,]%*%rhoInit1[,,r]-X[i,]%*%phiInit1[,,r])
+      Gam[i,r] = piInit1[r]*det(rhoInit1[,,r])*exp(-0.5*dotProduct)
+    }
+    sumGamI = sum(Gam[i,])
+    gamInit1[i,]= Gam[i,] / sumGamI
+  }
+  
+  miniInit = 10
+  maxiInit = 101
+  
+  new_EMG = EMGLLF(phiInit1,rhoInit1,piInit1,gamInit1,miniInit,maxiInit,1,0,X,Y,1e-6)
+  
+  new_EMG$phi
+  new_EMG$pi
+  LLFEessai = new_EMG$LLF
+  LLFinit1 = LLFEessai[length(LLFEessai)]
+
+
+b = which.max(LLFinit1)
+phiInit = phiInit1[,,,b]
+rhoInit = rhoInit1[,,,b]
+piInit = piInit1[b,]
+gamInit = gamInit1[,,b]
diff --git a/R/generateSampleInputs.R b/R/generateSampleInputs.R
index fb67b08..8edd031 100644
--- a/R/generateSampleInputs.R
+++ b/R/generateSampleInputs.R
@@ -1,9 +1,9 @@
 #' Generate a sample of (X,Y) of size n
-#' @param meanX matrix of group means for covariates (of size p*K)
-#' @param covX covariance for covariates (of size p*p*K)
+#' @param meanX matrix of group means for covariates (of size p)
+#' @param covX covariance for covariates (of size p*p)
 #' @param covY covariance for the response vector (of size m*m*K)
 #' @param pi	 proportion for each cluster
-#' @param beta regression matrix
+#' @param beta regression matrix, of size p*m*k
 #' @param n		sample size
 #'
 #' @return list with X and Y
@@ -12,20 +12,23 @@ generateXY = function(meanX, covX, covY, pi, beta, n)
 {
 	p = dim(covX)[1]
 	m = dim(covY)[1]
-	k = dim(covX)[3]
+	k = dim(covY)[3]
 
 	X = matrix(nrow=n,ncol=p)
 	Y = matrix(nrow=n,ncol=m)
+	class = matrix(nrow = n)
 
 	require(MASS) #simulate from a multivariate normal distribution
 	for (i in 1:n)
 	{
-		class = sample(1:k, 1, prob=pi)
-		X[i,] = mvrnorm(1, meanX[,class], covX[,,class])
-		Y[i,] = mvrnorm(1, X[i,] %*% beta[,,class], covY[,,class])
+		class[i] = sample(1:k, 1, prob=pi)
+		X[i,] = mvrnorm(1, meanX, covX)
+		print(X[i,])
+		print(beta[,,class[i]])
+		Y[i,] = mvrnorm(1, X[i,] %*% beta[,,class[i]], covY[,,class[i]])
 	}
 
-	return (list(X=X,Y=Y))
+	return (list(X=X,Y=Y, class = class))
 }
 
 #' Generate a sample of (X,Y) of size n with default values
diff --git a/R/gridLambda.R b/R/gridLambda.R
index e7946ae..35c412a 100644
--- a/R/gridLambda.R
+++ b/R/gridLambda.R
@@ -19,8 +19,8 @@ gridLambda = function(phiInit, rhoInit, piInit, gamInit, X, Y, gamma, mini, maxi
 	m = dim(phiInit)[2]
 	k = dim(phiInit)[3]
 
-	list_EMG = .Call("EMGLLF_core",phiInit,rhoInit,piInit,gamInit,mini,maxi,1,0,X,Y,tau)
-
+	#list_EMG = .Call("EMGLLF_core",phiInit,rhoInit,piInit,gamInit,mini,maxi,1,0,X,Y,tau)
+  list_EMG = EMGLLF(phiInit,rhoInit,piInit,gamInit,mini,maxi,1,0,X,Y,tau)
 	grid = array(0, dim=c(p,m,k))
 	for (i in 1:p)
 	{
diff --git a/R/initSmallEM.R b/R/initSmallEM.R
index 399f39f..541d7e1 100644
--- a/R/initSmallEM.R
+++ b/R/initSmallEM.R
@@ -62,8 +62,9 @@ initSmallEM = function(k,X,Y)
 		miniInit = 10
 		maxiInit = 11
 		
-		new_EMG = .Call("EMGLLF_core",phiInit1[,,,repet],rhoInit1[,,,repet],piInit1[repet,],
-			gamInit1[,,repet],miniInit,maxiInit,1,0,X,Y,1e-4)
+		#new_EMG = .Call("EMGLLF_core",phiInit1[,,,repet],rhoInit1[,,,repet],piInit1[repet,],
+#			gamInit1[,,repet],miniInit,maxiInit,1,0,X,Y,1e-4)
+		new_EMG = EMGLLF(phiInit1[,,,repet],rhoInit1[,,,repet],piInit1[repet,],gamInit1[,,repet],miniInit,maxiInit,1,0,X,Y,1e-4)
 		LLFEessai = new_EMG$LLF
 		LLFinit1[repet] = LLFEessai[length(LLFEessai)]
 	}
diff --git a/src/test/generate_test_data/EMGLLF.R b/src/test/generate_test_data/EMGLLF.R
index 7100f29..272eb6f 100644
--- a/src/test/generate_test_data/EMGLLF.R
+++ b/src/test/generate_test_data/EMGLLF.R
@@ -27,7 +27,6 @@ EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau)
   ps = matrix(0, m,k)
   nY2 = matrix(0, m,k)
   ps1 = array(0, dim=c(n,m,k))
-  nY21 = array(0, dim=c(n,m,k))
   Gam = matrix(0, n,k)
   EPS = 1E-15
   
@@ -58,8 +57,8 @@ EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau)
     ##########
     
     #pour pi
-    for (r in 1:k)
-      b[r] = sum(abs(phi[,,r]))
+    for (r in 1:k){
+      b[r] = sum(abs(phi[,,r]))}
     gam2 = colSums(gam)
     a = sum(gam %*% log(pi))
     
@@ -92,12 +91,9 @@ EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau)
         for (i in 1:n)
 				{
           ps1[i,mm,r] = Y2[i,mm,r] * sum(X2[i,,r] * phi[,mm,r])
-          nY21[i,mm,r] = Y2[i,mm,r]^2
         }
         ps[mm,r] = sum(ps1[,mm,r])
-        nY2[mm,r] = sum(nY21[,mm,r])
-
-#TODO: debug rho computation
+        nY2[mm,r] = sum(Y2[,mm,r]^2)
         rho[mm,mm,r] = (ps[mm,r]+sqrt(ps[mm,r]^2+4*nY2[mm,r]*gam2[r])) / (2*nY2[mm,r])
       }
     }
@@ -107,9 +103,9 @@ EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau)
 			{
         for (mm in 1:m)
 				{
-          S[j,mm,r] = -rho[mm,mm,r]*ps2[j,mm,r] +
-						(if(j>1) sum(phi[1:(j-1),mm,r] * Gram2[j,1:(j-1),r]) else 0) +
-						(if(j<p) sum(phi[(j+1):p,mm,r] * Gram2[j,(j+1):p,r]) else 0)
+          S[j,mm,r] = -rho[mm,mm,r]*ps2[j,mm,r] + sum(phi[-j,mm,r] * Gram2[j, setdiff(1:p, j),r])
+#						(if(j>1) sum(phi[1:(j-1),mm,r] * Gram2[j,1:(j-1),r]) else 0) +
+#						(if(j<p) sum(phi[(j+1):p,mm,r] * Gram2[j,(j+1):p,r]) else 0)
           if (abs(S[j,mm,r]) <= n*lambda*(pi[r]^gamma))
             phi[j,mm,r]=0
           else if(S[j,mm,r] > n*lambda*(pi[r]^gamma))
@@ -128,9 +124,8 @@ EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau)
 		{
       #precompute sq norms to numerically adjust their values
       sqNorm2 = rep(0,k)
-      for (r in 1:k)
-        sqNorm2[r] = sum( (Y[i,]%*%rho[,,r]-X[i,]%*%phi[,,r])^2 )
-      shift = 0.5*min(sqNorm2)
+      for (r in 1:k){
+        sqNorm2[r] = sum( (Y[i,]%*%rho[,,r]-X[i,]%*%phi[,,r])^2 )}
 
       #compute Gam(:,:) using shift determined above
       sumLLF1 = 0.0;
@@ -138,7 +133,7 @@ EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau)
 			{
 				#FIXME: numerical problems, because 0 < det(Rho[,,r] < EPS; what to do ?!
         #       consequence: error in while() at line 77
-				Gam[i,r] = pi[r] * exp(-0.5*sqNorm2[r] + shift) #* det(rho[,,r])
+				Gam[i,r] = pi[r] * exp(-0.5*sqNorm2[r])* det(rho[,,r])
         sumLLF1 = sumLLF1 + Gam[i,r] / (2*base::pi)^(m/2)
       }
       sumLogLLF2 = sumLogLLF2 + log(sumLLF1)
@@ -161,6 +156,7 @@ EMGLLF = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau)
 
     ite = ite+1
   }
-
-  return(list("phi"=phi, "rho"=rho, "pi"=pi, "LLF"=LLF, "S"=S))
+  
+  affec = apply(gam, 1,which.max)
+  return(list("phi"=phi, "rho"=rho, "pi"=pi, "LLF"=LLF, "S"=S, "affec" = affec ))
 }
-- 
2.44.0