mc_integrate_metropolis.cc

/************** A simple MC integrator **************/
#include<iostream>
#include<boost/random.hpp>

/************************************************
 In order to calculate $\int f(x) dx$ we write 
(all integrals are from 0 to $\infty$, but this doesn't matter much)
$$
  \int f(x) dx = \int w(x)dx \frac{ \int f/w w(x)dx }{ \int w(x)dx } \;.
$$

Now, for the latter fraction we use the Metropolis trick:
organize a random walk {x_\nu} using w(x) as an (unormalized!)
probability density. Then the stochastic sum
$$
  \frac{ \sum_\nu f(x_\nu) / w(x_\nu) }{ \sum_\nu 1 }
$$
converges to 
$$
  \frac{ \int f/w w(x)dx }{ \int w(x)dx }
$$
************************************************/

using namespace std ; 

double f(double);  // to be integrated
double w(double);  // probab. dens.
double int_w();    // \int w(x)dx

int main(){

  boost::mt19937 rng(11u); 
  boost::variate_generator< boost::mt19937&, boost::uniform_01<> > rndm(rng, boost::uniform_01<>()) ;

  double dt=0.01;
  double x=1.;  // "cnf"

  double sum=0., Z=0. ;
  size_t sweep=0;  
  for(;;){
   sweep++;
   for(int j=0; j<5000000 ; ++j ){
     //update
     double xnew=x+dt*(2.*rndm()-1);
     if(xnew<0)continue;
     double ratio = w(xnew)/w(x);
     // metropolis
     if( ratio >1. || rndm() <ratio ){ x=xnew;  }
     // measure
       sum += f(x)/w(x) ; 
       Z += 1.;   
   }   
   cout<<sweep<<":  "<<int_w()*sum/Z<<"\n";
  
   if(sweep==50){
     std::cout<<"************ therm done\n";
     sum=0.;
     Z=0;
   }

   // "print cnf"
   //cout<<x<<"\n";
  }  

}


// to be integrated
double f(double x){
  return exp(-1.*x);
}

// probab density
double w(double x){
  return exp(-2.*x);
}
double int_w(){
   return 0.5;
}