Weighted Random Number Generator algorithm

// WeightedPRNG.cpp : Definiert den Einstiegspunkt für die Konsolenanwendung.
//

#include "stdafx.h"
#include <vector>
#include <iostream>
#include <iomanip>

#define USE_CUSTOM_RNG

#ifdef USE_CUSTOM_RNG
#include "RandomNumberGenerator.h"
#include <windows.h>
#endif

typedef long double Number;

Number random_unit()
{
#ifdef USE_CUSTOM_RNG
    // this actually wraps around a boost Mersenne Twister PRNG implementation
    // which will give better distributed numbers
    // see http://www.boost.org/doc/libs/1_46_1/doc/html/boost_random.html
    static RandomNumberGenerator rng( (size_t)GetTickCount64() );

    return rng.random_double_unit();
#else
    // pretty bad resolution (16bit)
    // bad distribution pattern
    return rand() / (Number)RAND_MAX;
#endif
}

struct tWeightedItem
{
    tWeightedItem( const std::string& name, Number weight )
        : choses(0), name(name), weight(weight){};

    std::string name;
    Number weight;
    Number rel_weight;
    int choses;
};

typedef std::vector<tWeightedItem> WeightedItemList;

/// Actual weighted random choosing algorithm implementation
/// Works like cutting all items from paper with lengths
/// equalling their weight, then aligning them and throwing
/// a random dart at it:
///                                   v Dart
///  Strip  |-------|----|---------|---x-|------|--------|
///  Item#      #0     #1     #2      #3    #4       #5
/// Weight      8      5      10      6     7        9
WeightedItemList::iterator choose( WeightedItemList& wil )
{
    Number weight_sum = 0.0f;

    // sum up weights
    for( size_t n=0; n<wil.size(); n++ )
        weight_sum += wil[n].weight;

    // probe a random point within our total strip length,
    // random_double_unit will return a number in [0,1] (closed interval)
    Number random_choice = random_unit() * weight_sum;

    Number sum = 0.0f;

    // The loop will find the item where the probe number hits the
    // range of the weight
    for( WeightedItemList::iterator it=wil.begin(); it != wil.end(); ++it )
    {
        const tWeightedItem& wc = *it;

        // I'm not sure about the equality operators here, because some
        // numbers fall into two categories. But as the PRNG numbers are
        // inside [0,1] I have to include the leftmost and rightmost borders
        // of the number range which only works with this operator combination.
        // If your PRNG spits out numbers in (0,1] or [0,1) you might want
        // to adjust the condition.
        if( random_choice >= sum && random_choice <= (sum + wc.weight) )
            return it;

        sum += wc.weight;
    }

    // should never happen, if it does, the
    // random probe will
    return wil.end();
}

// Algorithm test function, checks if the weights assigned equal the
// actual stochastic occurences
int main(int argc, char* argv[])
{
    // number of stochastic samples to take
    // the higher, the more should the actual distribution
    // converge to the intended weight
    const size_t SAMPLES = 10000000;

    WeightedItemList items;

    // Some test probands
    items.push_back( tWeightedItem( "Hans", 1));
    items.push_back( tWeightedItem( "Frank", 7));
    items.push_back( tWeightedItem( "Joe", 2));
    items.push_back( tWeightedItem( "Lisa", 5));
    items.push_back( tWeightedItem( "Othello", 10));

    Number total = 0.0f;

    for( size_t n=0; n<items.size(); n++ )
    {
        total += items[n].weight;
    }

    for( size_t n=0; n<items.size(); n++ )
    {
        items[n].rel_weight = items[n].weight / total;
    }

    int fails = 0;

    for( size_t n=0; n<SAMPLES; n++ )
    {
        WeightedItemList::iterator chosen = choose( items );

        if( chosen != items.end() )
            chosen->choses++;
        else
            fails++;
    }

    std::cout << std::setw(15) << "Item" <<
        " | " << std::setprecision(5) << std::setw(10) << std::fixed << "Expected" <<
        " | " << std::setprecision(5) << std::setw(10) << std::fixed << "Real" <<
        " | " << std::setprecision(5) << std::setw(10) << std::fixed << "Difference" << std::endl;

    Number sum_expected = 0;
    Number sum_real = 0;
    Number sum_differences = 0;

    for( size_t n=0; n<items.size(); n++ )
    {
        tWeightedItem& wi = items[n];

        Number possibility = (wi.choses / (Number)SAMPLES);
        Number difference = fabs(possibility - wi.rel_weight);

        std::cout << std::setw(15) << wi.name.c_str() <<
            " | " << std::setprecision(5) << std::setw(10) << std::fixed << wi.rel_weight <<
            " | " << std::setprecision(5) << std::setw(10) << std::fixed << possibility <<
            " | " << std::setprecision(5) << std::setw(10) << std::fixed << difference;

        if( fabs(wi.rel_weight - (wi.choses / (Number)SAMPLES)) > 0.001f)
        {
            std::cout << " [!]";
        }

        sum_expected += wi.rel_weight;
        sum_real += possibility;
        sum_differences += difference;

        std::cout << std::endl;
    }

    std::cout << std::setw(15) << "Totals" <<
        " | " << std::setprecision(5) << std::setw(10) << std::fixed << sum_expected <<
        " | " << std::setprecision(5) << std::setw(10) << std::fixed << sum_real <<
        " | " << std::setprecision(5) << std::setw(10) << std::fixed << sum_differences << std::endl;

    std::cout << "Samples: " << SAMPLES << std::endl << "Fails: " << fails << std::endl;

    return 0;
}