def bivariate_normal(X, Y, sigmax=1.0, sigmay=1.0, mux=0.0, muy=0.0, sigmaxy=0.0):
 Xmu = (X - mux)
 Ymu = (Y - muy)
 rho = (sigmaxy / (sigmax * sigmay))
 z = ((((Xmu ** 2) / (sigmax ** 2)) + ((Ymu ** 2) / (sigmay ** 2))) - ((((2 * rho) * Xmu) * Ymu) / (sigmax * sigmay)))
 denom = ((((2 * np.pi) * sigmax) * sigmay) * np.sqrt((1 - (rho ** 2))))
 return (np.exp(((- z) / (2 * (1 - (rho ** 2))))) / denom)

def bivariate_normal ( X , Y , sigmax = 1.0 , sigmay = 1.0 , mux = 0.0 , muy = 0.0 , sigmaxy = 0.0 ) : Xmu = ( X - mux ) Ymu = ( Y - muy ) rho = ( sigmaxy / ( sigmax * sigmay ) ) z = ( ( ( ( Xmu ** 2 ) / ( sigmax ** 2 ) ) + ( ( Ymu ** 2 ) / ( sigmay ** 2 ) ) ) - ( ( ( ( 2 * rho ) * Xmu ) * Ymu ) / ( sigmax * sigmay ) ) ) denom = ( ( ( ( 2 * np . pi ) * sigmax ) * sigmay ) * np . sqrt ( ( 1 - ( rho ** 2 ) ) ) ) return ( np . exp ( ( ( - z ) / ( 2 * ( 1 - ( rho ** 2 ) ) ) ) ) / denom )
bivariate gaussian distribution for equal shape *x* .