b = rep(0, k)
pen = matrix(0, maxi, k)
X2 = array(0, dim=c(n,p,k))
- Y2 = array(0, dim=c(p,m,k))
+ Y2 = array(0, dim=c(n,m,k))
dist = 0
dist2 = 0
ite = 1
#calcul associé à Y et X
for(r in 1:k){
for(mm in 1:m){
- Y2[,mm,r] = sqrt(gam[,r]) .* Y[,mm]
+ Y2[,mm,r] = sqrt(gam[,r]) * Y[,mm] ##bon calcul ? idem pour X2 ??...
}
for(i in 1:n){
- X2[i,,r] = X[i,] .* sqrt(gam[i,r])
+ X2[i,,r] = X[i,] *sqrt(gam[i,r])
}
for(mm in 1:m){
ps2[,mm,r] = crossprod(X2[,,r],Y2[,mm,r])
}
for(j in 1:p){
for(s in 1:p){
- Gram2[j,s,r] = tcrossprod(X2[,j,r], X2[,s,r])
+ Gram2[j,s,r] = crossprod(X2[,j,r], X2[,s,r])
}
}
}
for(r in 1:k){
b[r] = sum(sum(abs(phi[,,r])))
}
- gam2 = sum(gam[1,]) #BIG DOUTE
- a = sum(gam*t(log(Pi)))
+ gam2 = colSums(gam)
+ a = sum(gam%*%(log(Pi)))
#tant que les props sont negatives
kk = 0
#t[m]la plus grande valeur dans la grille O.1^k tel que ce soit
#décroissante ou constante
- while((-1/n*a+lambda*((Pi.^gamma)*b))<(-1/n*gam2*t(log(Pi2))+lambda.*(Pi2.^gamma)*b) && kk<1000){
+ while((-1/n*a+lambda*((Pi^gamma)%*%t(b)))<(-1/n*gam2%*%t(log(Pi2))+lambda*(Pi2^gamma)%*%t(b)) && kk<1000){
Pi2 = Pi+0.1^kk*(1/n*gam2-Pi)
kk = kk+1
}
for(r in 1:k){
for(mm in 1:m){
for(i in 1:n){
- ps1[i,mm,r] = Y2[i,mm,r] * dot(X2(i,:,r), phi(:,mm,r))
+ ps1[i,mm,r] = Y2[i,mm,r] * (X2[i,,r]%*%(phi[,mm,r]))
nY21[i,mm,r] = (Y2[i,mm,r])^2
}
- ps[mm,r] = sum(ps1(:,mm,r));
- nY2[mm,r] = sum(nY21(:,mm,r));
+ ps[mm,r] = sum(ps1[,mm,r])
+ nY2[mm,r] = sum(nY21[,mm,r])
rho[mm,mm,r] = ((ps[mm,r]+sqrt(ps[mm,r]^2+4*nY2[mm,r]*(gam2[r])))/(2*nY2[mm,r]))
}
}
for(r in 1:k){
- for(j in 1:p){
+ p1 = p-1
+ for(j in 1:p1){
for(mm in 1:m){
- S[j,mm,r] = -rho[mm,mm,r]*ps2[j,mm,r] + dot(phi[1:j-1,mm,r],Gram2[j,1:j-1,r]) + dot(phi[j+1:p,mm,r],Gram2[j,j+1:p,r])
- if(abs(S(j,mm,r)) <= n*lambda*(Pi[r]^gamma)){
+ j1 = j-1
+ j2 = j+1
+ v1 = c(1:j1)
+ v2 = c(j2:p)
+ S[j,mm,r] = -rho[mm,mm,r]*ps2[j,mm,r] + phi[v1,mm,r]%*%(Gram2[j,v1,r]) + phi[v2,mm,r]%*%(Gram2[j,v2,r]) #erreur indice
+ if(abs(S[j,mm,r]) <= n*lambda*(Pi[r]^gamma)){
phi[j,mm,r]=0
}else{
if(S[j,mm,r]> n*lambda*(Pi[r]^gamma)){
sumPen = 0
for(r in 1:k){
- sumPen = sumPen + Pi[r].^gamma^b[r]
+ sumPen = sumPen + Pi[r]^gamma^b[r]
}
LLF[ite] = -(1/n)*sumLogLLF2 + lambda*sumPen
else
dist = (LLF[ite]-LLF[ite-1])/(1+abs(LLF[ite]))
- Dist1=max(max(max((abs(phi-Phi))./(1+abs(phi)))))
- Dist2=max(max(max((abs(rho-Rho))./(1+abs(rho)))))
- Dist3=max(max((abs(Pi-PI))./(1+abs(PI))))
- dist2=max([Dist1,Dist2,Dist3])
+ Dist1=max(max(max((abs(phi-Phi))/(1+abs(phi)))))
+ Dist2=max(max(max((abs(rho-Rho))/(1+abs(rho)))))
+ Dist3=max(max((abs(Pi-PI))/(1+abs(PI))))
+ dist2=max(c(Dist1,Dist2,Dist3))
ite=ite+1
}
- Pi = transpose(Pi)
- return(list(phi=phi, rho=rho, Pi=Pi, LLF=LLF, S=S))
-}
\ No newline at end of file
+ Pi = t(Pi)
+ return(list("phi"=phi, "rho"=rho, "pi"=Pi, "LLF"=LLF, "S"=S))
+}