From 31444abc970b7fe17463bcc916e95846272158db Mon Sep 17 00:00:00 2001 From: emilie <emilie@devijver.org> Date: Fri, 17 Mar 2017 15:07:20 +0100 Subject: [PATCH] few last things --- pkg/R/selectiontotale.R | 20 +- pkg/R/valse.R | 9 +- reports/.Rhistory | 512 +++++++++++++++++++++++++++ reports/bazar Emilie/simulatedData.R | 10 +- reports/essai17mars.R | 0 reports/simulData_17mars.R | 60 ++++ 6 files changed, 592 insertions(+), 19 deletions(-) create mode 100644 reports/.Rhistory create mode 100644 reports/essai17mars.R create mode 100644 reports/simulData_17mars.R diff --git a/pkg/R/selectiontotale.R b/pkg/R/selectiontotale.R index 7209fed..2cdac38 100644 --- a/pkg/R/selectiontotale.R +++ b/pkg/R/selectiontotale.R @@ -29,7 +29,7 @@ selectiontotale = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambd selectedVariables = list() Rho = list() Pi = list() - cpt = 0 + cpt = 1 #Pour chaque lambda de la grille, on calcule les coefficients for (lambdaIndex in 1:length(glambda)){ print(lambdaIndex) @@ -39,14 +39,16 @@ selectiontotale = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambd m = dim(phiInit)[2] #selectedVariables: list where element j contains vector of selected variables in [1,m] if (sum(params$phi) != 0){ - cpt = cpt+1 - selectedVariables[[cpt]] = sapply(1:p, function(j) { - #from boolean matrix mxk of selected variables obtain the corresponding boolean m-vector, - #and finally return the corresponding indices - c(seq_len(m)[ apply( abs(params$phi[j,,]) > thresh, 1, any ) ], rep(0, m-length(seq_len(m)[ apply( abs(params$phi[j,,]) > thresh, 1, any ) ] ) )) - }) - Rho[[cpt]] = params$rho - Pi[[cpt]] = params$pi + selectedVariables[[cpt]] = sapply(1:p, function(j) { + #from boolean matrix mxk of selected variables obtain the corresponding boolean m-vector, + #and finally return the corresponding indices + c(seq_len(m)[ apply( abs(params$phi[j,,]) > thresh, 1, any ) ], rep(0, m-length(seq_len(m)[ apply( abs(params$phi[j,,]) > thresh, 1, any ) ] ) )) + }) + if (length(unique(selectedVariables)) == length(selectedVariables)){ + Rho[[cpt]] = params$rho + Pi[[cpt]] = params$pi + cpt = cpt+1 + } } } list("selected"=selectedVariables,"Rho"=Rho,"Pi"=Pi) diff --git a/pkg/R/valse.R b/pkg/R/valse.R index 72e2d4d..9bc1f4a 100644 --- a/pkg/R/valse.R +++ b/pkg/R/valse.R @@ -16,7 +16,7 @@ #' @export #----------------------------------------------------------------------- valse = function(X,Y,procedure = 'LassoMLE',selecMod = 'DDSE',gamma = 1,mini = 10, - maxi = 50,eps = 1e-4,kmin = 2,kmax = 3, + maxi = 50,eps = 1e-4,kmin = 2,kmax = 2, rang.min = 1,rang.max = 10) { ################################## #core workflow: compute all models @@ -45,9 +45,9 @@ valse = function(X,Y,procedure = 'LassoMLE',selecMod = 'DDSE',gamma = 1,mini = 1 source('~/valse/pkg/R/gridLambda.R') grid_lambda <<- gridLambda(phiInit, rhoInit, piInit, gamInit, X, Y, gamma, mini, maxi, eps) - # if (length(grid_lambda)>50){ - # grid_lambda = grid_lambda[seq(1, length(grid_lambda), length.out = 50)] - # } + if (length(grid_lambda)>100){ + grid_lambda = grid_lambda[seq(1, length(grid_lambda), length.out = 100)] + } print("Compute relevant parameters") #select variables according to each regularization parameter #from the grid: A1 corresponding to selected variables, and @@ -56,7 +56,6 @@ valse = function(X,Y,procedure = 'LassoMLE',selecMod = 'DDSE',gamma = 1,mini = 1 params = selectiontotale(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,grid_lambda,X,Y,1e-8,eps) #params2 = selectVariables(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,grid_lambda[seq(1,length(grid_lambda), by=3)],X,Y,1e-8,eps) ## etrange : params et params 2 sont différents ... - selected <<- params$selected Rho <<- params$Rho Pi <<- params$Pi diff --git a/reports/.Rhistory b/reports/.Rhistory new file mode 100644 index 0000000..dc970cc --- /dev/null +++ b/reports/.Rhistory @@ -0,0 +1,512 @@ +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 +library("capushe", lib.loc="~/R/x86_64-pc-linux-gnu-library/3.3") +A = matrix(rnorm(25), ncol = 5) +A = rbind(A, matrix(0, ncol = 5, nrow = 2)) +A +A[rowSums(A)==0,]=c() +A[rowSums(A)==0,]=[] +A = A[rowSums(A)!=0,] +A +source('~/valse/pkg/R/valse.R') +setwd("~/valse/reports") +XY = cbind(X,Y) +X +p = 10 +q = 8 +k = 2 +D = 20 +meanX = rep(0,p) +covX = 0.1*diag(p) +covY = array(dim = c(q,q,k)) +covY[,,1] = 0.1*diag(q) +covY[,,2] = 0.2*diag(q) +beta = array(dim = c(p,q,2)) +beta[,,2] = matrix(c(rep(2,(D)),rep(0, p*q-D))) +beta[,,1] = matrix(c(rep(1,D),rep(0, p*q-D))) +n = 100 +pi = c(0.4,0.6) +data = generateXY(meanX,covX,covY, pi, beta, n) +source('~/valse/pkg/R/generateSampleInputs.R') +data = generateXY(meanX,covX,covY, pi, beta, n) +X = data$X +Y = data$Y +XY = cbind(X,Y) +data +data$class +affec = data$class +XY = cbind(X,Y) +for (r in 1:K){ +XY_class[[r]] = XY[affec == r, ] +} +K= 2 +for (r in 1:K){ +XY_class[[r]] = XY[affec == r, ] +} +XY_class= list() +for (r in 1:K){ +XY_class[[r]] = XY[affec == r, ] +} +for (r in 1:K){ +XY_class[[r]] = XY[affec == r, ] +meanPerClass[,r] = apply(XY_class[[r]], 2, mean) +} +meanPerClass= matrix() +for (r in 1:K){ +XY_class[[r]] = XY[affec == r, ] +meanPerClass[,r] = apply(XY_class[[r]], 2, mean) +} +apply(XY_class[[r]], 2, mean) +p +q +meanPerClass[r,] = apply(XY_class[[r]], 2, mean) +dim(XY) +meanPerClass= matrix(0, ncol = K, nrow = dim(XY)[2]) +for (r in 1:K){ +XY_class[[r]] = XY[affec == r, ] +meanPerClass[,r] = apply(XY_class[[r]], 2, mean) +} +matplot(meanPerClass) +matplot(meanPerClass, type='l') diff --git a/reports/bazar Emilie/simulatedData.R b/reports/bazar Emilie/simulatedData.R index 7eaefb9..6031822 100644 --- a/reports/bazar Emilie/simulatedData.R +++ b/reports/bazar Emilie/simulatedData.R @@ -10,12 +10,12 @@ library("mclust") T = seq(0,1.5,length.out = p) T2 = seq(0,3, length.out = 2*p) n = 100 - x1 = cos(2*pi*T) + 0.2*cos(4*2*pi*T) +2*c(rep(0,round(length(T)/7)),rep(1,round(length(T)*(1-1/7)))) + x1 = cos(2*base::pi*T) + 0.2*cos(4*2*base::pi*T) + 0.3*c(rep(0,round(length(T)/7)),rep(1,round(length(T)*(1-1/7))))+1 plot(T,x1) lines(T,x1) - - sigmaX = 0.085 - sigmaY = 0.1 + + sigmaX = 0.12 + sigmaY = 0.12 beta = list() p1= 0.5 beta[[1]] =diag(c(rep(p1,5),rep(1,5), rep(p1,5), rep(1, p-15))) @@ -94,7 +94,7 @@ for (ite in c(1:ITE)){ ########### ## k-means 1 ########### - mod1 = Mclust(t(XY[ite,,]),G = 2, mode='VII') + mod1 = Mclust(t(XY[ite,,]),G = 1:2, mode='VII') ARI1[ite] = adjustedRandIndex(mod1$classification, affec[[ite]]) Kmod1[ite] = mod1$G # ########### diff --git a/reports/essai17mars.R b/reports/essai17mars.R new file mode 100644 index 0000000..e69de29 diff --git a/reports/simulData_17mars.R b/reports/simulData_17mars.R new file mode 100644 index 0000000..cdb476c --- /dev/null +++ b/reports/simulData_17mars.R @@ -0,0 +1,60 @@ +simulData_17mars = function(ite){ + set.seed = 22021989+ite + + ########### + ## Modele + ########### + K = 2 + p = 48 + T = seq(0,1.5,length.out = p) + T2 = seq(0,3, length.out = 2*p) + n = 100 + x1 = cos(2*base::pi*T) + 0.2*cos(4*2*base::pi*T) + 0.3*c(rep(0,round(length(T)/7)),rep(1,round(length(T)*(1-1/7))))+1 + sigmaX = 0.12 + sigmaY = 0.12 + beta = list() + p1= 0.5 + beta[[1]] =diag(c(rep(p1,5),rep(1,5), rep(p1,5), rep(1, p-15))) + p2 = 1 + beta[[2]] = diag(c(rep(p2,5),rep(1,5), rep(p2,5), rep(1, p-15))) + ARI1 = ARI2 = ARI3 = 0 + + ########### + ## Data + Projection + ########### + require(wavelets) + XY = array(0, dim = c(2*p,n)) + Xproj = array(0, dim=c(48,n)) + Yproj = array(0, dim=c(48,n)) + XYproj = array(0, dim=c(96,n)) + x = x1 + matrix(rnorm(n*p, 0, sigmaX), ncol = n) + affec = sample(c(1,2), n, replace = TRUE) + y = x + xy = matrix(0,ncol=n, nrow= 2*p) + for (i in c(1:n)){ + y[,i] = x[,i] %*% beta[[affec[i]]] + rnorm(p, 0, sigmaY) + xy[,i] = c(x[,i],y[,i]) + XY[,i] = xy[,i] - mean(xy[,i]) + Dx = dwt(x[,i], filter='haar')@W + Dx = rev(unlist(Dx)) + Dx = Dx[2:(1+3+6+12+24)] + Ax = dwt(x[,i], filter='haar')@V + Ax = rev(unlist(Ax)) + Ax = Ax[2:(1+3)] + Dy = dwt(y[,i], filter='haar')@W + Dy = rev(unlist(Dy)) + Dy = Dy[2:(1+3+6+12+24)] + Ay = dwt(y[,i], filter='haar')@V + Ay = rev(unlist(Ay)) + Ay = Ay[2:(1+3)] + Xproj[,i] = c(Ax,Dx) + Yproj[,i] = c(Ay,Dy) + XYproj[,i] = c(Ax,Dx,Ay,Dy) + } + + res_valse = valse(x,y) + res_valse_proj = valse(Xproj, Yproj) + + save(res_valse,file=paste("./Out/Res_",ite, ".RData",sep="")) + save(res_valse_proj,file=paste("./Out/ResProj_",ite, ".RData",sep="")) +} -- 2.44.0