51485a7d |
1 | constructionModelesLassoMLE = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma, |
2 | X,Y,seuil,tau,selected, parallel = FALSE) |
46a2e676 |
3 | { |
51485a7d |
4 | if (parallel) { |
5 | #TODO: parameter ncores (chaque tâche peut aussi demander du parallélisme...) |
6 | cl = parallel::makeCluster( parallel::detectCores() / 4 ) |
7 | parallel::clusterExport(cl=cl, |
8 | varlist=c("phiInit","rhoInit","gamInit","mini","maxi","X","Y","seuil","tau"), |
9 | envir=environment()) |
10 | #Pour chaque lambda de la grille, on calcule les coefficients |
11 | out = parLapply( seq_along(glambda), function(lambda) |
12 | { |
13 | n = dim(X)[1] |
14 | p = dim(phiInit)[1] |
15 | m = dim(phiInit)[2] |
16 | k = dim(phiInit)[3] |
17 | |
18 | #TODO: phiInit[selected] et X[selected] sont bien sûr faux; par quoi remplacer ? |
19 | #lambda == 0 c'est normal ? -> ED : oui, ici on calcule le maximum de vraisembance, donc on ne pénalise plus |
20 | res = EMGLLF(phiInit[selected],rhoInit,piInit,gamInit,mini,maxi,gamma,0.,X[selected],Y,tau) |
21 | |
22 | #comment évaluer la dimension à partir du résultat et de [not]selected ? |
23 | #dimension = ... |
24 | |
25 | #on veut calculer la vraisemblance avec toutes nos estimations |
26 | densite = vector("double",n) |
27 | for (r in 1:k) |
28 | { |
29 | delta = Y%*%rho[,,r] - (X[selected]%*%res$phi[selected,,r]) |
30 | densite = densite + pi[r] * |
31 | det(rho[,,r])/(sqrt(2*base::pi))^m * exp(-tcrossprod(delta)/2.0) |
32 | } |
33 | llh = c( sum(log(densite[,lambda])), (dimension+m+1)*k-1 ) |
34 | list("phi"=res$phi, "rho"=res$rho, "pi"=res$pi, "llh" = llh) |
35 | }) |
36 | parallel::stopCluster(cl) |
37 | out |
38 | } |
39 | else { |
40 | #Pour chaque lambda de la grille, on calcule les coefficients |
41 | n = dim(X)[1] |
42 | p = dim(phiInit)[1] |
43 | m = dim(phiInit)[2] |
44 | k = dim(phiInit)[3] |
45 | L = length(selected) |
46 | phi = list() |
47 | phiLambda = array(0, dim = c(p,m,k)) |
48 | rho = list() |
49 | pi = list() |
50 | llh = list() |
51 | |
3f62d540 |
52 | out = lapply( seq_along(selected), function(lambda) |
53 | { |
64b28e3e |
54 | print(lambda) |
51485a7d |
55 | sel.lambda = selected[[lambda]] |
56 | col.sel = which(colSums(sel.lambda)!=0) |
64b28e3e |
57 | if (length(col.sel)>0){ |
58 | res_EM = EMGLLF(phiInit[col.sel,,],rhoInit,piInit,gamInit,mini,maxi,gamma,0.,X[,col.sel],Y,tau) |
59 | phiLambda2 = res_EM$phi |
60 | rhoLambda = res_EM$rho |
61 | piLambda = res_EM$pi |
62 | for (j in 1:length(col.sel)){ |
63 | phiLambda[col.sel[j],,] = phiLambda2[j,,] |
51485a7d |
64 | } |
64b28e3e |
65 | |
66 | dimension = 0 |
67 | for (j in 1:p){ |
68 | b = setdiff(1:m, sel.lambda[,j]) |
69 | if (length(b) > 0){ |
70 | phiLambda[j,b,] = 0.0 |
71 | } |
72 | dimension = dimension + sum(sel.lambda[,j]!=0) |
73 | } |
74 | |
75 | #on veut calculer la vraisemblance avec toutes nos estimations |
76 | densite = vector("double",n) |
77 | for (r in 1:k) |
78 | { |
79 | delta = Y%*%rhoLambda[,,r] - (X[, col.sel]%*%phiLambda[col.sel,,r]) |
80 | densite = densite + piLambda[r] * |
81 | det(rhoLambda[,,r])/(sqrt(2*base::pi))^m * exp(-tcrossprod(delta)/2.0) |
82 | } |
83 | llhLambda = c( sum(log(densite)), (dimension+m+1)*k-1 ) |
84 | list("phi"= phiLambda, "rho"= rhoLambda, "pi"= piLambda, "llh" = llhLambda) |
51485a7d |
85 | } |
51485a7d |
86 | } |
3f62d540 |
87 | ) |
88 | return(out) |
51485a7d |
89 | } |
c3bc4705 |
90 | } |