*! 1.1.0 GML 08 January 2004 * 1.0.0 GML 10 January 2002 pro def bnormpdf version 8 syntax varlist(min=2 max=2) [if] [in], Dens(string) [ M1(real 0) M2(real 0) /* */ S1(real 1) S2(real 1) Rho(real 0) REPLACE DOUble ] tokenize `varlist' local x1 "`1'" local x2 "`2'" marksample touse cap assert `s1' > 0 & `s2' > 0 if _rc~=0 { di as error "args to s1() and s2() must be > 0" error _rc } cap assert abs(`rho') < 1 if _rc~=0 { di as error "arg to rho() must be between -1 and 1" error _rc } if "`dens'"!=`""' { local varnum : word count `dens' cap assert `varnum' == 1 if _rc~=0 { di as error "more than one newvar name in dens( ) option" error _rc } } if `"`replace'"'!=`""' { cap drop `dens' } else { cap confirm new variable `dens' if _rc~=0 { di as error "variable `dens' already exists," di as error `" use "replace" option to replace `dens'"' error _rc } } /* it seems Stata matrix size limitations won't allow calculating this using matrix notation , i.e. can't calc the quadratic form: (X - mu)Sigma_inverse(X - mu)' ? The following form is from Bickel and Doksum, p.23-4 */ #delimit ; gen `double' `dens' = (1.0/(2.0*_pi*`s1'*`s2'*sqrt(1.0 - `rho'^2))) * exp( (-1/(2*(1.0-`rho'^2))) * ( (((`x1'-`m1')/`s1')^2) - (2 *`rho'* ((`x1'-`m1')/`s1') * ((`x2'-`m2')/`s2')) + (((`x2'-`m2')/`s2')^2) ) ) if `touse' ; #delimit cr end