normalwishartdistribution.py 文件源码

def wishartrand(self, nu, sigma, 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] = np.sqrt(random.gammavariate(0.5*(nu - D + 1), 2.0))
        return np.dot(np.dot(np.dot(C, a), a.T), C.T)