Efficient sampling from mixture models

The changes in this pull request are the following

####  Metropolis Hastings Step Method ===
A new step method, called MetropolisHastings which allows arbitrary (non-symmetric) proposal distributions. Implemented wrappers for the existing Proposals as well as explicit support for Categorical proposal distributions. This can easily be extended to cover more proposals (how about HamiltonianMC Proposals, or Swendsen-Wang ?). 

The implementation could be made more efficient (some potentially complex terms get added to the acceptance logp, and then subtracted again later on), but it's hard to do that in a general manner. Maybe one could teach Theano to optimize this out. 
 
####  Efficient sampling from Mixture Models 
Using that Metropolis Hastings step method, it's easy to create an efficient sampler for Mixture models. So I wrote a first implementation for a Multivariate Normal Mixture distribution. Seems to work well. This could easily be extended to allow flexible mixture kernels.

#### Small Examples
For both, the new Step Method and the Mixture Model, I provided examples. 

#### Small fixes

I noticed that the "bound" function from distributions/dist_math.py fails to work correctly in some cases. Especially in the logp of the Categorical distribution, it failed not to evaluate at disallowed values, which
could then lead to index out of bounds exceptions. I fixed that by introducing a new "bound_scalar" function
which uses theano.ifelse.ifelse(..), which has the advantage over theano.switch(..) of not evaluating inactive parts of the expression
graph.

Another small fix had to do with TransformedDistributions. In one example I used a transformed Dirichlet distribution where one of the hyperparameters was a 1. This led to an invalid mode for the distribution. That's no problem, unless you transform that Dirichlet, in which case the mode becomes the starting point for the sampler. Since it's invalid (contains NaN values), that leads to problems.

#### Potential uses
The GMMs can be used as flexible prior densities. They can approximate any function.

Another use case for efficient GMMs is to use them to represent posterior densities. Since it's possible to
efficiently sample from products of mixtures of gaussians, these can be used to enable distributed scalable MCMC, both in the form of embarassingly parallel MCMC, as well as in distributed Nonparametric Belief Propagation (NBP) inference algorithms. I plan to write more on that later. See also my PyPGMc project on GitHub.
Hi,

All of these changes sound reasonable, I will change the code accordingly
once I find some time. Here are a few more things I will probably implement
next:

1.) Add a method to Model to calculate the logp of a specific set of
variables (or to exclude a specific set of variables)
2.) Allow GenericProposal instances to indicate whether they take care of
"their" variables, i.e that they sample from them exactly. In that case,
the Metropolis Hastings Sampler doesn't need to evaluate the logp of these
variables. Evaluating them can be expensive for MvGaussianMixture, since
one way or another, either the sampling or the logp evaluation involves
quite a few expensive calculations, like a Matrix inversion etc.
3.) Make the MvGaussianMixture parameterized via covariance matrices
instead of precision, since covariance is what's required for sampling (and
thus reduces the computational effort). Precision matrices won't ever have
to enter the stage this way.
4.) Implement an InverseWishart Distribution (so we have a proper prior for
these covariance matrices)

any thoughts on these ?

best,

Kai







2014-06-09 3:16 GMT+02:00 Chris Fonnesbeck <notifications@github.com>:

> In pymc/step_methods/metropolis_hastings.py:
>
> > +
> > +        return newpoint
> > +
> > +    def _ensure_tuned(self):
> > +        if self.tune and not self.steps_until_tune:
> > +            # Tune scaling parameter
> > +            acc_rate = self.accepted / float(self.tune_interval)
> > +            self.acceptance_rate.append(acc_rate)
> > +            self.scaling = mh_tune(self.scaling, acc_rate)
> > +            for p in self.proposals:
> > +                p.tune_stepsize(self.scaling)
> > +            # Reset counter
> > +            self.steps_until_tune = self.tune_interval
> > +            self.accepted = 0
> > +
> > +def mh_tune(scale, acc_rate):
>
> Can we not use metropolis.tune here instead? No need to copy and paste
> identical code.
>
> —
> Reply to this email directly or view it on GitHub
> <https://github.com/pymc-devs/pymc/pull/547/files#r13527552>.
>
>