bivariate_normal_pdf := ( x_1, x_2, mu_1, sigma_1, mu_2, sigma_2, rho ) ->
  exp( 
       ( - 1 / ( 2 * ( 1 - rho^2 ) ) ) * 
         ( 
           ( ( x_1 - mu_1 ) / sigma_1 )^2
           - 2 * rho * ( ( x_1 - mu_1 ) / sigma_1 ) *
           ( ( x_2 - mu_2 ) / sigma_2 )
           + ( ( x_2 - mu_2 ) / sigma_2 )^2
         ) 
      ) / 
  ( 2 * Pi * sqrt( 1 - rho^2 ) * sigma_1 * sigma_2 );

# mu_1 is the mean of x_1, and mu_2 is the mean of x_2

# sigma_1 is the standard deviation (sd) of x_1, and
# sigma_2 is the standard deviation (sd) of x_2

# rho is the correlation between x_1 and x_2

plot3d( 
  bivariate_normal_pdf( x_1, x_2, 0, 1, 0, 1, 0 ),
  x_1 = -3 .. 3, x_2 = -3 .. 3 );

# x_1 and x_2 both have mean 0 and sd 1,
# and are uncorrelated (correlation 0)

plot3d( 
  bivariate_normal_pdf( x_1, x_2, 1, 1, -2, 1, 0 ),
  x_1 = -3 .. 3, x_2 = -3 .. 3 );

# here i've kept the correlation at 0 and the sds at 1,
# and set the mean of x_1 to +1 and the mean of x_2 to -2,
# but i forgot that this affects the range over which
# the density should be plotted

# we'll see later that ( mean +/- 3 * sd ) captures
# most of the probability for a normal random variable

plot3d( 
  bivariate_normal_pdf( x_1, x_2, 1, 1, -2, 1, 0 ),
  x_1 = -2 .. 4, x_2 = -5 .. 1 );

# now the plot is identical to the first plot above
# except that x_1 has been shifted to the right by 1 unit
# and x_2 has been shifted to the left by 2 units

# the means evidently locate the center of the distribution

plot3d( 
  bivariate_normal_pdf( x_1, x_2, 0, 2, 0, 3, 0 ),
  x_1 = -3 .. 3, x_2 = -3 .. 3 );

# here i've kept the correlation at 0 and the means at 0,
# and set the sd of x_1 to 2 and the sd of x_2 to 3,
# but i forgot again that this affects the range over which
# the density should be plotted

# the sd values affect the spread of the distribution

plot3d( 
  bivariate_normal_pdf( x_1, x_2, 0, 2, 0, 3, 0 ),
  x_1 = -6 .. 6, x_2 = -9 .. 9 );

# now i'll leave the means at 0 and the sd values at 1
# and change the correlation from 0 to +0.9

plot3d( 
  bivariate_normal_pdf( x_1, x_2, 0, 1, 0, 1, 0.9 ),
  x_1 = -3 .. 3, x_2 = -3 .. 3 );

# the larger the absolute value of rho is, the easier
# it becomes to predict one variable from another

# here's another way to visualize 2-dimensional distributions:

plot3d( 
  bivariate_normal_pdf( x_1, x_2, 0, 1, 0, 1, 0 ),
  x_1 = -3 .. 3, x_2 = -3 .. 3 );

with( plots );

contourplot( 
  bivariate_normal_pdf( x_1, x_2, 0, 1, 0, 1, 0 ),
  x_1 = -3 .. 3, x_2 = -3 .. 3 );

contourplot3d( 
  bivariate_normal_pdf( x_1, x_2, 0, 1, 0, 1, 0 ),
  x_1 = -3 .. 3, x_2 = -3 .. 3 );

contourplot( 
  bivariate_normal_pdf( x_1, x_2, 0, 1, 0, 1, 0.9 ),
  x_1 = -3 .. 3, x_2 = -3 .. 3 );

contourplot3d( 
  bivariate_normal_pdf( x_1, x_2, 0, 1, 0, 1, 0.9 ),
  x_1 = -3 .. 3, x_2 = -3 .. 3 );