Gibbs sampling parameters of Normal-InverseWishart distribution, does NIW have a prior itself?

I'll try to lend what experience I have, but I can't promise it'll be what you need. My background in this is in RBMs, so if that's not what you're tinkering with I hope this at least makes sense.

If you haven't already, read through the paper Exponential Family Harmoniums with an Application to Information Retrieval by Welling, Rosen-Zvi, and Hinton. It demonstrates that an RBM can be used with any exponential family distribution.

You want to try to construct an energy function that, when solved for P(v_i|h) or P(h_j|v), the two are either exponential family, or can be otherwise demonstrated to work if it isn't in the exponential family (eg, sigmoid, tanh, RLUs).

I haven't done the math for the Wishart, Inverse Wishart, or Normal Inverse Wishart distributions, but I have worked through designing an RBM that can use the Normal, Bernoulli, Binomial, Multinomial, Negative Binomial, Poisson, and while I haven't been successful I've at least attempted Pareto, Geometric, Exponential, Gamma, Chi-Squared, Beta, and Weibull.

The ones I succeeded on are 'easy' in the sense that I didn't run into any significant hurdles finding a representation that seemed useful with little in the way of fudging. The ones I didn't succeed on usually had restricted support or other complicating factors with their parameters.

For example, the exponential distribution only has support on [0,inf], and the inverse scale parameter (lambda) is required to be positive as well. Since the exponential distribution has only a single parameter, it's the only parameter that could contain the linear combinations W*tr(h) or v*W that would ultimately be part of the energy function, but any linear combination could potentially violate the constraint lambda>0 -- which could make things annoying when it comes time to sample from the distribution.

There are ways to get around that sort of problem, but by the time I got to this point I felt I had learned enough about the process (and I had no need to go on), so I moved on to other projects.

You'll run into similar issues with the Inverse-Wishart. I'm not at all familiar with that distribution, but from the wikipedia page on it it I'm going to assume the univariate version of the IW distribution needs parameters psi>0 and m>0 (since psi is p x p, and it's univariate, so p=1 => m>1-1). This situation is similar to the Pareto distribution in that all the parameters are restricted, and the support is restricted.

With all of that out of the way (my apologies for the verbosity), this is essentially the process I went through to construct an energy function and derive conditional distributions for doing Gibbs sampling:

• Construct a generic energy function of the form:
• E(v,h) = EVH(v,h) + EV(v) + EH(h)
• Supposing P(v|h) is the target distribution, EH(h) can be completely disregarded -- it will fall out while solving for P(v_i|h). Further, since the weight matrix is meant to link v to h, EV(v) can only depend on the visible bias -- though you can't disregard v just yet (that ends up being the problem with using the Poisson distribution; when sampling v\_i, it requires knowledge of v, hence Salakhutdinov & Hinton 'cheated' by breaking the conditional independence assumption).
• Minus the EH(h) term, E(v,h) will take a form that is nearly identical to ln( P(v_i|h) ) when you're through -- assuming P(v_i|h) is a true exponential-family distribution.
• It then becomes a matter of choosing a form for EVH and EV that, when you start working out P(v|h) = P(v,h) / P(h), will make it easy to reduce it from P(v|h) to P(v_i|h) and additionally have P(v_i|h) be an exponential-family distribution.

If you want some more details, when I get home I can post a small example of how I came up with an energy function for using the Poisson distribution.

Lastly: take a look at Alexander Krizhevsky's non-thesis masters thesis Learning Multiple Layers of Features from Tiny Images. I got some very valuable information about the process I outlined above, both from his paper as well as from email conversations. He doesn't outline the energy function derivation -- that he just states outright. But he does derive the conditionals, and explains that his training procedure required him to learn two sets of parameters -- the weights/biases, and a standard deviation for each visible unit.

I don't recall if his paper says he had to train them in separate phases as Yann LeCunn recommends doing when optimizing for multiple goals or constraints, but that's the impression I have -- you optimize the weights/bias holding the standard deviations constant, then hold those constant and optimize the standard deviations; rinse, repeat.

With the number of parameters the NIW has, I'm guessing you'll find yourself in a similar situation.