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\$
  • 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\$ – French 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\$ – French Boiethios Aug 31 '18 at 14:16
  • \$\begingroup\$ Beside this, I do not see anything else in your code :) \$\endgroup\$ – French 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

Your Answer

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

Browse other questions tagged or ask your own question.