wishart.py 文件源码

def wishrnd(sigma, v_0, C=None):
    """Return a sample from a Wishart distribution."""
    if C == None:
        C = np.linalg.cholesky(sigma)
    D = sigma.shape[0]
    a = np.zeros((D, D), dtype=np.float32)
    for r in xrange(D):
        if r != 0:
            a[r, :r] = np.random.normal(size=(r,))
        a[r, r] = math.sqrt(random.gammavariate(0.5*(v_0 - D + 1), 2.0))
    return np.dot(np.dot(np.dot(C, a), a.T), C.T)