Hey Rick,
thank you for the fast reply. I'am ashamed about point two and three, sorry for that. With respect to the mahalanobis distance: I had notice that there is an inbuilt function for calculating this distance in IML. However, because I will avoid that a singular covariance matrix is optained during the process of optimization (which will produce error message in the log), I have adopt the mahalanobis distance to reject such solutions (see adoptions of me in the modul).
Underneath is the code with your suggestions included. However, runig this code will produce an error message in the log "NEWRAP Optimization cannot be completed." and a warning message "Optimization routine cannot improve the function value.". The solution mue = {-5.7 4.97} and cove = {0.05 -0.83, -0.83 12.55} is far away from the sample estimates mu = {3 4} and cov = {8 1.1, 1.1 11.1}.
I have comparde my objective function of the MVN with your pdf function of the MVN (http://blogs.sas.com/content/iml/2012/07/05/compute-multivariate-normal-denstity.html) and I get the same result (when I log your solution). So I assume it is not a problem of the objective function. I also compare my gradient function (GradLogLik) with the results of finit difference approximation and the results are the same. Hence, I assume it is also not a problem of the gradient function. Do you have any further suggestions?
proc iml;
/*** Write a module that computes the log-likelihood ***/
/** Mahalanobis module from Rick Wicklin (2012) **/
/** Cut and paste the definition of the Mahalanobis module. **/
/** See http://blogs.sas.com/content/iml/2012/02/22/how-to-compute-mahalanobis-distance-in-sas/ **/
start mahalanobis(x, center, cov);
n = nrow(x);
m = nrow(center);
d2 = j(n,m);
U = root(cov, "NoError");
if all(u=.) then return(1E10); /* adoption of me: if cov is singular, this will force NLP to try another estimate of cov (see Gu & Yung, 2013: A SAS/IML program using the Kalman filter for estimating state space models) */
Ur = round(trace(ginv(U)*U));
if Ur < nrow(U) then return(1E10); /* adoption of me: if U is singular, this will force NLP to try another estimate of cov */
do k = 1 to m;
v = x - center[k,];
z = trisolv(2, U, v`);
w = trisolv(1, U, z);
d2[,k] = (v # w`)[,+];
end;
return (sqrt(d2));
finish;
start LogLikMVN(param) global(x,nmu);
mu = t(param[1:nmu]);
cov1 = t(param[(nmu+1):ncol(param)]);
cov = sqrvech(cov1);
rowx = nrow(x);
eins = j(1,rowx,1);
d = mahalanobis(x, mu, cov);
if d=1E10 then return(1E10); /* if cov is singular, this will force NLP to try another estimate of cov */
dd = d#d;
dsd = eins * dd;
logdet1 = logabsdet(cov);
logdet = choose(logdet1[2],logdet1[1] * logdet1[2], 0);
f = - (rowx*nmu/2 * log(2*constant("PI"))) - (rowx/2 * logdet) - (1/2 * dsd);
return( f );
finish;
/*** Write a module that computes the gradient of the log-likelihood function ***/
start GradLogLik(param) global (x,nmu);
mu = t(param[1:nmu]);
cov1 = t(param[(nmu+1):ncol(param)]);
cov = sqrvech(cov1);
rowx = nrow(x);
colx = ncol(x);
eins = j(1,rowx,1);
xdif = x - mu;
xdifs = t(xdif[+,]);
invcov=inv(cov);
dfdmu = t(invcov*xdifs);
B = j(colx,colx,0);
do i=1 to rowx;
xred = xdif[i,];
A = t(xred)*xred;
B = B + A;
end;
dfdcov1 = -0.5 * (invcov - invcov*B*invcov);
dfdcov = t(vech(dfdcov1));
return ( dfdmu || dfdcov );
finish;
/*** Call an optimization routine ***/
nmu = {2};
x = {0 4,
2 2,
1 4,
6 3,
2 1,
1 0,
3 8,
5 1,
9 7,
1 10};
con = { . . 0 . 0,
. . . . .};
p = {2 3 7 1 10};
opt = {0,5};
call nlpnra(rc, result, "LogLikMVN", p, opt, con) grd="GradLogLik";
quit;
