function lnLc = LVsarmadouble(rho,n,n2,y,X,Wcan1,Wus1,Wcan2,Wus2)

lambda = rho(3:4);

A = eye(n2) - rho(1)*Wcan1 - rho(2)*Wus1;
% p = colamd(A);
% [l,u] = lu(A(:,p));
% detrho = sum(log(abs(diag(u))));
detrho = log(det(A));

B = eye(n2) + lambda(1)*Wcan2 + lambda(2)*Wus2;
% q = colamd(B);
% [m,v] = lu(B(:,q));
% detlambda = sum(log(abs(diag(v))));
detlambda = log(det(B));

Bi = inv(B);
Ay = A*y;
BAy = Bi*Ay;
BX = Bi*X;
beta = inv(BX'*BX)*BX'*BAy;
e = BAy - BX*beta;

lnLc = - detrho + detlambda + (n2/2)*log((e'*e)/n2);

