possible bug in wishart sampler

[GoogleCodeExporter]

I randomly came across your code as I was trying to confirm the correctness of 
my own Wishart sampler, and I think you might have a bug. Instead of:

  random.gammavariate(0.5*(self.dof-d+1),2.0)

I think it should be:

  random.gammavariate(0.5*(self.dof-r),2.0)

Anyway, sorry to bother you, and definitely let me know if I'm wrong. :) Thanks.

Original issue reported on code.google.com by amcna...@gmail.com on 2 May 2013 at 7:05

[GoogleCodeExporter]

Hi,
I also believe that this is a bug (but would like some expert confirmation).

For one, suggested fix corresponds to other sources, such as the Julia 
implementation, 
https://github.com/JuliaStats/Distributions.jl/blob/master/src/matrix/wishart.jl
 , and the R implementation, 
https://svn.r-project.org/R/trunk/src/library/stats/src/rWishart.c
(and Wikipedia, but unfortunately I don't have access to [Smith,1972] to check).

The expected value of samples from the Wishart W ~ Wish(S,d) should be E[W] = 
S*d
Without the suggested fix, this test case however does not work:

# ----------------
import numpy
import gcp.wishart
# create a sample mean
d = 10
Ss = numpy.random.rand(d,2)
S = Ss.T.dot(Ss)

W = gcp.wishart.Wishart(2)
W.setDof(d)
W.setScale(S)

# should give approximately the same output:
print numpy.array([W.sample() for _ in xrange(10000)]).mean(axis=0)
print S * d
# ----------------

Thank you

Original comment by julian.k...@gmail.com on 13 Aug 2013 at 10:30

[GoogleCodeExporter]

Here's the Smith 1972 paper. I hope you find it helpful.

Original comment by amcna...@gmail.com on 13 Aug 2013 at 10:32

Attachments:

• smith-1972-wishart-variate-generator.pdf

[GoogleCodeExporter]

Thank you, that seems to confirm it: diagonal value A_{i,i} should be sampled 
from Xi^2_{dof-i+1}, where 1 <= i <= d is using one-based indexing (although 
now I am unsure why that paper doesn't say to take the square root of those 
samples even though that seems correct ... but it is late now).

Original comment by julian.k...@gmail.com on 13 Aug 2013 at 11:37