--- /dev/null
+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)
--- /dev/null
+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
--- /dev/null
+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
--- /dev/null
+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 ]
--- /dev/null
+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))
+}
+
+
--- /dev/null
+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
--- /dev/null
+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
--- /dev/null
+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
--- /dev/null
+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
projects / valse.git / commitdiff