NAME

ransampl - random number sampling

SYNOPSIS

#include <ransampl.h>

ransampl_ws* ransampl_alloc( int n, double(*rng)() );

void ransampl_set( ransampl_ws *ws, double *p );

int ransampl_draw( ransampl_ws *ws );

void ransampl_free( ransampl_ws *ws );

DESCRIPTION

These routines allow to draw a random index according to given probabilities p[0], .., p[n-1].

The implementation uses the alias method of Walker and Vose in a form given by Schwarz. Using precomputed tables, the cost for a single drawing is only O(1). Therefore this is the recommended method for M drawings if M>>n>>1.

To prepare, a workspace must be allocated by ransampl_alloc, and some tables must be precomputed by ransampl_set. The supplied probabilities p are not requested to be normalized. Along with the workspace allocation, a floating-point valued uniform random-number generator function rng must be supplied.

Random samples are then drawn by the function ransampl_draw.

When done with the drawing, ransampl_free deallocates the workspace.

EXAMPLE

Draw representative inhabitants according to the given population numbers of nine federal states of Austria. This example is also contained in the source distribution, demo/sampling1.c.

    int main()
    {
        const int M=1000000;
        int i, m;

        // Discrete probability distribution example:
        const int n = 9;
        // states of Austria
        const char* names[] = {
            "Wien", "Niederoesterreich", "Oberoesterreich",
            "Steiermark", "Tirol", "Kaernten",
            "Salzburg", "Vorarlberg", "Burgenland" };
        // inhabitants in millions [www.statistik.at, 2011 data]
        double p[] = { 1.721573, 1.614661, 1.415020, 1.211506, .711161,
                       .558056, .532713, .370833, .285377 };

        // Allocate workspace, and precompute tables:
        printf( "Precomputing tables ...\n" );
        ransampl_ws* ws = ransampl_alloc( n, &random_number_generator );
           // To make this work, either define and implement
           // 'double random_number_generator()', or replace the pointer
           // &random_number_generator by a pointer to whatever uniform
           // random number generator function, typically from a library.
        ransampl_set( ws, p );

    #ifdef RANSAMPL_SHALL_REVEAL_ITS_INTERNALS
        // Inspect tables:
        printf( "  %-3s  %-3s  %-9s\n", "i", "alias", "prob" );
        for ( int i=0; i<n; ++i )
            printf( "  %3i  %3i  %9.7f\n",
                    i, ws->alias[i], ws->prob[i] );
    #endif

        // Draw M random samples;
        // accumulate statistic in histogram 'cumul':
        printf( "Drawing %i samples ...\n", M );
        double cumul[n];
        for ( i=0; i<n; ++i )
            cumul[i] = 0;
        for ( m=0; m<M; ++m ) {
            i = ransampl_draw( ws );
            cumul[i] += 1;
        }

        // Print given probability and obtained frequency:
        printf( "Result (input->output):\n");
        double sum = 0;
        for ( int i=0; i<n; ++i )
            sum += p[i];
        printf( "  %-18s  %-9s  %-9s  %-9s\n",
                "state", "N (Mio.)", "rel", "sim" );
        for ( int i=0; i<n; ++i )
            printf( "  %-18s  %9.7f  %9.7f  %9.7f\n",
                    names[i], p[i], p[i]/sum, ((double)cumul[i])/M );

        // Free workspace and terminate:
        ransampl_free( ws );
        return 0;
    }

COPYING

Software: FreeBSD License

Documentation: Creative Commons Attribution Share Alike

NAME

SYNOPSIS

DESCRIPTION

EXAMPLE

COPYING

SEE ALSO