2 /* Metropolis_Gaussian_Random.c */
3 /* Demonstration of random numbers */
4 /* Metropolis Algorithm */
5 /* Gaussian Random Numbers distribution */
9 /* Metropolis_Gaussian_Demo 10000 > data_file */
11 /* Michel Vallieres */
18 /* probability profile of the random numbers */
19 double weight( double x, double y )
21 return exp( - ( x*x + y*y ) );
25 /* Metropolis algorithm to generate Random */
26 /* Numbers in the plane with probability profile */
28 void metropolis ( double *x, double *y, double (*weight)(double,double), double delta )
30 double xt, yt, xx, yy, eta, ratio;
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 */
40 *x = xt; /* selected numbers */
46 /* demo program to illustrate Gaussina random numbers in the plane */
47 int main( int argc, char *argv[] )
50 double x, y, x_old, y_old, delta;
51 /* argument list -> N */
55 printf( "Syntax: Metropolis_Gaussian_Random N \n" );
59 /* Metropolis delta - size of jump */
61 x = 0.5; /* initila position in plane */
63 count = 0; /* counter for selected numbers */
64 printf( " %f %f\n", x, y );
67 for ( i=0 ; i<N ; i++ )
71 metropolis ( &x, &y, weight, delta ); /* Metropolis Algorithm */
72 if ( x != x_old && y != y_old )
74 printf( " %f %f\n", x, y );
76 /* Acceptance ratio - should be > ~0.5 */
77 fprintf( stderr, "\n Acceptance ratio: %f\n\n", (double)count/(double)N );