0

I am trying to generate some uniform real numbers for a Monte Carlo integration but the routine I build was returning some really strange values. Upon closer inspection I notices that Boost was returning some crazy looking random numbers e.g.:

     temp = -0.185276
     temp = -0.864523
     temp = -0.0942081
     temp = -0.164991
     temp = -0.873013
     temp = -0.0311322
     temp = -0.0866241
     temp = -0.778966
     temp = -0.367641
     temp = -0.691833
     temp = 5.66499e-310
     temp = 9.42007e-311
     temp = 6.29821e-310
     temp = 5.80603e-310
     temp = 8.82973e-311
     temp = 6.73679e-310
     temp = 6.35094e-310
     temp = 1.53691e-310
     temp = 4.39696e-310
     temp = 2.14277e-310

Whilst these numbers are technically still reals generated between the bounds -1 and 1 I would prefer it if they weren't quite so small!

My implementation of the call to boost is in a function which is called multiple times (for different bounding values) as follows:

// Define Boost typedefs
typedef boost::mt19937 Engine;
typedef boost::uniform_real<double> Distribution;
typedef boost::variate_generator <Engine, Distribution> Generator;

int main (void) {
    ...
    Integral = MCRecursion(...);
    ...
    return 0;
}

double MCRecursion (int Count, double Lower, double Upper, double (*Integrand)(double)) {

// Define Boost objects
Engine       Eng;
Distribution Dist (Lower, Upper);
Generator    RandomGen (Eng, Dist);

Eng.seed(time(0));

// Variables for Monte Carlo sample sums
double Sum = 0.0;
double temp;

for (int i = 0; i < Count; i++) {
    temp = RandomGen();
    std::cout << "         temp = " << temp << std::endl;
    Sum += Integrand(temp);
}

return (Upper - Lower) * Sum / Count;
}

I assume the problem is something with my implementation but I can't find any errors. Any and all help appreciated! Cheers, Jack

EDIT

Code for calling MCRecursion:

The Code I am writting runs a Monte Carlo on the entire domain I am interested in [Lower, Upper] and then looks again at the left half of the whole domain and the right half of the domain. e.g. if we were integrating f(x) between -a and a I calculate the full integral using:

double FullDomain = MCRecursion (1e5, LowerBound, UpperBound, f);
double Centre = (Upper + Lower) / 2.0;
double LeftHalf = MCRecursion (1e5, LowerBound, Centre, f);
double RightHalf = MCRecursion (1e5, Centre, UpperBound, f);

and I then look at the uncertainty by calculating: double difference = fabs(FullDomain - LeftHalf - Righthalf); to see if more samples is 'worth it' in some sense Jack

JMzance
  • 1,704
  • 4
  • 30
  • 49
  • Well that seems fine to me. If you want your values to have a guaranteed minimal absolute value, you have to create a distribution on [eps, 1] and another one for choosing the sign. Are the values causing any problems? – stefan Dec 30 '13 at 21:12
  • 1
    check out http://stackoverflow.com/questions/15747194/why-is-boosts-random-number-generation-on-a-normal-distribution-always-giving – tinkertime Dec 30 '13 at 21:14
  • Yes @stefan because a load of the numbers are ~0.0 whatever I function I integrate is hugely oversampled near zero so let's say I was integrating f(x)=x then my Monte Carlo sum would be significantly less than one would expect. – JMzance Dec 30 '13 at 21:16
  • @JackMedley You chose a `uniform_real` distribution, so if boost implements it correctly, it shouldn't be mostly around 0.0. Well as I've said: if you want to establish a minimum on the absolute value you need to either ignore values which are too small or mirror a range [eps, 1]. – stefan Dec 30 '13 at 21:18
  • 1
    Can you please show the code calling `MCRescursion` and make sure the output is actually output of your posted code? The problem is not reproducable for me, I made a test case and got evenly distributed random numbers. –  Dec 30 '13 at 21:19
  • 1
    I guess another way of saying what @Nabla's saying is "please construct a [minimal test-case](http://sscce.org)" ;) – Oliver Charlesworth Dec 30 '13 at 21:20
  • @JackMedley It is still not a minimal test-case as I cannot copy-paste it and compile it. You have two sets of bounds: `Upper`/`Lower` and `UpperBound`/`LowerBound`, is this intended? What values did they have in your test run? Is the output from one or two calls to `MCRecursion`? –  Dec 30 '13 at 21:29
  • Yes @zahir I have a little 'if' statement that throws an error if something weird happens with my bounds – JMzance Dec 30 '13 at 21:29
  • Hey @Nabla here is the full thing....should compile fine. – JMzance Dec 30 '13 at 21:31
  • http://pastebin.com/WycKapx4 – JMzance Dec 30 '13 at 21:32
  • They are called difference things in different functions but they DO change, they focus in on smaller and smaller sections of the full domain in an attempt to look for regions of high variance – JMzance Dec 30 '13 at 21:33
  • @JackMedley Please post code in the question. Also have you read the link Oli Charlesworth posted. This is not a minimal test-case. –  Dec 30 '13 at 21:33
  • That said, your code throws several warnings. You use `Right` uninitialized. Also I made some debug output. At one point your code calls `MCSample` with `Lower=0` and `Upper=6.91473e-310` –  Dec 30 '13 at 21:36
  • Hey @Nabla check out my pastebin link above for a minimal working example, that's the full thing (still fairly small file). – JMzance Dec 30 '13 at 21:37
  • Hey @Nabla thanks, I tried printing out these values but didnt get that value, could you let me know what line you're on? – JMzance Dec 30 '13 at 21:41

1 Answers1

1

Based on the pastebin the questioner posted in the comments:

This is not a problem with the random library but rather a simple programming error. Compiling the code throws the warning:

../src/Test.cpp: In function ‘double AdaptiveMCRecursion(double, double, double, double, int, double, double (*)(double))’:
../src/Test.cpp:100:72: warning: ‘Right’ is used uninitialized in this function [-Wuninitialized]
         double Right = MCSample (Count, Central, Right,   Integrand);

So all the behaviour from that line on is basically undefined. Especially it results in calling the function MCSample with an undetermined Upper parameter. So your result is not unexpected. You are actually lucky the program runs at all.