function [lamout,Z,B,S,lnlik] = box_coxd(X,Y,lamin,guessin);
%
%  Routine to compute MLE of lambda in a multivariate set-up.  Note: each
%  dependent variable is transformed using its own value of lambda.
%
%  If three three arguments are present, this routine computes the
%  profile likelihood using the values in lamin as the transformation
%  parameters.
%
%  If two input arguments are present, then the routine computes the
%  MLE of lambda. 
%
%  If four input arguments are present, then the routine computes the
%  MLE of lambda, and uses the third parameter (lamin) as the
%  initial guess.  The value of the fourth parameter is of no
%  importance. All that matters is whether or not the fourth
%  parameter is present or absent.
%
%  The output arguments are the MLE of lambda, the transformed data matrix,
%  the MLE of B, the MLE of Sigma, and the log likelihood.
%
%
maxy=max(max(Y));
Y=Y/maxy;
C=pinv(X'*X)*X';
[n,d]=size(Y);
nd=n*d;
Ly=log(Y);
lgbar=mean(Ly)';
IM=eye(n)-X*pinv(X'*X)*X';
if nargin==3
lamout=lamin;
  Z=[];
for j=1:d;
  yj=Y(:,j);
  lyj=Ly(:,j);
  if abs(lamin(j)) <= eps;
    zj=lyj;
    else;
      zj=yj.^lamin(j);
      zj=(zj-1)/lamin(j);
   end;
Z=[Z zj];
end
B=C*Z;
A=Z'*IM*Z;
S=A/n;
lnlik=-nd/2-n*log(det(A/n))/2+(lamin-1)'*n*lgbar;
lnlik=lnlik-n*d*log(maxy);
return
end


lam=ones(d,1);
if nargin == 4
  lam=lamin;
end;
del=10;
lik=-1.e20;
while del >=1.e-5;
Z=[];
V1=[];
V2=[];
for j=1:d;
  yj=Y(:,j);
  lyj=Ly(:,j);
  if abs(lam(j)) <= eps;
    zj=lyj;
    vj=lyj.^2/2;
    v2j=lyj.^3/3;
    else;
      zj=yj.^lam(j);
      zj=(zj-1)/lam(j);
      vj=lyj.*zj+lyj/lam(j)-zj/lam(j);
      v2j=lyj.*vj+(zj-lyj-lam(j)*vj)/lam(j)^2;
   end;
   V1=[V1 vj];
   V2=[V2 v2j];
   Z=[Z zj];
end;
BZ=C*Z;
BV=C*V1;
XpZ=X'*Z;
XpV=X'*V1;
A=Z'*Z-BZ'*XpZ;
AV=V1'*V1-BV'*XpV;
Ai=inv(A);
S=A/n;
% disp([lik lam'])
newlik=-nd/2-n*log(det(S))/2+(lam-1)'*n*lgbar;
M0=V1'*Z-BV'*XpZ;
M1=M0*Ai;
f=-diag(M1)+lgbar;
M2=Ai*(Z'*V2-BZ'*(X'*V2));
M3=M1*M0';
F=-(diag(diag(M2))+AV.*Ai-M1'.*M1-M3.*Ai);
newlam=lam-pinv(F)*f;
% disp([f'*pinv(F)*f newlik])
lam=newlam;
del=abs(lik-newlik);
lik=newlik;
end
lnlik=lik;
B=BZ;
lamout=lam;
lnlik=lnlik-n*d*log(maxy);