The Standard BVN Problem

The standard bivariate normal distribution function is given by

$$\Phi({\bf b},\rho)=\frac{1}{2\pi\sqrt{1-\rho^{2}}}\int_{-\in... ...nfty}^{b_2} e^{-(x^{2}-2\rho xy +y^{2})/(2(1-\rho^{2}))}dydx,$$

where ${\bf b}= (b_1,b_2)$. Early work on BVN computations (see Andel, 1974, and Terza and Welland, 1990) studied the bivariate normal probability defined by

$$L(h,k,\rho) = \frac{1}{2\pi\sqrt{1-\rho^{2}}} \int_h^\infty\int_k^{\infty}e^{-(x^2-2\rho xy +y^2)/(2(1-\rho^{2}))}dydx,$$ (1)

which is related to the standard bivariate normal distribution function by $\Phi({\bf b},\rho)=L(-b_1,-b_2,\rho)$. In this section the discussion will be focused on methods for the computation of $L(h,k,\rho)$, in order to allow consistent references to earlier work. Drezner and Wesolowsky (1990) studied the formula

$$L(h,k,\rho) = \Phi(-h)\Phi(-k) + \frac{1}{2\pi} \int_0^\rho \frac{1}{\sqrt{1-r^2}}e^{-\frac{h^2+k^2-2rhk}{2(1-r^2)}}dr,$$

and used numerical integration for computation of BVN probabilities. The formula derived by Plackett (1954) for the correlation coefficient partial derivative of the bivariate normal distribution can be written

$$\frac{\partial L(h,k,r)}{\partial r} = \frac{e^{-\frac{h^2-2r hk+k^{2}}{2(1-r^2)}}}{2\pi\sqrt{1-r^2}}.$$ (2)

The integration of this formula for $r$ between $0$ and $\rho$ produces the formula studied by Drezner and Wesolowsky.