function [l, lvec] = maxbivnorm(theta,data,weight)

l=1;
mu_x  = theta(l); l=l+1;
mu_y  = theta(l); l=l+1;
sig_x = exp(theta(l)); l=l+1;
sig_y = exp(theta(l)); l=l+1;
rho   = (exp(theta(l))-1)/(exp(theta(l))+1);

gamma_x = 1-rho*sig_x/sig_y;
gamma_y = 1-rho*sig_y/sig_x;

if gamma_x == 0 
    mu_x_tild = mu_x-mu_y;
    sig_x_tild = sig_x*(1-rho^2)^0.5;
else
    mu_x_tild = (mu_x-rho*sig_x/sig_y*mu_y)/gamma_x;
    sig_x_tild = sig_x*(1-rho^2)^0.5/abs(gamma_x);
end

if gamma_y == 0 
    mu_y_tild = mu_y-mu_x;
    sig_y_tild = sig_y*(1-rho^2)^0.5;
else
    mu_y_tild = (mu_y-rho*sig_y/sig_x*mu_x)/gamma_y;
    sig_y_tild = sig_y*(1-rho^2)^0.5/abs(gamma_y);
end 

l = (1/sig_x)*normpdf((mu_x-data)/sig_x).*(1-normcdf((mu_y_tild-data)/sig_y_tild)) + ...
   (1/sig_y)*normpdf((mu_y-data)/sig_y).*(1-normcdf((mu_x_tild-data)/sig_x_tild));
lvec = l;
l = -sum(log(l).*weight);
end