src/monte_carlo/Metropolis_Gaussian_Demo.c

   1
   2 /*   Metropolis_Gaussian_Random.c                                   */
   3 /*   Demonstration of random numbers                                */
   4 /*   Metropolis Algorithm                                           */
   5 /*   Gaussian Random Numbers distribution                           */
   6 /*                                                                  */
   7
   8 /*   Syntax:                                                        */
   9 /*   Metropolis_Gaussian_Demo 10000 > data_file                     */
  10
  11                                                         /* Michel Vallieres */
  12
  13 #include <stdlib.h>
  14 #include <stdio.h>
  15 #include <math.h>
  16
  17
  18 /* probability profile of the random numbers */
  19 double weight( double x, double y )
  20 {
  21   return exp( - ( x*x + y*y ) );
  22 }
  23
  24
  25 /* Metropolis algorithm to generate Random       */
  26 /* Numbers in the plane with probability profile */
  27 /* weight(x,y)                                   */
  28 void metropolis ( double *x, double *y, double (*weight)(double,double), double delta )
  29 {
  30   double xt, yt, xx, yy, eta, ratio;
  31
  32   xx = (double)rand()/(double)RAND_MAX;          /* new random numbers 0 to 1 */
  33   yy = (double)rand()/(double)RAND_MAX;
  34   xt = *x + delta * ( 2*xx - 1.0 );              /* in range -1 to 1 */
  35   yt = *y + delta * ( 2*yy - 1.0 );
  36   ratio = weight( xt, yt ) / weight( *x, *y );   /* ration of probability profiles */
  37   eta = (double)rand()/(double)RAND_MAX;         /* random # in 0-1  */
  38   if ( ratio > eta )                             /* Metropolis selection */
  39     {
  40       *x = xt;                                   /* selected numbers */
  41       *y = yt;
  42     }
  43 }
  44
  45
  46 /* demo program to illustrate Gaussina random numbers in the plane */
  47 int main( int argc, char *argv[] )
  48 {
  49   int    i, count, N;
  50   double x, y, x_old, y_old, delta;
  51                                                  /* argument list -> N */
  52   if ( argc != 2 )
  53     {
  54       printf( "\n" );
  55       printf( "Syntax: Metropolis_Gaussian_Random  N \n" );
  56       exit(1);
  57     }
  58   N = atoi( argv[1] );
  59                                                 /*  Metropolis delta - size of jump */
  60   delta = 0.2;
  61   x = 0.5;                                      /* initila position in plane */
  62   y = 0.6;
  63   count = 0;                                    /* counter for selected numbers */
  64   printf( " %f %f\n", x, y );
  65
  66
  67   for ( i=0 ; i<N ; i++ )
  68     {
  69       x_old = x;
  70       y_old = y;
  71       metropolis ( &x, &y, weight, delta );     /* Metropolis Algorithm */
  72       if ( x != x_old && y != y_old )
  73           count++;
  74       printf( " %f %f\n", x, y );
  75     }
  76                                                 /* Acceptance ratio - should be > ~0.5 */
  77   fprintf( stderr, "\n Acceptance ratio: %f\n\n", (double)count/(double)N );
  78
  79 }