Compute random values distributed according to the specified bivariate
    distribution.

    Inputs:

    * x0: the center of the x distribution (i.e. its intended mean)
    * y0: the center of the y distribution
    * sx: standard deviation (not variance) of x var
    * sy: standard deviation (not variance) of y var
    * cxy: covariance (not correlation coefficient) of x and y
    * size (optional): the number of values to compute

    Returns: array of shape (size, 2); or just (2, ), if size was not
      specified.

    The bivariate parameters of the generated data are approximately
    recoverable by calling 'databiv(retval)'.

def bivrandom (x0, y0, sx, sy, cxy, size=None):
    """Compute random values distributed according to the specified bivariate
    distribution.

    Inputs:

    * x0: the center of the x distribution (i.e. its intended mean)
    * y0: the center of the y distribution
    * sx: standard deviation (not variance) of x var
    * sy: standard deviation (not variance) of y var
    * cxy: covariance (not correlation coefficient) of x and y
    * size (optional): the number of values to compute

    Returns: array of shape (size, 2); or just (2, ), if size was not
      specified.

    The bivariate parameters of the generated data are approximately
    recoverable by calling 'databiv(retval)'.

    """
    from numpy.random import multivariate_normal as mvn
    p0 = np.asarray ([x0, y0])
    cov = np.asarray ([[sx**2, cxy],
                       [cxy, sy**2]])
    return mvn (p0, cov, size)