Sample Probability Calculations

Sample Probability Calculations

We can calculate probabilities for the Bivariate Normal distribution using either the joint density

or the joint distribution function. We will denote the density by f and the distribution function by F .

The next bit of code will draw the joint distribution function for = 0 and = 1. Warning : Maple doesn't seem to like evaluating the Doubleint at , so don't change the value of .

> restart:

> with(plots,display,textplot3d): with(student):

> mu[1]:=0; mu[2]:=0; sigma[1]:=1; sigma[2]:=1; rho:=0:

> f(s,t):=exp((-1/(2*(1-rho^2)))*(((s-mu[1])/sigma[1])^2-2*rho*(s-mu[1])*(t-mu[2])/(sigma[1]*sigma[2])+((t-mu[2])/sigma[2])^2))/(2*Pi*sigma[1]*sigma[2]*sqrt(1-rho^2)):

The joint distribution function appears below.

> F(x,y):=Doubleint(f(s,t),s=-infinity..x,t=-infinity..y);

> plot3d(value(F(x,y)),x=-3..3, y=-3..3, axes=frame, title="Bivariate Normal Joint Distribution, rho=0");

Now, try and find the value of the yellow shaded volume, i.e. everything under the surface where x < 0.5 and y < 1.5.

> f(x,y):=exp((-1/(2*(1-rho^2)))*(((x-mu[1])/sigma[1])^2-2*rho*(x-mu[1])*(y-mu[2])/(sigma[1]*sigma[2])+((y-mu[2])/sigma[2])^2))/(2*Pi*sigma[1]*sigma[2]*sqrt(1-rho^2)):

> pic[1]:=plot3d(f(x,y),x=-3.5..3.5,y=-3.5..3.5,axes=framed,color=black,style=line,title="Volume under the joint density"):

> pic[2]:=plot3d(f(x,y),x=-3.5..0.5,y=-3.5..1.5,color=yellow, filled=true,style=PATCHNOGRID):

> display([pic[1],pic[2]]);

> value(subs(x=0.5,y=1.5,F(x,y)));

If you wish, you can also find the value of the volume by integrating the density function.

> value(Doubleint(f(x,y),x=-infinity..0.5,y=-infinity..1.5));

Now try and find the value of the shaded volume, i.e. the volume over < 0.9 and 0 < y < 1.

> pic[1]:=plot3d(f(x,y),x=-3.5..3.5,y=-3.5..3.5,axes=framed,color=black,style=line,title="volume under the joint density"):

> pic[2]:=plot3d(f(x,y),x=-1..0.9,y=0..1,color=yellow, filled=true,style=PATCHNOGRID):

> display([pic[1],pic[2]]);

Using the joint distribution function, we would find the volume over < 0.9 and < 1, and subtract the volume where x < -1 and y > 1 or y < 0.

> value(subs(x=0.9,y=1,F(x,y)));

Make sure you understand why the code is written the way it is below. Hint: You don't want to double count a certain volume.

> evalf(%-(value(subs(x=-1,y=1,F(x,y))))-(value(subs(x=0.9,y=0,F(x,y)))-value(subs(x=-1,y=0,F(x,y)))));

An easier way to find the volume in this case is simply to integrate over the relevant region.

> value(Doubleint(f(x,y),x=-1..0.9,y=0..1));