MetropolisHastings object

So the rationale for creating this is as follows. There is a function `MarkovChain.metropolisHastings` which does a very basic MarkovChain sampler. It doesn't incorporate important features like burnin (to skip a bunch of initial samples, which do not represent your posterior) or dropping a chunk of the samples (to reduce autocorrelation). 

The benchmarks also show that `ArbitraryMetropolisHastings` is 20% faster than the `MarkovChain.metropolisHastings` (provided we set burnin=0 and dropCount=0). This is probably due to the fact that:

1. If `a > 1`, I don't call Uniform(0,1).draw()
2. The monadic way of writing code creates a lot of anonymous function objects, which can hurt performance.

Please go through the code **carefully** since I don't know how to write unit tests for an MCMC algorithm. 

I've also created a variation which puts the generation work into a background thread. This can be handy to either reduce latency, or to use multiple CPUs in the event that a large amount of work needs to be done after the sample is generated. 
No. You have to be able to also compute the probability/density of
proposing the sample x' you just proposed, given the old sample x, and vice
versa. It's not enough to be able to sample.

-- David


On Wed, Jul 16, 2014 at 3:59 AM, Chris Stucchio <notifications@github.com>
wrote:

> In math/src/main/scala/breeze/stats/mcmc/MetropolisHastings.scala:
>
> > +
> > +trait TracksStatistics { self:MetropolisHastings[_] =>
> > +  /* This obviously can't track stats for you, so you need to build it in.
> > +   * Tracking stats appears to have minimal performance implication, probably
> > +   * because incrementing a few longs is far cheaper than generating random
> > +   * numbers.
> > +   */
> > +  def total: Long
> > +  def acceptanceCount: Long
> > +  def aboveOneCount: Long
> > +  def rejectionCount: Long = total - acceptanceCount
> > +  def aboveOneFrac: Double = aboveOneCount.toDouble / total.toDouble
> > +  def rejectionFrac: Double = rejectionCount.toDouble / total.toDouble
> > +}
> > +
> > +abstract class BaseMetropolisHastings[T](val logLikelihood: T => Double, init: T, burnIn: Long = 0, dropCount: Int = 0)(implicit val rand:RandBasis=Rand) extends MetropolisHastings[T] with Process[T] with TracksStatistics {
>
> That's proposalDraw. Not implemented on the abstract base class,
> downstream classes implement it. They can do it in the slow way x =>
> rnd(x).draw() or via some faster method with no object creation.
>
> —
> Reply to this email directly or view it on GitHub
> <https://github.com/scalanlp/breeze/pull/277/files#r14990895>.
>