From: emilie <emilie@devijver.org> Date: Thu, 16 Mar 2017 14:49:30 +0000 (+0100) Subject: update EMGLLF.R and some few details X-Git-Url: https://git.auder.net/variants/current/doc/scripts/common.css?a=commitdiff_plain;h=f227455a1604906b255ef366d64c10a93e796983;p=valse.git update EMGLLF.R and some few details --- 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 )) }