Bo-Cox for multivariate normalization

[robmio]

The proposed problem, for a single column matrix, has been solved with the Newton Raphson method. In the program below, "gij" is the gradient, "Hij" the Hessian matrix.
I find it difficult to extend the program to matrices with more than one column (multivariate matrices). Is there perhaps a command in HP Prime to optimize the function described in the first post?

Code:


#cas
BoxCox_Newton_Cas(r,λo,iter):=
BEGIN
LOCAL y, z, u, v, x, Ii, Hh, sσ2;
LOCAL dmsr, j, gij, Hij;
y:=mat2list(r);
dmsr:=rowDim(r);
x:=MAKEMAT(1,dmsr,1);
Hh:=x*inv(TRN(x)*x)*TRN(x);
Ii:=IDENMAT(dmsr);
FOR j FROM 1 TO iter DO
IF λo≠0 THEN
z:=(y.^λo-1)./λo;
ELSE
z:=LN(y);
END;
IF λo≠0 THEN
u:=(λo*y.^λo*LN(y)-y.^λo+1)./λo.^2;
ELSE
u:=LN(y).^2/2;
END;
IF λo≠0 THEN
v:=(2-((λo*LN(y)).^2-2*λo*LN(y)+2)*
y.^λo)./λo.^3;
ELSE
v:=−LN(y).^3/3;
END;
z:=list2mat(z,1);
sσ2:=(TRN(z)*(Ii-Hh)*z)/dmsr;
u:=list2mat(u,1);
gij:=−(TRN(u)*(Ii-Hh)*z)*inv(sσ2)+
ΣLIST(LN(y));
v:=list2mat(v,1);
Hij:=(TRN(v)*(Ii-Hh)*z-TRN(u)*(Ii-Hh)*
u)*inv(sσ2)+2/dmsr*((TRN(u)*(Ii-Hh)*
z)*inv(sσ2)).^2; 
λo:=λo-inv(Hij)*gij;
λo:=GET(GET(λo,1),1);
END;
RETURN {λo,gij};
END;
#end