Random multivariate normal variates.

def rmv_normal_chol(mu, sig, size=1):
    """
    Random multivariate normal variates.
    """
    mu_size = np.shape(mu)

    if size == 1:
        out = np.random.normal(size=mu_size)
        try:
            flib.dtrmm_wrap(sig, out, 'L', 'N', 'L', 1.)
        except:
            out = np.dot(sig, out)
        out += mu
        return out
    else:
        if not hasattr(size, '__iter__'):
            size = (size,)
        tot_size = np.prod(size)
        out = np.random.normal(size=(tot_size,) + mu_size)
        for i in xrange(tot_size):
            try:
                flib.dtrmm_wrap(sig, out[i, :], 'L', 'N', 'L', 1.)
            except:
                out[i, :] = np.dot(sig, out[i,:])
            out[i, :] += mu
        return out.reshape(size + mu_size)