[SciPy-user] Monte Carlo package

Robert Kern robert.kern at gmail.com
Wed Jan 25 13:38:49 CST 2006


David Huard wrote:
> I mean an arbitrary user-defined distribution. I want to generate
> samples using a Gibbs sampler, and one of the conditional distribution
> has a weird shape involving a sum of logs... Anyway, there is no chance
> that I can find an analytical expression for this cdf and I wondered if
> there was a routine somewhere that would let me define the function and
> would return samples drawn from this distribution. I thought about that
> since the intsampler of montecarlo does something similar for discrete
> distributions.

Well, Gibbs sampling is only useful when the conditional distributions can be
sampled quickly. Otherwise, I would recommend using a full Metropolis-Hastings
algorithm with Gibbs-type jumps for the conditionals that you can sample quickly.

Now, there are semi-black-box algorithms for sampling nearly-arbitrary
univariate distributions, and I do intend to get around to implementing them for
scipy someday. For your problem, I'm not sure they would be more useful than
simply doing M-H since the parameters of your conditional will be changing each
time. _Automatic Nonuniform Random Variate Generation_ is an excellent book
which covers these algorithms:

  http://statistik.wu-wien.ac.at/projects/arvag/monograph/index.html

> I realize that if it was so simple, we wouldn't need MCMC algorithms...
> Anyway, I think I'll try to tweak PyMC to mix Gibbs sampling with
> Metropolis jumps. I guess that's possible, isn't it ?

Gibbs sampling is simply a special case of M-H sampling with an acceptance of 1,
so yes, I think so.

-- 
Robert Kern
robert.kern at gmail.com

"In the fields of hell where the grass grows high
 Are the graves of dreams allowed to die."
  -- Richard Harter



More information about the SciPy-user mailing list