Skip to content

Instantly share code, notes, and snippets.

@jfrelinger
Created May 8, 2012 18:55
Show Gist options
  • Save jfrelinger/2638485 to your computer and use it in GitHub Desktop.
Save jfrelinger/2638485 to your computer and use it in GitHub Desktop.
wishart and inverse wishart sampler
import numpy as np
import numpy.random as npr
from numpy.linalg import inv, cholesky
from scipy.stats import chi2
def invwishartrand_prec(nu,phi):
return inv(wishartrand(nu,phi))
def invwishartrand(nu, phi):
return inv(wishartrand(nu, inv(phi)))
def wishartrand(nu, phi):
dim = phi.shape[0]
chol = cholesky(phi)
#nu = nu+dim - 1
#nu = nu + 1 - np.arange(1,dim+1)
foo = np.zeros((dim,dim))
for i in range(dim):
for j in range(i+1):
if i == j:
foo[i,j] = np.sqrt(chi2.rvs(nu-(i+1)+1))
else:
foo[i,j] = npr.normal(0,1)
return np.dot(chol, np.dot(foo, np.dot(foo.T, chol.T)))
if __name__ == '__main__':
npr.seed(1)
nu = 5
a = np.array([[1,0.5,0],[0.5,1,0],[0,0,1]])
#print invwishartrand(nu,a)
x = np.array([ invwishartrand(nu,a) for i in range(20000)])
nux = np.array([invwishartrand_prec(nu,a) for i in range(20000)])
print x.shape
print np.mean(x,0),"\n", inv(np.mean(nux,0))
#print inv(a)/(nu-a.shape[0]-1)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment