Inputparametters.R + SelecModel

diff --git a/InputParameters/basicInitParameters.R b/InputParameters/basicInitParameters.R
new file mode 100644 (file)
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 (file)
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 (file)
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 (file)
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 (file)
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 (file)
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 (file)
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 (file)
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 (file)
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