- phi[j,mm,r]=0
- else if(S[j,mm,r] > n*lambda*(pi[r]^gamma))
- phi[j,mm,r] = (n*lambda*(pi[r]^gamma)-S[j,mm,r]) / Gram2[j,j,r]
- else
- phi[j,mm,r] = -(n*lambda*(pi[r]^gamma)+S[j,mm,r]) / Gram2[j,j,r]
- }
- }
- }
-
- ##########
- #Etape E #
- ##########
-
- sumLogLLH2 = 0
- for (i in 1:n)
+ phi[j,mm,r]=0
+ else if(S[j,mm,r] > n*lambda*(pi[r]^gamma))
+ phi[j,mm,r] = (n*lambda*(pi[r]^gamma)-S[j,mm,r]) / Gram2[j,j,r]
+ else
+ phi[j,mm,r] = -(n*lambda*(pi[r]^gamma)+S[j,mm,r]) / Gram2[j,j,r]
+ }
+ }
+ }
+
+ ##########
+ #Etape E #
+ ##########
+
+ # Precompute det(rho[,,r]) for r in 1...k
+ detRho = sapply(1:k, function(r) det(rho[,,r]))
+
+ sumLogLLH = 0
+ for (i in 1:n)