From 35b42a4bd37b162a3d579693b2b5fa4913a52ed5 Mon Sep 17 00:00:00 2001 From: Ben_G Date: Wed, 16 Nov 2016 14:33:53 +0100 Subject: [PATCH] Inputparametters.R + SelecModel --- InputParameters/basicInitParameters.R | 29 ++++++++++ InputParameters/generateIO.R | 25 ++++++++ InputParameters/generateIOdefault.R | 25 ++++++++ InputParameters/gridLambda.R | 29 ++++++++++ InputParameters/initSmallEM.R | 80 ++++++++++++++++++++++++++ SelectModel/selectionindice.R | 24 ++++++++ SelectModel/selectionmodele.R | 38 ++++++++++++ SelectModel/suppressionmodelesegaux.R | 18 ++++++ SelectModel/suppressionmodelesegaux2.R | 26 +++++++++ 9 files changed, 294 insertions(+) create mode 100644 InputParameters/basicInitParameters.R create mode 100644 InputParameters/generateIO.R create mode 100644 InputParameters/generateIOdefault.R create mode 100644 InputParameters/gridLambda.R create mode 100644 InputParameters/initSmallEM.R create mode 100644 SelectModel/selectionindice.R create mode 100644 SelectModel/selectionmodele.R create mode 100644 SelectModel/suppressionmodelesegaux.R create mode 100644 SelectModel/suppressionmodelesegaux2.R diff --git a/InputParameters/basicInitParameters.R b/InputParameters/basicInitParameters.R new file mode 100644 index 0000000..b70edab --- /dev/null +++ b/InputParameters/basicInitParameters.R @@ -0,0 +1,29 @@ +basic_Init_Parameters = function(n,p,m,k){ + phiInit = array(0, dim=c(p,m,k)) + + piInit = (1.0/k)*rep.int(1,k) + + rhoInit = array(0, dim=c(m,m,k)) + + for(i in 1:k){ + rhoInit[,,i] = diag(m) + } + + gamInit = 0.1*array(1, dim=c(n,k)) + + R = sample(1:k,n, replace= TRUE) + + for(i in 1:n){ + gamInit[i,R[i]] = 0.9 + } + gamInit = gamInit/sum(gamInit[1,]) + + + return(list(phiInit, rhoInit, piInit, gamInit)) +} + +n= 10 +p = 10 +m = 5 +k = 5 +list_param = basic_Init_Parameters(n,p,m,k) diff --git a/InputParameters/generateIO.R b/InputParameters/generateIO.R new file mode 100644 index 0000000..9e84af5 --- /dev/null +++ b/InputParameters/generateIO.R @@ -0,0 +1,25 @@ +library(MASS) #simulate from a multivariate normal distribution + +generateIO = function(meanX, covX, covY, pi, beta, n){ #don't need meanX + size_covX = dim(covX) + p = size_covX[1] + k = size_covX[3] + + size_covY = dim(covY) + m = size_covY[1] + + Y = matrix(0,n,m) + BX = array(0, dim=c(n,m,k)) + + for(i in 1:n){ + for(r in 1:k){ + BXir = rep(0,m) + for(mm in 1:m){ + Bxir[[mm]] = X[i,] %*% beta[,mm,r] + } + Y[i,]=Y[i,] + pi[[r]] * mvrnorm(1,BXir, covY[,,r]) + } + } + + return(list(X,Y)) +} \ No newline at end of file diff --git a/InputParameters/generateIOdefault.R b/InputParameters/generateIOdefault.R new file mode 100644 index 0000000..1d5c160 --- /dev/null +++ b/InputParameters/generateIOdefault.R @@ -0,0 +1,25 @@ +generateIOdefault = function(n, p, m, k){ + rangeX = 100 + meanX = rangeX*(1-matrix(runif(k*p),ncol = p)) + + covX = array(0, dim=c(p,p,k)) + covY = array(0, dim=c(m,m,k)) + + for(r in 1:k){ + covX[,,r] = diag(p) + covY[,,r] = diag(m) + } + + pi = (1/k) * rep(1,k) + + beta = array(0, dim=c(p,m,k)) + + for(j in 1:p){ + nonZeroCount = ceiling(m * runif(1)) + beta[j,1:nonZeroCount,] = matrix(runif(nonZeroCount*k),ncol = k) + } + + generate = generateIO(meanX, covX, covY, pi, beta, n) + + return(list(generate[[1]],generate[[2]])) +} \ No newline at end of file diff --git a/InputParameters/gridLambda.R b/InputParameters/gridLambda.R new file mode 100644 index 0000000..87e9299 --- /dev/null +++ b/InputParameters/gridLambda.R @@ -0,0 +1,29 @@ +gridLambda = function(phiInit, rhoInit, piInit, gamInit, X, Y, gamma, mini, maxi, tau){ + n = nrow(X) + p = dimension(phiInit)[1] + m = dimension(phiInit)[2] + k = dimension(phiInit)[3] + list_EMG = EMGLLF(phiInit,rhoInit,piInit,gamInit,mini,maxi,1,0,X,Y,tau) + #.C("EMGLLF", phiInit = phiInit, rhoInit = rhoInit, ...) + phi = list_EMG[[1]] + rho = list_EMG[[2]] + pi = list_EMG[[3]] + S = list_EMG[[5]] + + grid = array(0, dim=c(p,m,k)) + for(i in 1:p){ + for(j in 1:m){ + grid[i,j,] = abs(S[i,j,]) / (n*pi^gamma) + } + } + grid = unique(grid) + grid = grid[grid <=1 ] + + return(grid) +} + + +#test pour voir si formatage à la fin de grid ok +grid= array(mvrnorm(5*5*2,1,1), dim=c(5,5,2)) +grid = unique(grid) +grid = grid[grid<= 1 ] diff --git a/InputParameters/initSmallEM.R b/InputParameters/initSmallEM.R new file mode 100644 index 0000000..8f3c86b --- /dev/null +++ b/InputParameters/initSmallEM.R @@ -0,0 +1,80 @@ +library(MASS) #generalized inverse of matrix Monroe-Penrose + +vec_bin = function(X,r){ + Z = c() + indice = c() + j=1 + for(i in 1:length(X)){ + if(X[i] == r){ + Z[i] = 1 + indice[j] = i + j=j+1 + } + else{ + Z[i] = 0 + } + } + return(list(Z,indice)) +} + +initSmallEM = function(k,X,Y,tau){ + n = nrow(Y) + m = ncol(Y) + p = ncol(X) + + betaInit1 = array(0, dim=c(p,m,k,20)) + sigmaInit1 = array(0, dim = c(m,m,k,20)) + phiInit1 = array(0, dim = c(p,m,k,20)) + rhoInit1 = array(0, dim = c(m,m,k,20)) + piInit1 = matrix(0,20,k) + gamInit1 = array(0, dim=c(n,k,20)) + LLFinit1 = list() + + + for(repet in 1:20){ + clusters = hclust(dist(y)) #default distance : euclidean + clusterCut = cutree(clusters,k) + Zinit1[,repet] = clusterCut #retourne les indices (à quel cluster indiv_i appartient) d'un clustering hierarchique (nb de cluster = k) + + for(r in 1:k){ + Z = Zinit1[,repet] + Z_bin = vec_bin(Z,r) + Z_vec = Z_bin[[1]] #vecteur 0 et 1 aux endroits où Z==r + Z_indice = Z_bin[[2]] #renvoit les indices où Z==r + + betaInit1[,,r,repet] = ginv(t(x[Z_indice,])%*%x[Z_indice,])%*%t(x[Z_indice,])%*%y[Z_indice,] + sigmaInit1[,,r,repet] = diag(m) + phiInit1[,,r,repet] = betaInit1[,,r,repet]/sigmaInit1[,,r,repet] + rhoInit1[,,r,repet] = solve(sigmaInit1[,,r,repet]) + piInit1[repet,r] = sum(Z_vec)/n + } + + for(i in 1:n){ + for(r in 1:k){ + dotProduct = (y[i,]%*%rhoInit1[,,r,repet]-x[i,]%*%phiInit1[,,r,repet]) %*% (y[i,]%*%rhoInit1[,,r,repet]-x[i,]%*%phiInit1[,,r,repet]) + Gam[i,r] = piInit1[repet,r]*det(rhoInit1[,,r,repet])*exp(-0.5*dotProduct) + } + sumGamI = sum(gam[i,]) + gamInit1[i,,repet]= Gam[i,] / sumGamI + } + + miniInit = 10 + maxiInit = 11 + + new_EMG = EMGLLF(phiInit1[,,,repet],rhoInit1[,,,repet],piInit1[repet,],gamInit1[,,repet],miniInit,maxiInit,1,0,x,y,tau) + ##.C("EMGLLF", phiInit = phiInit, rhoInit = rhoInit, ...) + LLFEessai = new_EMG[[4]] + LLFinit1[[repet]] = LLFEessai[[length(LLFEessai)]] + } + + b = which.max(LLFinit1) + + phiInit = phiInit1[,,,b] + rhoInit = rhoInit1[,,,b] + piInit = piInit1[b,] + gamInit = gamInit1[,,b] + + return(list(phiInit, rhoInit, piInit, gamInit)) +} + + diff --git a/SelectModel/selectionindice.R b/SelectModel/selectionindice.R new file mode 100644 index 0000000..6782b2e --- /dev/null +++ b/SelectModel/selectionindice.R @@ -0,0 +1,24 @@ +selectionindice = function(phi, seuil){ + dim_phi = dim(phi) + pp = dim_phi[1] + m = dim_phi[2] + + A = matrix(0, pp, m) + B = matrix(0, pp, m) + + for(j in 1:pp){ + cpt1 = 0 + cpt2 = 0 + for(mm in 1:m){ + if(max(phi[j,mm,])> seuil){ + cpt1 = cpt1 + 1 + A[j,cpt] = mm + } + else{ + cpt2 = cpt2+1 + B[j, cpt2] = mm + } + } + } + return(list(A,B)) +} \ No newline at end of file diff --git a/SelectModel/selectionmodele.R b/SelectModel/selectionmodele.R new file mode 100644 index 0000000..b92c4e4 --- /dev/null +++ b/SelectModel/selectionmodele.R @@ -0,0 +1,38 @@ +vec_bin = function(X,r){ + Z = c() + indice = c() + j=1 + for(i in 1:length(X)){ + if(X[i] == r){ + Z[i] = 1 + indice[j] = i + j=j+1 + } + else{ + Z[i] = 0 + } + } + return(list(Z,indice)) +} + +selectionmodele = function(vraisemblance){ + D = vraimsemblance[,2] + D1 = unique(D) + + indice = rep(1, length(D1)) + + #select argmax MLE + if(length(D1)>2){ + for(i in 1:length(D1)){ + A = c() + for(j in 1:length(D)){ + if(D[[j]]==D1[[i]]){ + a = c(a, vraimsemblance[j,1]) + } + } + b = max(a) + indice[i] = which.max(vec_bin(vraimsemblance,b)[[1]]) #retourne le premier indice du vecteur binaire où u_i ==1 + } + } + return(list(indice,D1)) +} \ No newline at end of file diff --git a/SelectModel/suppressionmodelesegaux.R b/SelectModel/suppressionmodelesegaux.R new file mode 100644 index 0000000..a558efb --- /dev/null +++ b/SelectModel/suppressionmodelesegaux.R @@ -0,0 +1,18 @@ +suppressionmodelesegaux = function(B1,B2,glambda,rho,pi){ + ind = c() + for(j in 1:length(glambda)){ + for(ll in 1:(l-1)){ + if(B1[,,l] == B1[,,ll]){ + ind = c(ind, l) + } + } + } + ind = unique(ind) + B1 = B1[,,-ind] + glambda = glambda[-ind] + B2 = B2[,,-ind] + rho = rho[,,,-ind] + pi = pi[,-ind] + + return(list(B1,B2,glambda,ind,rho,pi)) +} \ No newline at end of file diff --git a/SelectModel/suppressionmodelesegaux2.R b/SelectModel/suppressionmodelesegaux2.R new file mode 100644 index 0000000..4574afc --- /dev/null +++ b/SelectModel/suppressionmodelesegaux2.R @@ -0,0 +1,26 @@ +suppressionmodelesegaux2 = function(B1,rho,pi){ + ind = c() + dim_B1 = dim(B1) + B2 = array(0,dim=c(dim_B1[1],dim_B1[2],dim_B1[3])) + nombreLambda=dim_B1[[2]] + glambda = rep(0,nombreLambda) + + #for(j in 1:nombreLambda){ + # for(ll in 1:(l-1)){ + # if(B1[,,l] == B1[,,ll]){ + # ind = c(ind, l) + # } + # } + #} + #ind = unique(ind) + #B1 = B1[,,-ind] + #rho = rho[,,,-ind] + #pi = pi[,-ind] + + suppressmodel = suppressionmodelesegaux(B1,B2,glambda,rho,pi) + B1 = suppressmodel[[1]] + ind = suppressmodel[[4]] + rho = suppressmodel[[5]] + pi = suppressmodel[[6]] + return(list(B1,ind,rho,pi)) +} \ No newline at end of file -- 2.44.0