X-Git-Url: https://git.auder.net/?a=blobdiff_plain;f=R%2FinitSmallEM.R;h=c24fca91162425ba01480f727c0bda026872c9e0;hb=23fff6d03a3b19a2ee5f70b71808d39ae9b4ef8a;hp=5044a38743f564c71cd8f2fb66536791412bbdd7;hpb=b4899af94061fca34163bdef9ae3ff2155038bb7;p=valse.git diff --git a/R/initSmallEM.R b/R/initSmallEM.R index 5044a38..c24fca9 100644 --- a/R/initSmallEM.R +++ b/R/initSmallEM.R @@ -12,11 +12,13 @@ initSmallEM = function(k,X,Y,tau) n = nrow(Y) m = ncol(Y) p = ncol(X) - + + Zinit1 = array(0, dim=c(n,20)) #doute sur la taille 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)) + Gam = matrix(0, n, k) piInit1 = matrix(0,20,k) gamInit1 = array(0, dim=c(n,k,20)) LLFinit1 = list() @@ -25,18 +27,17 @@ initSmallEM = function(k,X,Y,tau) require(mclust) # K-means with selection of K for(repet in 1:20) { - clusters = Mclust(matrix(c(X,Y),nrow=n),k) #default distance : euclidean + clusters = Mclust(X,k) #default distance : euclidean #Mclust(matrix(c(X,Y)),k) Zinit1[,repet] = clusters$classification for(r in 1:k) { Z = Zinit1[,repet] Z_bin = vec_bin(Z,r) - Z_vec = Z_bin$Z #vecteur 0 et 1 aux endroits o? Z==r + Z_vec = Z_bin$vec #vecteur 0 et 1 aux endroits o? Z==r Z_indice = Z_bin$indice #renvoit les indices o? Z==r - betaInit1[,,r,repet] = - ginv(t(x[Z_indice,])%*%x[Z_indice,])%*%t(x[Z_indice,])%*%y[Z_indice,] + 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]) @@ -47,19 +48,18 @@ initSmallEM = function(k,X,Y,tau) { 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]) + dotProduct = 3 * (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,]) + sumGamI = sum(Gam[i,]) gamInit1[i,,repet]= Gam[i,] / sumGamI } miniInit = 10 maxiInit = 11 - new_EMG = .Call("EMGLLF",phiInit1[,,,repet],rhoInit1[,,,repet],piInit1[repet,], - gamInit1[,,repet],miniInit,maxiInit,1,0,x,y,tau) + new_EMG = .Call("EMGLLF_core",phiInit1[,,,repet],rhoInit1[,,,repet],piInit1[repet,], + gamInit1[,,repet],miniInit,maxiInit,1,0,X,Y,tau) LLFEessai = new_EMG$LLF LLFinit1[repet] = LLFEessai[length(LLFEessai)] }