# Bivariate distribution ... graphing example

Blog Post created by Dan_Patterson on Nov 12, 2016

A stats ditty... so I don't forget and some of you may be interested in graphics using MatPlotLib

Code for the above...written verbosely so you get the drift.

``"""Script:    bivariate_normal.pyPath:      F:\A0_Main\Author:    Dan.Patterson@carleton.caCreated:   2015-06-06Modified:  2015-06-06  (last change date)Purpose:   To examine the affect of parameters on the bivariate distributionRequires:  numpy and matplot libNotes:   see help on (np.random.normal)  x = numpy.array([1, 2, 3])         # put in the values for the x values  y = numpy.array([10, 20, 30])   # ditto for y  XX, YY = numpy.meshgrid(x, y)   # make your 'mesh' which is the result  ZZ = XX + YY                    # in this case the sum of X and Y  ZZ => array([[11, 12, 13],               [21, 22, 23],               [31, 32, 33]])  >>> (X,Y) = meshgrid(x,y)     # actually Y  YY, XX = numpy.mgrid[10:40:10, 1:4]  # Y,X  ZZ = XX + YY # These are equivalent to the output of meshgrid  YY, XX = numpy.ogrid[10:40:10, 1:4] #  ZZ = XX + YY # These are equivalent to the atleast_2d example  which  XX, YY = numpy.atleast_2d(x, y)  YY = YY.T # transpose to allow broadcasting  ZZ = XX + YYReferences: many"""import numpy as npimport matplotlib.pyplot as pltnp.set_printoptions(edgeitems=3,linewidth=75,precision=2,                    suppress=True,threshold=10)# (1) make a floating point grid and get some stats#     100x100 grid cells numbered from top leftXs = np.arange(0., 100.1, 1) # as floatsYs = np.arange(0., 100.1, 1)Xc = Xs.mean();   Yc = Ys.mean()Xmin = Xs.min();  Xmax = Xs.max()Ymin = Ys.min();  Ymax = Ys.max()# (2) ....now do some work.... -----------------------------------------|X,Y = np.meshgrid(Xs,Ys)   # X, Y as a meshgridXX = X**2                  # square the Xout = 2*X + Y - 3.0        # do the spatial mathZ = X**2 + Y**2            # getting it now?        # (3) .... calculate the pdf -------------------------------------------|#   normal pdf/ellipserho = -0.5s_x = 60.0s_y = 45.0   # 4*3 ratio Xc,Yc = (50,50) stds at 2Std leveld_x = ((X-Xc)/s_x)d_y = ((Y-Yc)/s_y)rho2 = (1.0-rho**2)m_rsq = np.sqrt(rho2)lower = 2.0*(1-rho**2)                    # rho = 0 lower = 2upper = d_x**2 + d_y**2 - 2.0*rho*d_x*d_y # if rho= 0 drop last termleft = (1.0/(2.0*np.pi*s_x*s_y*m_rsq))f_xy = left**(upper/lower)# (4) .... make a figure -----------------------------------------------|fig = plt.figure(facecolor='white')ax = fig.add_subplot(1, 1, 1)plt.axis('equal')plt.axis([Xmin, Xmax, Ymax, Ymin]) # decreasing Yplt.set_cmap('Blues')cont = plt.contourf(Xs, Ys, f_xy, origin='upper')plt.title("Bivariate normal distribution")plt.xlabel("X ==>")plt.ylabel("<== Y")cbar = plt.colorbar(cont)cbar.ax.set_ylabel('values')lbl_r = "rho = {}\n2 std x\n1 std y".format(abs(rho))  # reverse rho since plotting -ve Yplt.text(0,20,lbl_r)#plt.show()plt.close()``

Check out their code gallery for a multitude of options on the MatPlotLib Home Page.