7
\$\begingroup\$

the following is an implementation of the standard Metropolis Hastings Monte Carlo sampler. You can read more about it here.

At the end I am going to give you a link to the Rust playground, so you can test the code yourself!

First I start with some traits, that define probability functions:

type Float=f64;
/*Standard trait for distributions. Some applications only
require the ability to sample from a distribution without the need
for an explicit probability density function (pdf).*/
pub trait Distribution
{
    type T;
    fn sample<R:Rng+?Sized>(&self, rng:&mut R)->Self::T;
}

/*Standard trait for distributions. Some applications only
require the ability to sample from a distribution without the need
for an explicit probability density function (pdf).*/
pub trait PDF: Distribution
{
    fn pdf(&self,x:&Self::T)->Float;
}

/*When Rust offers optional function arguments, the conditional 
traits can be depraciated.*/
pub trait ConditionalDistribution
{
    type T;
    fn csample<R:Rng+?Sized>(&self, rng:&mut R,x_prev:&Self::T)->Self::T;
}

/*Conditional probability density function. This is especially useful, if the
expression for the conditional pdf does only depend on the shape.*/
pub trait ConditionalPDF: ConditionalDistribution
{
    fn cpdf(&self,x:&Self::T,x_prev:&Self::T)->Float;
}

The trait declarations are modelled after the ones in the rand crate.

Then comes the sampler itself:

#[doc = "Markov Chain Metropolis Hastings algorithm. Numerically evaluate the value of an integral, where the integrand is proportional to
a probability distribution."]

#[allow(non_snake_case)]
#[allow(dead_code)]

struct MetropolisHastings<D,F>
where 
    D:ConditionalDistribution+ConditionalPDF,
    D::T:AsRef<[Float]>+std::fmt::Debug+Clone,
    F:Fn(&[Float])->Float
{
    log_f:F,
    log_proposal:D,
    initial:D::T,
    burnIn:usize,
    isSymmetric:bool,
}

impl<D,F> MetropolisHastings<D,F>
where 
    D:ConditionalDistribution+ConditionalPDF,
    D::T:AsRef<[Float]>+std::fmt::Debug+Clone,
    F:Fn(&[Float])->Float
{
    fn new(log_f:F,log_proposal:D,x0:D::T)->MetropolisHastings<D,F>
    {
        MetropolisHastings{
        log_f,
        log_proposal,
        initial:x0,
        burnIn:(250 as Float) as usize,
        isSymmetric:false,
    }

}

    fn sample(&self, n:usize)->Vec<Float>
    {
        let dim=1;//change!
        let mut samples=Vec::with_capacity(n*dim);
        let mut rng=rand::thread_rng();
        let mut x=self.initial.clone();

        for i in 0..self.burnIn+n as usize{
            let x_p=self.log_proposal.csample(&mut rng,&x);
            let acc_rat=(self.log_f)(x_p.as_ref())+self.log_proposal.cpdf(&x_p,&x)-(self.log_f)(x.as_ref())-self.log_proposal.cpdf(&x,&x_p);
            let r=f64::min(0.0,acc_rat);
            let u=rng.gen::<Float>();

            //the polynomial is a truncated Taylor series to u.ln(x) at x=0.5             
            if -0.7*2.0*(u-0.5)-2.0*(u-0.5)*(u-0.5) <r
            { 
                    x=x_p.clone();
            }

            if i>=self.burnIn
            {
                samples.extend_from_slice(x_p.as_ref());
            } 
        } 
        samples
    }
}

It is required, that the user inputs the densities in log scale. Because Rust's f64::ln() function takes too long, I approximated this function in the intervall [0,1] with a truncated Taylor series at x=0.5. If the input is multidimensional, the output will be a flattened vector of size numSamples*dimension.

Questions: With a burn-in of 250 and 100.000.000 million samples, it takes 2.3 seconds on my machine. Questions:

  • Can I make the code any faster or better?
  • Do you like my traits?
  • Can I infer the dimension of the input automatically instead of fixing dim=1 in fn sample()?

Here you can execute the bare code yourself:

https://play.rust-lang.org/?gist=9e9ac875d3b6071d9b873b5384bef8a8&version=stable&mode=debug&edition=2015

Here you can execute the code with an added example: https://play.rust-lang.org/?gist=c37ee575716e37e1668b5a6f62f8e8cd&version=stable&mode=debug&edition=2015

\$\endgroup\$
5
  • 1
    \$\begingroup\$ Don't take it personally, but your code formatting is ugly :( You should use the official Rust formatting if you want to share your code. \$\endgroup\$ – Boiethios Aug 31 '18 at 13:59
  • \$\begingroup\$ Is the non snake case name? \$\endgroup\$ – Theodor Johnson Aug 31 '18 at 14:13
  • \$\begingroup\$ Use rustfmt. The weirder are your brackets and the lack of spaces. \$\endgroup\$ – Boiethios Aug 31 '18 at 14:16
  • \$\begingroup\$ Beside this, I do not see anything else in your code :) \$\endgroup\$ – Boiethios Aug 31 '18 at 14:21
  • 3
    \$\begingroup\$ Hello! I installed rustfmt and also fixed the formatting to the best of my abilities! Thanks for the comment! \$\endgroup\$ – Theodor Johnson Aug 31 '18 at 14:36
0
\$\begingroup\$

Wow, this is old. I've gone through it and made some very small format changes and changes to make it compile/build/run; See playground.

Highlights:

Documentation commments can be added using attributes

#[doc = "Markov Chain Metropolis Hastings algorithm. Numerically evaluate the value of an integral, where the integrand is proportional to
a probability distribution."]

or using ///-comment blocks:

/// Markov Chain Metropolis Hastings algorithm. Numerically evaluate the
/// value of an integral, where the integrand is proportional to
/// a probability distribution.
\$\endgroup\$
2
  • 4
    \$\begingroup\$ Please heed How do I write a Good Answer?: your post provides no insight about the code presented for review. \$\endgroup\$ – greybeard Mar 31 at 10:17
  • 1
    \$\begingroup\$ Thanks for the fixes! It has been appreciated even almost three years later! \$\endgroup\$ – Theodor Johnson Apr 2 at 13:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.