From: emilie Date: Thu, 13 Apr 2017 14:24:26 +0000 (+0200) Subject: fix initialization and made some update X-Git-Url: https://git.auder.net/%7B%7B%20asset%28%27mixstore/css/user/current/%7B%7B?a=commitdiff_plain;h=4d9db27f0d1749e5577038dedbc5f4d0826f2772;p=valse.git fix initialization and made some update --- diff --git a/pkg/R/computeGridLambda.R b/pkg/R/computeGridLambda.R index 8ec4d66..9d06aed 100644 --- a/pkg/R/computeGridLambda.R +++ b/pkg/R/computeGridLambda.R @@ -33,5 +33,5 @@ computeGridLambda = function(phiInit, rhoInit, piInit, gamInit, X, Y, grid[i,j,] = abs(list_EMG$S[i,j,]) / (n*list_EMG$pi^gamma) } grid = unique(grid) - grid + sort(grid) } diff --git a/pkg/R/initSmallEM.R b/pkg/R/initSmallEM.R index 5dcafb8..bd5ce17 100644 --- a/pkg/R/initSmallEM.R +++ b/pkg/R/initSmallEM.R @@ -13,21 +13,21 @@ initSmallEM = function(k,X,Y, fast=TRUE) n = nrow(Y) m = ncol(Y) p = ncol(X) - - Zinit1 = array(0, dim=c(n,20)) - 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)) + nIte = 20 + Zinit1 = array(0, dim=c(n,nIte)) + betaInit1 = array(0, dim=c(p,m,k,nIte)) + sigmaInit1 = array(0, dim = c(m,m,k,nIte)) + phiInit1 = array(0, dim = c(p,m,k,nIte)) + rhoInit1 = array(0, dim = c(m,m,k,nIte)) Gam = matrix(0, n, k) - piInit1 = matrix(0,20,k) - gamInit1 = array(0, dim=c(n,k,20)) + piInit1 = matrix(0,nIte,k) + gamInit1 = array(0, dim=c(n,k,nIte)) LLFinit1 = list() #require(MASS) #Moore-Penrose generalized inverse of matrix - for(repet in 1:20) + for(repet in 1:nIte) { - distance_clus = dist(X) + distance_clus = dist(cbind(X,Y)) tree_hier = hclust(distance_clus) Zinit1[,repet] = cutree(tree_hier, k) @@ -62,13 +62,12 @@ initSmallEM = function(k,X,Y, fast=TRUE) miniInit = 10 maxiInit = 11 - new_EMG = EMGLLF(phiInit1[,,,repet], rhoInit1[,,,repet], piInit1[repet,], + init_EMG = EMGLLF(phiInit1[,,,repet], rhoInit1[,,,repet], piInit1[repet,], gamInit1[,,repet], miniInit, maxiInit, gamma=1, lambda=0, X, Y, eps=1e-4, fast) - LLFEessai = new_EMG$LLF + LLFEessai = init_EMG$LLF LLFinit1[repet] = LLFEessai[length(LLFEessai)] } - - b = which.max(LLFinit1) + b = which.min(LLFinit1) phiInit = phiInit1[,,,b] rhoInit = rhoInit1[,,,b] piInit = piInit1[b,] diff --git a/pkg/R/main.R b/pkg/R/main.R index 634c273..238160c 100644 --- a/pkg/R/main.R +++ b/pkg/R/main.R @@ -123,36 +123,39 @@ valse = function(X, Y, procedure='LassoMLE', selecMod='DDSE', gamma=1, mini=10, print(tableauRecap) tableauRecap = tableauRecap[which(tableauRecap[,4]!= Inf),] - modSel = capushe::capushe(tableauRecap, n) - indModSel <- - if (selecMod == 'DDSE') - as.numeric(modSel@DDSE@model) - else if (selecMod == 'Djump') - as.numeric(modSel@Djump@model) - else if (selecMod == 'BIC') - modSel@BIC_capushe$model - else if (selecMod == 'AIC') - modSel@AIC_capushe$model - mod = as.character(tableauRecap[indModSel,1]) - listMod = as.integer(unlist(strsplit(mod, "[.]"))) - modelSel = models_list[[listMod[1]]][[listMod[2]]] + return(tableauRecap) - ##Affectations - Gam = matrix(0, ncol = length(modelSel$pi), nrow = n) - for (i in 1:n){ - for (r in 1:length(modelSel$pi)){ - sqNorm2 = sum( (Y[i,]%*%modelSel$rho[,,r]-X[i,]%*%modelSel$phi[,,r])^2 ) - Gam[i,r] = modelSel$pi[r] * exp(-0.5*sqNorm2)* det(modelSel$rho[,,r]) - } - } - Gam = Gam/rowSums(Gam) - modelSel$affec = apply(Gam, 1,which.max) - modelSel$proba = Gam - - if (plot){ - print(plot_valse(X,Y,modelSel,n)) - } - - return(modelSel) + # modSel = capushe::capushe(tableauRecap, n) + # indModSel <- + # if (selecMod == 'DDSE') + # as.numeric(modSel@DDSE@model) + # else if (selecMod == 'Djump') + # as.numeric(modSel@Djump@model) + # else if (selecMod == 'BIC') + # modSel@BIC_capushe$model + # else if (selecMod == 'AIC') + # modSel@AIC_capushe$model + # + # mod = as.character(tableauRecap[indModSel,1]) + # listMod = as.integer(unlist(strsplit(mod, "[.]"))) + # modelSel = models_list[[listMod[1]]][[listMod[2]]] + # + # ##Affectations + # Gam = matrix(0, ncol = length(modelSel$pi), nrow = n) + # for (i in 1:n){ + # for (r in 1:length(modelSel$pi)){ + # sqNorm2 = sum( (Y[i,]%*%modelSel$rho[,,r]-X[i,]%*%modelSel$phi[,,r])^2 ) + # Gam[i,r] = modelSel$pi[r] * exp(-0.5*sqNorm2)* det(modelSel$rho[,,r]) + # } + # } + # Gam = Gam/rowSums(Gam) + # modelSel$affec = apply(Gam, 1,which.max) + # modelSel$proba = Gam + # + # if (plot){ + # print(plot_valse(X,Y,modelSel,n)) + # } + # + # return(modelSel) }