# Metropolis Monte Carlo Sampler in Rust

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!

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)]

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 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

• 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. – Boiethios Aug 31 '18 at 13:59
• Is the non snake case name? – Theodor Johnson Aug 31 '18 at 14:13
• Use rustfmt. The weirder are your brackets and the lack of spaces. – Boiethios Aug 31 '18 at 14:16
• Beside this, I do not see anything else in your code :) – Boiethios Aug 31 '18 at 14:21
• Hello! I installed rustfmt and also fixed the formatting to the best of my abilities! Thanks for the comment! – Theodor Johnson Aug 31 '18 at 14:36

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.

• Please heed How do I write a Good Answer?: your post provides no insight about the code presented for review. – greybeard Mar 31 at 10:17
• Thanks for the fixes! It has been appreciated even almost three years later! – Theodor Johnson Apr 2 at 13:16