Bivariate Normal Distribution

The bivariate normal distribution is the statistical distribution with probability density function

(1)

where

(2)

and

(3)

is the correlation of and (Kenney and Keeping 1951, pp. 92 and 202-205; Whittaker and Robinson 1967, p. 329) and is the covariance.

The probability density function of the bivariate normal distribution is implemented as MultinormalDistribution[mu1, mu2, sigma11, sigma12, sigma12, sigma22] in the Wolfram Language package MultivariateStatistics` .

The marginal probabilities are then

			(4)
			(5)

and

			(6)
			(7)

(Kenney and Keeping 1951, p. 202).

Let and be two independent normal variates with means and for , 2. Then the variables and defined below are normal bivariates with unit variance and correlation coefficient :

			(8)
			(9)

To derive the bivariate normal probability function, let and be normally and independently distributed variates with mean 0 and variance 1, then define

			(10)
			(11)

(Kenney and Keeping 1951, p. 92). The variates and are then themselves normally distributed with means and , variances

			(12)
			(13)

and covariance

(14)

The covariance matrix is defined by

(15)

where

(16)

Now, the joint probability density function for and is

(17)

but from (◇) and (◇), we have

(18)

As long as

(19)

this can be inverted to give

			(20)
			(21)

Therefore,

x_1^2+x_2^2=([sigma_(22)(y_1-mu_1)-sigma_(12)(y_2-mu_2)]^2)/((sigma_(11)sigma_(22)-sigma_(12)sigma_(21))^2)
+([-sigma_(21)(y_1-mu_1)+sigma_(11)(y_2-mu_2)]^2)/((sigma_(11)sigma_(22)-sigma_(12)sigma_(21))^2),

(22)

and expanding the numerator of (22) gives

(23)

(x_1^2+x_2^2)(sigma_(11)sigma_(22)-sigma_(12)sigma_(21))^2
=(y_1-mu_1)^2(sigma_(21)^2+sigma_(22)^2)-2(y_1-mu_1)(y_2-mu_2)(sigma_(11)sigma_(21)+sigma_(12)sigma_(22))+(y_2-mu_2)^2(sigma_(11)^2+sigma_(12)^2)
=sigma_2^2(y_1-mu_1)^2-2(y_1-mu_1)(y_2-mu_2)(rhosigma_1sigma_2)+sigma_1^2(y_2-mu_2)^2
=sigma_1^2sigma_2^2[((y_1-mu_1)^2)/(sigma_1^2)-(2rho(y_1-mu_1)(y_2-mu_2))/(sigma_1sigma_2)+((y_2-mu_2)^2)/(sigma_2^2)].

(24)

Now, the denominator of (◇) is

(25)

			(26)
			(27)
			(28)

can be written simply as

(29)

and

(30)

Solving for and and defining

(31)

gives

			(32)
			(33)

But the Jacobian is

			(34)
			(35)
			(36)

(37)

and

(38)

where

(39)

Q.E.D.

The characteristic function of the bivariate normal distribution is given by

			(40)
			(41)

where

(42)

and

(43)

Now let

			(44)
			(45)

Then

(46)

where

			(47)
			(48)

Complete the square in the inner integral

$int_(-infty)^inftyexp{-1/(2(1-rho^2))1/(sigma_1^2)[u^2-(2rhosigma_1w)/(sigma_2)u]}e^(t_1u)du =int_(-infty)^inftyexp{-1/(2sigma_1^2(1-rho^2))[u-(rho_1sigma_1w)/(sigma^2)]^2}{1/(2sigma_1^2(1-rho^2))((rho_1sigma_1w)/(sigma_2))^2}e^(it_1u)du.$

(49)

Rearranging to bring the exponential depending on outside the inner integral, letting

(50)

and writing

(51)

gives

$phi(t_1,t_2)=N^'int_(-infty)^inftye^(it_2w)exp[-1/(2sigma_2^2(1-rho^2))w^2]exp[(rho^2)/(2sigma_2^2(1-rho^2))w^2]int_(-infty)^inftyexp[-1/(2sigma_2^2(1-rho^2))v^2]{cos[t_1(v+(rhosigma_1w)/(sigma_2))]+isin[t_1(v+(rhosigma_1w)/(sigma_2))]}dvdw.$

(52)

Expanding the term in braces gives

[cos(t_1v)cos((rhosigma_1wt_1)/(sigma_2))-sin(t_1v)sin((rhosigma_1w)/(sigma_2t_1))]+i[sin(t_1v)cos((rhosigma_1w)/(sigma_2t_1))+cos(t_1v)sin((rhosigma_1wt_1)/(sigma_2))]
=[cos((rhosigma_1wt_1)/(sigma_2))+isin((rhosigma_1wt_1)/(sigma_2))][cos(t_1v)+isin(t_1v)]=exp((irhosigma_1w)/(sigma_2)t_1)[cos(t_1v)+isin(t_1v)].

(53)

But is odd, so the integral over the sine term vanishes, and we are left with

phi(t_1,t_2)=N^'int_(-infty)^inftye^(it_2w)exp[-(w^2)/(2sigma_2^2)]exp[(rho^2w^2)/(2sigma_2^2(1-rho^2))]exp[(irhosigma_1wt_1)/(sigma_2)]dwint_(-infty)^inftyexp[-(v^2)/(2sigma_1^2(1-rho^2))]cos(t_1v)dv
=N^'int_(-infty)^inftyexp[iw(t_2+t_1(rho(sigma_1)/(sigma_2)))]exp[-(w^2)/(2sigma_2^2)]dwint_(-infty)^inftyexp[-(v^2)/(2sigma_1^2(1-rho^2))]cos(t_1v)dv.

(54)

Now evaluate the Gaussian integral

			(55)
			(56)

to obtain the explicit form of the characteristic function,

$phi(t_1,t_2)=(e^(i(t_1mu_1+t_2mu_2)))/(2pisigma_1sigma_2sqrt(1-rho^2)){sigma_2sqrt(2pi)exp[-1/4(t_2+rho(sigma_1)/(sigma_2)t_1)^22sigma_2^2]}{sigma_1sqrt(2pi(1-rho^2))exp[-1/4t_1^22sigma_1^2(1-rho^2)]} =e^(i(t_1mu_1+t_2mu_2))exp{-1/2[t_2^2sigma_2^2+2rhosigma_1sigma_2t_1t_2+rho^2sigma_1^2t_1^2+(1-rho^2)sigma_1^2t_1^2]} =exp[i(t_1mu_1+t_2mu_2)-1/2(sigma_1^2t_1^2+2rhosigma_1sigma_2t_1t_2+sigma_2^2t_2^2)].$