8
\$\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, 2018 at 13:59
  • \$\begingroup\$ Is the non snake case name? \$\endgroup\$ Aug 31, 2018 at 14:13
  • \$\begingroup\$ Use rustfmt. The weirder are your brackets and the lack of spaces. \$\endgroup\$
    – Boiethios
    Aug 31, 2018 at 14:16
  • \$\begingroup\$ Beside this, I do not see anything else in your code :) \$\endgroup\$
    – Boiethios
    Aug 31, 2018 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\$ Aug 31, 2018 at 14:36

1 Answer 1

-2
\$\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, 2021 at 10:17
  • 1
    \$\begingroup\$ Thanks for the fixes! It has been appreciated even almost three years later! \$\endgroup\$ Apr 2, 2021 at 13:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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