prepare structure for R package
[valse.git] / R / initSmallEM.R
1 library(MASS) #generalized inverse of matrix Monroe-Penrose
2
3 vec_bin = function(X,r){
4 Z = c()
5 indice = c()
6 j=1
7 for(i in 1:length(X)){
8 if(X[i] == r){
9 Z[i] = 1
10 indice[j] = i
11 j=j+1
12 }
13 else{
14 Z[i] = 0
15 }
16 }
17 return(list(Z,indice))
18 }
19
20 initSmallEM = function(k,X,Y,tau){
21 n = nrow(Y)
22 m = ncol(Y)
23 p = ncol(X)
24
25 betaInit1 = array(0, dim=c(p,m,k,20))
26 sigmaInit1 = array(0, dim = c(m,m,k,20))
27 phiInit1 = array(0, dim = c(p,m,k,20))
28 rhoInit1 = array(0, dim = c(m,m,k,20))
29 piInit1 = matrix(0,20,k)
30 gamInit1 = array(0, dim=c(n,k,20))
31 LLFinit1 = list()
32
33
34 for(repet in 1:20){
35 clusters = hclust(dist(y)) #default distance : euclidean
36 clusterCut = cutree(clusters,k)
37 Zinit1[,repet] = clusterCut #retourne les indices (à quel cluster indiv_i appartient) d'un clustering hierarchique (nb de cluster = k)
38
39 for(r in 1:k){
40 Z = Zinit1[,repet]
41 Z_bin = vec_bin(Z,r)
42 Z_vec = Z_bin[[1]] #vecteur 0 et 1 aux endroits où Z==r
43 Z_indice = Z_bin[[2]] #renvoit les indices où Z==r
44
45 betaInit1[,,r,repet] = ginv(t(x[Z_indice,])%*%x[Z_indice,])%*%t(x[Z_indice,])%*%y[Z_indice,]
46 sigmaInit1[,,r,repet] = diag(m)
47 phiInit1[,,r,repet] = betaInit1[,,r,repet]/sigmaInit1[,,r,repet]
48 rhoInit1[,,r,repet] = solve(sigmaInit1[,,r,repet])
49 piInit1[repet,r] = sum(Z_vec)/n
50 }
51
52 for(i in 1:n){
53 for(r in 1:k){
54 dotProduct = (y[i,]%*%rhoInit1[,,r,repet]-x[i,]%*%phiInit1[,,r,repet]) %*% (y[i,]%*%rhoInit1[,,r,repet]-x[i,]%*%phiInit1[,,r,repet])
55 Gam[i,r] = piInit1[repet,r]*det(rhoInit1[,,r,repet])*exp(-0.5*dotProduct)
56 }
57 sumGamI = sum(gam[i,])
58 gamInit1[i,,repet]= Gam[i,] / sumGamI
59 }
60
61 miniInit = 10
62 maxiInit = 11
63
64 new_EMG = EMGLLF(phiInit1[,,,repet],rhoInit1[,,,repet],piInit1[repet,],gamInit1[,,repet],miniInit,maxiInit,1,0,x,y,tau)
65 ##.C("EMGLLF", phiInit = phiInit, rhoInit = rhoInit, ...)
66 LLFEessai = new_EMG[[4]]
67 LLFinit1[[repet]] = LLFEessai[[length(LLFEessai)]]
68 }
69
70 b = which.max(LLFinit1)
71
72 phiInit = phiInit1[,,,b]
73 rhoInit = rhoInit1[,,,b]
74 piInit = piInit1[b,]
75 gamInit = gamInit1[,,b]
76
77 return(list(phiInit, rhoInit, piInit, gamInit))
78 }
79
80