gyak2_rejection_sampling.c

Go to the documentation of this file.

/*      gyak2_rejection_sampling.c
 *
 *      created by: racze
 *      date: 2011-09-25
 *
 *      This program is an example for rejection sampling.
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_rng.h>
void Generate_Corrected_Pareto(double* pointer, size_t length, double alpha, double delta, double xmin, const gsl_rng *rng);

int Hi2Lo(const void* a, const void* b);

int main(void){
        size_t N = 10000;
        size_t i;
        double data[N];

        /*Initialize for the GSL random number generator.*/
        const gsl_rng_type * T;
        gsl_rng * r;
     
        gsl_rng_env_setup();
        /*Choose rng type.*/
        T = gsl_rng_default;
        /*Allocate memory to the rng.*/
        r = gsl_rng_alloc (T);
        /*Set the seed of the random number generator.*/
        gsl_rng_set(r, 123456);

        Generate_Corrected_Pareto(data, N, 2.0, 1.0, 10.0, r);

/*      for(i = 0; i < N; ++i){
                fprintf(stdout, "%g\n" , data[i]);
        }*/
        
        /*Free the memory allocated to the rng.*/
        gsl_rng_free(r);

        /*
 *      HOMEWORK: Check whether rejection sampling actually works:
 *      1) analytically
 *      and/or
 *      2) by constructing the empirical distribution function of a larger sample (say N=1000)
 *      and comparing it to the analytic function
 */

        qsort(data, N, sizeof(double), Hi2Lo);

        for(i = 0; i < N; ++i){
                fprintf(stdout, "%g\t%g\n", data[i], ((double) i+1)/((double) N));
        }


        return 0;       
}


void Generate_Corrected_Pareto(double* data, size_t length, double alpha, double delta, double xmin, const gsl_rng *rng){
        size_t index = 0;
        short accepted;
        double x, u;
        while(index < length){
                accepted = 0;
                fprintf(stderr, "------------- index = %zu -------------", index);
                while(!accepted){
                        /*1) Generate a random number distributed according to (x/xmin)^(-alpha)*/
                        u = gsl_rng_uniform(rng);
                        x = xmin*pow(u, -1/alpha);
                        /*2) Generate another random number to decide whether 
 *                         to accept or reject the previously generated value*/
                        u = gsl_rng_uniform(rng);
                        accepted = (u <  1.0 - pow(x/xmin, -delta));
                        /*accepted = 1;*/
                        fprintf(stderr, "accepted = %d\n", accepted);
                }
                data[index++] = x;
        }
}

int Hi2Lo(const void* a, const void* b){
        double vala = *(double*) a;
        double valb = *(double*) b;

        if(vala > valb){
                return -1;
        }else if(vala < valb){
                return 1;
        }else{
                return 0;
        }
}