function [lam_out,Z,B,S,lnlik] = box_coxdh(X,Y,lam_in,grp,lam_guess)
%
%  Routine to compute MLE of lambda in a multivariate set-up.  Note: each
%  dependent variable is transformed using a distinct value for lambda.
%
%  If a value for lam_in (i.e., lambda input) is given, then this
%  routine computes the profile likelihood using the value in
%  lam_in as the transformation parameter.
%
%  If lam_guess is empty, then iteration begins with lam = 1.
%  Otherwise, iteration begins with lam = lam_guess.
%
%  The input variable grp identifies a group structure 
%  to allow for heterogeneous covariance matrices.  The variable
%  grp is an n x 1 vector whose components designate the group
%  membership for each case.  The X matrix is partitioned according
%  to the group structure and each group is allowed its own
%  covariance matrix.  If there is no group structure, then set
%  grp = [] or set grp to an n-vector of ones (i.e., all cases are 
%  in group 1).
%
%  The function name is meant to be a reminder that Box-Cox parameters are
%  estimated in a d-dimensional multivariate setting, where
%  covariance matrices may be heterogeneous, and a separate lambda
%  is estimated for each variable.
%
[n,d]=size(Y);
if isempty(lam_guess) == 1
  lam_guess=ones(d,1);
end


const=max(max(Y));
Y=Y/const;
Z0=log(Y);
nd=n*d;
if isempty(grp) == 1
grp=ones(n,1);
end;

g=max(grp);
for j=1:g
  ii=find(grp==j);
  YY{j}=Y(ii,:);
  XX{j}=X(ii,:);
  nn(j)=length(ii);
  IM{j}=eye(nn(j))-ppo(XX{j});
  C{j}=pinv(XX{j}'*XX{j})*XX{j}';
  ZZ0{j}=Z0(ii,:);
  lnG{j}=mean(log(YY{j}))';
end;
if isempty(lam_in)==0
lam_out=lam_in;
lnlik=0;
for j=1:g
  for l1=1:d 
  if abs(lam_in(l1)) <= 1.e-10
    Z{j}(:,l1)=log(YY{j}(:,l1));
  else
    Z{j}(:,l1)=(YY{j}(:,l1).^lam_in(l1)-1)/lam_in(l1);
  end
end
B{j}=C{j}*Z{j};
A{j}=Z{j}'*IM{j}*Z{j};
S{j}=A{j}/nn(j);
lnlik=lnlik-nn(j)*d/2-nn(j)*log(det(S{j}))/2+nn(j)*(lam_in-1)'*lnG{j};
end
lnlik=lnlik-n*d*log(const);
return
end

lam=lam_guess;
del=1;
lnlik=-1.e10;

while del >1.e-5
  newlnlik=0;
  h=zeros(d,1);
  H=zeros(d,d);
for j=1:g
  for l1=1:d
  if abs(lam(l1)) <=1.e-10
    Z{j}(:,l1)=ZZ0{j}(:,l1);
    Z1{j}(:,l1)=Z{j}(:,l1).^2/2;
    Z2{j}(:,l1)=Z{j}(:,l1).^3/3;
   else
     Y0{j}(:,l1)=YY{j}(:,l1).^lam(l1);
     Z{j}(:,l1)=(Y0{j}(:,l1)-1)/lam(l1);
     Y1{j}(:,l1)=Y0{j}(:,l1).*ZZ0{j}(:,l1);
     Z1{j}(:,l1)=Y1{j}(:,l1)/lam(l1)-Z{j}(:,l1)/lam(l1);
     Z2{j}(:,l1)=Y1{j}(:,l1).*ZZ0{j}(:,l1)/lam(l1)-2*Z1{j}(:,l1)/lam(l1);
   end
   end
A00{j}=Z{j}'*IM{j}*Z{j};
A10{j}=Z1{j}'*IM{j}*Z{j};
A11{j}=Z1{j}'*IM{j}*Z1{j};
A20{j}=Z2{j}'*IM{j}*Z{j};
S{j}=A00{j}/nn(j);
Si{j}=inv(S{j});
newlnlik=newlnlik -nn(j)*d/2-nn(j)*log(det(S{j}))/2+nn(j)*(lam-1)'*lnG{j};
L1{j}=Si{j}*A10{j}';
L11{j}=A10{j}*Si{j}*A10{j}';
L2{j}=Si{j}*A20{j}';

for l1=1:d
  h(l1)=h(l1)-L1{j}(l1,l1)+nn(j)*lnG{j}(l1);
  for l2=1:l1
    H(l1,l2)=H(l1,l2)+L1{j}(l1,l2)*L1{j}(l2,l1)/nn(j) ...
	+L11{j}(l1,l2)*Si{j}(l1,l2)/nn(j) ...
       -Si{j}(l1,l2)*A11{j}(l1,l2) -(l1==l2)*L2{j}(l1,l2);
   H(l2,l1)=H(l1,l2);
  end
end
B{j}=C{j}*Z{j};
end
% disp(['h'])
% disp(h')
% disp(['H'])
% disp(H)
% eig(H)'
% disp('newlnlik')
% disp(newlnlik)

% disp([lam del])
% lam-h/H
lam=lam-H\h;
% disp('newlam')
% disp(lam')
del = abs(lnlik-newlnlik);
% pause
% disp([lam newlnlik])
% pause
lnlik=newlnlik;
end
lam_out=lam;
lnlik=lnlik-n*d*log(const);