X-Git-Url: https://git.auder.net/?a=blobdiff_plain;f=R%2FinitSmallEM.R;h=541d7e1cb922b590d1893eea5bba3472eeb1e4ed;hb=f227455a1604906b255ef366d64c10a93e796983;hp=6dd74576eb9d06f610dc12073fbcff3209c6ceda;hpb=c3bc47052f3ccb659659c59a82e9a99ea842398d;p=valse.git diff --git a/R/initSmallEM.R b/R/initSmallEM.R index 6dd7457..541d7e1 100644 --- a/R/initSmallEM.R +++ b/R/initSmallEM.R @@ -3,17 +3,18 @@ #' @param k number of components #' @param X matrix of covariates (of size n*p) #' @param Y matrix of responses (of size n*m) -#' @param tau threshold to stop EM algorithm #' #' @return a list with phiInit, rhoInit, piInit, gamInit #' @export -initSmallEM = function(k,X,Y,tau) +#' @importFrom methods new +#' @importFrom stats cutree dist hclust runif +initSmallEM = function(k,X,Y) { n = nrow(Y) m = ncol(Y) p = ncol(X) - Zinit1 = array(0, dim=c(n,20)) + 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)) @@ -34,8 +35,13 @@ initSmallEM = function(k,X,Y,tau) { Z = Zinit1[,repet] Z_indice = seq_len(n)[Z == r] #renvoit les indices où Z==r - - betaInit1[,,r,repet] = ginv(crossprod(X[Z_indice,])) %*% crossprod(X[Z_indice,], Y[Z_indice,]) + if (length(Z_indice) == 1) { + betaInit1[,,r,repet] = ginv(crossprod(t(X[Z_indice,]))) %*% + crossprod(t(X[Z_indice,]), Y[Z_indice,]) + } else { + betaInit1[,,r,repet] = ginv(crossprod(X[Z_indice,])) %*% + crossprod(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]) @@ -56,7 +62,9 @@ initSmallEM = function(k,X,Y,tau) miniInit = 10 maxiInit = 11 - new_EMG = .Call("EMGLLF_core",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,1e-4) + new_EMG = EMGLLF(phiInit1[,,,repet],rhoInit1[,,,repet],piInit1[repet,],gamInit1[,,repet],miniInit,maxiInit,1,0,X,Y,1e-4) LLFEessai = new_EMG$LLF LLFinit1[repet] = LLFEessai[length(LLFEessai)] }