*! version 1.0.0  10aug1996
program define bipr_l2
	version 5.0
	local b "`1'"   /* beta */
	local f "`2'"   /* value of L */
	local g "`3'"   /* gradient */
	local h "`4'"   /* Hessian */

	mac shift 4
	local options "LASTIT FIRSTIT FAST(string)"
	parse "`*'"

	if "`lastit'"!="" { exit }

	tempvar qi1 qi2 wi1 wi2 gi1 gi2 rhoi phi2 t1 t2

	qui gen byte `qi1' = 2*($S_y1!=0) - 1 if $S_mlwgt==1
	qui gen byte `qi2' = 2*($S_y2!=0) - 1 if $S_mlwgt==1

	tempname b1 b2 rho
	matsplit `b' `b1' `b2' `rho'  /* Also defines S_L1 S_L2 */

	if `rho' >= .99 {
		scalar `rho' = .99
	}
	if `rho' <= -.99 {
		scalar `rho' = -.99
	}

	mat score double `wi1' = `b1' if $S_mlwgt==1
	mat score double `wi2' = `b2' if $S_mlwgt==1

	qui replace `wi1' = `qi1'*`wi1' if $S_mlwgt==1
	qui replace `wi2' = `qi2'*`wi2' if $S_mlwgt==1

	qui gen double `rhoi' = `qi1'*`qi2'*`rho' if $S_mlwgt==1

	/* Calculate log likelihood */

	qui gen double `phi2' = binorm(`wi1',`wi2',`rhoi') if $S_mlwgt==1
	qui gen double `t1' = log(`phi2') if $S_mlwgt==1
	qui summ `t1', meanonly
	scalar `f' = _result(18)
	if `fast' == 0  { 
		exit 
	}
	
	drop `t1'

	local nr = rowsof(matrix(`b'))
	local nc = colsof(matrix(`b'))
	matrix `g' = J(1,`nc',0)
	matrix `h' = J(`nc',`nc',0)

	qui gen double `gi1' = normd(`wi1')*/*
		*/ normprob((`wi2'-`rhoi'*`wi1')/sqrt(1-`rhoi'*`rhoi')) /*
		*/ if $S_mlwgt==1
	qui gen double `gi2' = normd(`wi2')*/*
		*/ normprob((`wi1'-`rhoi'*`wi2')/sqrt(1-`rhoi'*`rhoi')) /*
		*/ if $S_mlwgt==1

	/* Calculate gradient */

	qui gen double `t1' = `qi1'*`gi1'/`phi2' if $S_mlwgt==1
	if "$SC_1" != "" {
		qui replace $SC_1 = `t1' if $S_mlwgt==1
	}
	
	qui gen double `t2' = .
	local i 1
	while `i' <= $S_L1 {
		local xx : word `i' of $S_X1
		qui replace `t2' = `t1'*`xx' if $S_mlwgt == 1
		qui summ `t2', meanonly
		mat `g'[1,`i'] = _result(18)
		local i = `i'+1
	}
	qui summ `t1', meanonly
	mat `g'[1,`i'] = _result(18)
	local j = `i'+1

	qui replace `t1' = `qi2'*`gi2'/`phi2' if $S_mlwgt == 1
	if "$SC_2" != "" {
		qui replace $SC_2 = `t1' if $S_mlwgt
	}
	
	local i 1
	while `i' <= $S_L2 {
		local xx : word `i' of $S_X2
		qui replace `t2' = `t1'*`xx' if $S_mlwgt == 1
		qui summ `t2', meanonly
		mat `g'[1,`j'] = _result(18)
		local j = `j'+1
		local i = `i'+1
	}
	qui summ `t1', meanonly
	mat `g'[1,`j'] = _result(18)
	local j = `j'+1

	qui replace `t1' = exp(-.5*(`wi1'^2+`wi2'^2-2*`rhoi'*`wi1'*`wi2')/ /*
	*/ (1-`rhoi'^2))/(2*_pi*sqrt(1-`rhoi'^2)) if $S_mlwgt == 1

	qui replace `t2' = `qi1'*`qi2'*`t1'/`phi2' if $S_mlwgt == 1
	if "$SC_3" != "" {
		qui replace $SC_3 = `t2' if $S_mlwgt == 1
	}
	
	qui summ `t2', meanonly
	mat `g'[1,`j'] = _result(18)
	
	/* Calculate Hessian */

	tempvar di vi1 vi2

	qui gen double `di' = 1/sqrt(1-`rhoi'^2) if $S_mlwgt == 1
	qui gen double `vi1' =  `di'*(`wi2'-`rhoi'*`wi1') if $S_mlwgt == 1
	qui gen double `vi2' =  `di'*(`wi1'-`rhoi'*`wi2') if $S_mlwgt == 1

	qui replace `gi1' = normd(`wi1')*normprob(`vi1') if $S_mlwgt == 1
	qui replace `gi2' = normd(`wi2')*normprob(`vi2') if $S_mlwgt == 1

	qui replace `t1' = `di'*normd(`wi1')*normd(`vi1') if $S_mlwgt == 1

	local i 1
	while `i' <= $S_L1+1 {
		local xi : word `i' of $S_X1
		if "`xi'" == "" {
			local xi = 1 
		}
		
		local j = `i'
		while `j' <= $S_L1+1 {
			local xj : word `j' of $S_X1
			if "`xj'" == "" {
				local xj = 1 
			}
			qui replace `t2' = `xi'*`xj'*(-`wi1'*`gi1'/`phi2' - /*
				*/ `rhoi'*`t1'/`phi2' - `gi1'^2/`phi2'^2) /*
				*/ if $S_mlwgt == 1
			qui summ `t2', meanonly
			mat `h'[`i',`j'] = _result(18)
			mat `h'[`j',`i'] = _result(18)
			local j = `j'+1
		}

		local k 1
		while `k' <= $S_L2+1 {
			local xk : word `k' of $S_X2
			if "`xk'" == "" {
				local xk = 1 
			}
			qui replace `t2' = `qi1'*`qi2'*`xi'*`xk'* /*
				*/ (`t1'/`phi2' - `gi1'*`gi2'/`phi2'^2) /*
				*/ if $S_mlwgt == 1
			qui summ `t2', meanonly
			mat `h'[`i',`j'] = _result(18)
			mat `h'[`j',`i'] = _result(18)
			local k = `k' + 1
			local j = `j'+1
		}

		qui replace `t2' = `qi2'*`xi'*`t1'/`phi2'* /*
			*/ (`rhoi'*`di'*`vi1' - `wi1'-`gi1'/`phi2') /*
			*/ if $S_mlwgt == 1
		qui summ `t2', meanonly
		mat `h'[`i',`j'] = _result(18)
		mat `h'[`j',`i'] = _result(18)
		local i = `i'+1
	}

	local k = $S_L1+2
	local i 1
	while `i' <= $S_L2+1 {
		local xi : word `i' of $S_X2
		if "`xi'" == "" {
			local xi = 1 
		}
		local m = `k'
		local j = `i'
		while `j' <= $S_L2+1 {
			local xj : word `j' of $S_X2
			if "`xj'" == "" {
				local xj = 1 
			}
			qui replace `t2' = `xi'*`xj'*(-`wi2'*`gi2'/`phi2' - /*
				*/ `rhoi'*`t1'/`phi2' - `gi2'^2/`phi2'^2) /*
				*/ if $S_mlwgt == 1
			qui summ `t2', meanonly
			mat `h'[`k',`m'] = _result(18)
			mat `h'[`m',`k'] = _result(18)
			local m = `m'+1	
			local j = `j'+1
		}

		qui replace `t2' = `qi1'*`xi'*`t1'/`phi2'*/*
			*/ (`rhoi'*`di'*`vi2' - `wi2'-`gi2'/`phi2') /*
			*/ if $S_mlwgt == 1
		qui summ `t2', meanonly
		mat `h'[`k',`m'] = _result(18)
		mat `h'[`m',`k'] = _result(18)
		local i = `i'+1
		local k = `k'+1
	}

	local k = $S_L1+$S_L2+3
	qui replace `t2' = `t1'/`phi2'*(`di'^2*`rhoi'*(1-`di'^2*(`wi1'^2+/*
	*/ `wi2'^2-2*`rhoi'*`wi1'*`wi2'))+`di'*`wi1'*`wi2'-`t1'/`phi2') /*
	*/ if $S_mlwgt == 1
	qui summ `t2', meanonly
	mat `h'[`k',`k'] = _result(18)

	tempname minus
	scalar `minus' = -1
	mat `h' = `h'*`minus'
end

program define matsplit
	local mat "`1'"
	local mat1 "`2'"
	local mat2 "`3'"
	local mat3 "`4'"

	local j : coleq(`mat')
	local k : word count `j'
	local km1 = `k'-1
	parse "`j'", parse(" ")
	local i 2
	while "``i''" == "`1'" {
		local i = `i'+1
	}
	local im1 = `i'-1
	mat `mat1' = `mat'[1,1..`im1']
	mat `mat2' = `mat'[1,`i'..`km1']
	scalar `mat3' = `mat'[1,`k']
	global S_L1 = `im1'-1
	global S_L2 = `km1'-`i'
end