Expectation maximization algorithm for normalized Gaussian mixtures

Input: data set 𝒜 ⊂ (ℝ^n), initial centers (Subscript[z, 1]),⃛,(Subscript[z, k]) ∈ (ℝ^n);
Output: partition Π={(Subscript[π, 1]),⃛,(Subscript[π, k])}, centers (Subscript[c, 1]),⃛,(Subscript[c, k]) ∈ (ℝ^n), covariance matrices (∑_j)∈(ℝ^(n×n)), j=1⃛k, number of iterations It;
Submodules: dm[], MWF[], MFmin[], Hdist[]

(* It is used dm distance-like function *)
EM[w_,A_,z_,Ind_:0]:=Module[{m=Length[A],n=Length[A[[1]]],k=Length[z],It=0,MaxIt=30,eps=.005,z0},
				(* Initialization *)
z0=z;Sigma0=Table[IdentityMatrix[n],{i,k}];P0=Table[1/k,{i,k}];lam=Table[0,{i,m}];
p[u_,mu_,Cov_]:=Exp[-(u-mu).Inverse[Cov].(u-mu)/2]/(Sqrt[2Pi Det[Cov]]);
				(* Iterative process *)
While[It=It+1;
PI=Table[{},{j,k}];
Do[
lam[[i]]=Table[(P0[[j]] p[A[[i]],z0[[j]],Sigma0[[j]]])/Sum[p[A[[i]],z0[[s]],Sigma0[[s]]]P0[[s]],{s,k}],{j,k}]
,{i,m}];
Do[
imin=Ordering[lam[[i]],-1][[1]];
PI[[imin]]=Append[PI[[imin]],A[[i]]],
{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI];
If[Ind≠0,
slPI=ListPlot[Table[PI[[j]],{j,k}],PlotStyle->Table[Hue[j/k],{j,k}]];
sl0=ListPlot[{z0},PlotStyle->{Brown,Opacity[.8],PointSize[.025]}];
Do[
rad=Max[Table[dm[z0[[j]],PI[[j,s]],Sigma0[[j]]],{s,Length[PI[[j]]]}]];
slcon[j]=ContourPlot[dm[z0[[j]],{x,y},Sigma0[[j]]]==75rad/100,{x,xG[[1]],xG[[2]]},{y,yG[[1]],yG[[2]]}],{j,k}];
sl=Show[Table[slcon[j],{j,k}]];
Print[Show[slPI,sl,sl0,Ticks->{{-.5,.5},{-.5,.5}},AxesOrigin->{0,0},ImageSize->Small]];
];
z1=Table[Sum[lam[[i,j]]A[[i]],{i,m}]/Sum[lam[[i,j]],{i,m}],{j,k}];
Chop[Hdist[z0,z1]]>eps&&It≤ MaxIt,
Sigma1=Table[Sum[lam[[i,j]]KroneckerProduct[(A[[i]]-z0[[j]]),((A[[i]]-z0[[j]]))],{i,m}]/Sum[lam[[i,j]],{i,m}],{j,k}];
P1=Table[Sum[lam[[i,j]],{i,m}]/m,{j,k}];
If[Ind≠0,
Print[It,". ITERACIJA: \nCenters: ",z1,"\nCov: ",Table[Sigma1[[j]]//MatrixForm,{j,Length[z0]}],"\nP: ",P1,"\nHdist= ",Hdist[z0,z1]];
];
P0=P1;Sigma0=Sigma1;z0=z1;
];
{PI,z1,Sigma1,It,P1}
]