DQNOK wrote:
Do you know the algorithm for using the rand() function (which generates a uniform PDF) to generate values within a specified PDF?
I answered my own question. For a PDF called 'y(x)' (where x is the random varaible), and a computer library function rand() that returns a random value in the range [0, 1) from a uniform PDF, we seek to map rand's return value (call it rx) to the correct x. Effectively, we are looking for x such that the area under the curve y(x) LEFT OF x is equal to the area under the uniform PDF left of rx. But since rand produces values within the unit interval [0,1), rx IS the area under the curve left of rx. So, we seek x such that: rx = integral from -infinity to x of y(x) Exactly HOW we evaluate the integral (i.e. solve for x) is an implementation detail that depends on how y(x) is defined. If the integral is already given via a formula, we just use a standard non-linear scalar solver, like a Newton-Raphson. It not, then perhaps I would just use a stepping approach, something like:ndx = getBestStartingIndex(rx); x = X_STARTS[ndx]; // X_STARTS[ndx] <= rx < X_STARTS[ndx+1] area = AREAS[ndx]; // precomputed areas based on x. old_y = y(x); while( area < rx ) { x += dx; new_y = y(x); area += dx * (old_y + new_y)/2 ; //trapezoid rule old_y = new_y; } return x;
Which can be polished for better precision, time-performance etc.
David --------- Empirical studies indicate that 20% of the people drink 80% of the beer. With C++ developers, the rule is that 80% of the developers understand at most 20% of the language. It is not the same 20% for different people, so don't count on them to understand each other's code. http://yosefk.com/c++fqa/picture.html#fqa-6.6 ---------
modified on Thursday, April 17, 2008 9:28 AM