With these five macros, the multivariate Gaussian model can be estimated with code as demonstrated below. Note that this is not a Gaussian mixture model. I want to introduce the use of these macros to construct and solve the simpler problem first of just a multivariate Gaussian model. Again, I appeal to the sashelp.iris data.

title "Multivariate (k=4) normal model with one cluster";
proc nlmixed data=sashelp.iris(where=(PetalWidth10));
parms mu1 2 mu2 34 mu3 50 mu4 14
s11 2 s12 s13 s14 .5
s22 20 s23 s24 .5
s33 20 s34 .5
s44 10;

/* List variables assumed to have multivariate Gaussian distribution in array */
array _X_ {4} PetalWidth SepalWidth SepalLength PetalLength;

/* List parameters mu in array and construct row vectors and column vectors (X-mu) and T(X-mu) */
array _mu {4} mu1 mu2 mu3 mu4;
array _XminMu {4,1} _temporary_;
array _XminMuT {1,4} _temporary_;
do i=1 to 4;
_xminMu{i,1} = _x_{i} - _mu{i};
_xminMuT{1,i} = _x_{i} - _mu{i};
end;

/* Construct variance matrix. Note that upper triangular part of variance matrix is filled already. */
array _V {4,4} S11 S12 S13 S14
S21 S22 S23 S24
S31 S32 S33 S34
S41 S42 S43 S44;
/* Complete the lower triangular part of the variance matrix */
S21 = S12;
S31 = S13;
S41 = S14;
S32 = S23;
S42 = S24;
S43 = S34;

/* Define matrices holding the Cholesky decomposition of variance L, transpose of L */
/* inverse of T(L), transpose of inv(T(L)), and inverse of the variance matrix */
array _chol {4,4} _temporary_;
array _cholT {4,4} _temporary_;
array _cholInvT {4,4} _temporary_;
array _cholInv {4,4} _temporary_;
array _Vinv {4,4} _temporary_;

/* Now, compute all of these elements */
%mat_chol(input=_V, cholesky=_chol)
%mat_transpose(input=_chol, transpose=_cholT);
%mat_LowerDiagInv(LowDiagInput=_cholT, LowDiagInv=_cholInvT, dim=4)
%mat_transpose(input=_cholInvT, transpose=_cholInv)
%mat_mult(left=_cholInv, right=_cholInvT, result=_Vinv)
%mat_det_chol(cholesky=_chol, det=det)

/* Compute MahalanobisDSq form T(X-mu)*inv(V)*(X-mu) for the cluster */
array _quadLeft {1,4} _temporary_;
array _MahalanobisDSq {1,1} _temporary_;
%mat_mult(left=_XminMuT, right=_Vinv, result=_quadleft)
%mat_mult(left=_quadleft, right=_XminMu, result=_MahalanobisDSq);

/* Compute log-likelihood for the cluster */
LL = -0.5*(log(Det) + _MahalanobisDSq{1,1} + 4*log(2*constant('pi')));

/* Maximize the log-likelihood */
model LL ~ general(LL);
run;