443,730 Members | 1,559 Online Need help? Post your question and get tips & solutions from a community of 443,730 IT Pros & Developers. It's quick & easy.

# function to generate 1 or 0 according to a given probability

 P: n/a hi everyone, im trying to create a function that generates a 1 or a 0, with the probability of a 1 being generated equal to X (which is passed to the function as a parameter). any ideas? thanks in anticipation, Mayank Nov 15 '05 #1
11 Replies

 P: n/a Mayank Kaushik wrote: hi everyone, im trying to create a function that generates a 1 or a 0, with the probability of a 1 being generated equal to X (which is passed to the function as a parameter). any ideas? thanks in anticipation, /* BEGIN new.c */ #include #include int X_prob(double x); int main(void) { unsigned count; printf("\n0.95 "); for (count = 0; count != 70; ++count) { printf("%d", X_prob(0.95)); } printf("\n\n0.50 "); for (count = 0; count != 70; ++count) { printf("%d", X_prob(0.5)); } printf("\n\n0.05 "); for (count = 0; count != 70; ++count) { printf("%d", X_prob(0.05)); } putchar('\n'); return 0; } int X_prob(double x) { return RAND_MAX * x >= rand(); } /* END new.c */ -- pete Nov 15 '05 #2

 P: n/a "Mayank Kaushik" wrote in news:11*********************@f14g2000cwb.googlegro ups.com: im trying to create a function that generates a 1 or a 0, with the probability of a 1 being generated equal to X (which is passed to the function as a parameter). Assuming rand generates uniformly distributed values, this is trivial. I suspect that the question might be off-topic in clc: comp.programming might be a better place to ask this question. #include #include #include int success(double prob) { return prob >= ((double) rand())/RAND_MAX; } int main(void) { int i; srand(time(NULL)); for(i = 0; i < 10; ++i) { printf("%s\n", success(0.38) ? "success" : "failure"); } return 0; } Sinan -- A. Sinan Unur <1u**@llenroc.ude.invalid> (reverse each component and remove .invalid for email address) Nov 15 '05 #3

 P: n/a A. Sinan Unur wrote: "Mayank Kaushik" wrote in news:11*********************@f14g2000cwb.googlegro ups.com: im trying to create a function that generates a 1 or a 0, with the probability of a 1 being generated equal to X (which is passed to the function as a parameter). Assuming rand generates uniformly distributed values, this is trivial. I suspect that the question might be off-topic in clc: comp.programming might be a better place to ask this question. #include #include #include int success(double prob) { return prob >= ((double) rand())/RAND_MAX; } What if prob == 0.0 and the call to rand() returns 0? Nov 15 '05 #4

 P: n/a pete wrote: int X_prob(double x) { return RAND_MAX * x >= rand(); } Dingo's post made me think that either 0.0 or 1.0 would probably be handled better as a special case. return x >= 1.0 ? 1 : RAND_MAX * x > rand(); -- pete Nov 15 '05 #5

 P: n/a thanks guys, i didnt know about RAND_MAX, i stand enlightened. Mayank UT-Austin Nov 15 '05 #6

 P: n/a pete writes: pete wrote: int X_prob(double x) { return RAND_MAX * x >= rand(); } Dingo's post made me think that either 0.0 or 1.0 would probably be handled better as a special case. return x >= 1.0 ? 1 : RAND_MAX * x > rand(); The multiplication isn't quite right; fencepost error. The probabilities will be a little off. Fixing up the fencepost error gets rid of the need for special casing the boundary conditions: int X_prob( double x ){ assert( 0 <= x && x <= 1 ); return (RAND_MAX + 1.0) * x - 0.5 > rand(); } Now any probability less than 1 / (2 * (RAND_MAX + 1.0)) will always yield zero, and any probability greater than 1 - (1 / (2 * (RAND_MAX + 1.0))) will always yield one, which is the best that can be done under the circumstances. (That is, "best" in sense of delivering the closest probability possible, given that the function doesn't save any state and calls rand() only once.) Nov 15 '05 #7

 P: n/a Tim Rentsch wrote: pete writes: pete wrote: int X_prob(double x) { return RAND_MAX * x >= rand(); } Dingo's post made me think that either 0.0 or 1.0 would probably be handled better as a special case. return x >= 1.0 ? 1 : RAND_MAX * x > rand(); The multiplication isn't quite right; fencepost error. The probabilities will be a little off. Fixing up the fencepost error gets rid of the need for special casing the boundary conditions: int X_prob( double x ){ assert( 0 <= x && x <= 1 ); return (RAND_MAX + 1.0) * x - 0.5 > rand(); } Now any probability less than 1 / (2 * (RAND_MAX + 1.0)) will always yield zero, and any probability greater than 1 - (1 / (2 * (RAND_MAX + 1.0))) will always yield one, which is the best that can be done under the circumstances. (That is, "best" in sense of delivering the closest probability possible, given that the function doesn't save any state and calls rand() only once.) Thank you. -- pete Nov 15 '05 #8

 P: n/a "Dingo" wrote in news:1131495313.121949.181570 @g49g2000cwa.googlegroups.com: A. Sinan Unur wrote: .... int success(double prob) { return prob >= ((double) rand())/RAND_MAX; } What if prob == 0.0 and the call to rand() returns 0? Good call. I see others have already posted improved versions of such a routine. Sinan Nov 15 '05 #9

 P: n/a pete wrote: pete wrote: int X_prob(double x) { return RAND_MAX * x >= rand(); } Dingo's post made me think that either 0.0 or 1.0 would probably be handled better as a special case. return x >= 1.0 ? 1 : RAND_MAX * x > rand(); How is every falling into this trap? What if 0 < x < 0.5 / RAND_MAX ? You can read more about a similar problem here: http://www.azillionmonkeys.com/qed/random.html There there is a function defined there called drand() which will solve the problem far more precisely and for a much larger range, and will usually be sufficient in practice: return drand() < x; But in general, to get an exact probability, you may have to call a discrete PRNG many times depending on your desired precision (unless your probability is a multiple of a power of a negative power of 2, in which case you can guarantee the exact probability.) -- Paul Hsieh http://www.pobox.com/~qed/ http://bstring.sf.net/ Nov 15 '05 #10

 P: n/a In article <11**********************@g47g2000cwa.googlegroups .com>, we******@gmail.com writes: But in general, to get an exact probability, you may have to call a discrete PRNG many times depending on your desired precision (unless your probability is a multiple of a power of a negative power of 2, in which case you can guarantee the exact probability.) A few months back someone posted to sci.crypt an algorithm for deriving a random value with exact probability while consuming only the minimal number of (pseudo)random bits. I think it assumed bits were unbiased, but that can be corrected deterministically if the bias is known. (If the bias is unknown, it can be corrected using eg von Neumann's method, but that's non-deterministic and may require discarding an arbitrary number of bits.) I didn't look at the algorithm closely and I don't recall the details, but it may be of interest. -- Michael Wojcik mi************@microfocus.com Q: What is the derivation and meaning of the name Erwin? A: It is English from the Anglo-Saxon and means Tariff Act of 1909. -- Columbus (Ohio) Citizen Nov 15 '05 #11

 P: n/a Michael Wojcik wrote: we******@gmail.com writes: But in general, to get an exact probability, you may have to call a discrete PRNG many times depending on your desired precision (unless your probability is a multiple of a power of a negative power of 2, in which case you can guarantee the exact probability.) A few months back someone posted to sci.crypt an algorithm for deriving a random value with exact probability while consuming only the minimal number of (pseudo)random bits. Yeah, I don't think its so hard that I need to appeal to the experts in sci.crypt. I have augmented my page to include a trivial way to transform rand(), to a binary output rand with a settable bias (see randbias()). In general translating one bias to another can be seen geometrically -- just draw a number line and plot out a length equal to the probability you want, mark off tiled slots corresponding to the probability of each output for your source rand(), get a sample and if the sample corresponds to a slot that doesn't straddle the position where your desired probability is return according the best weighted estimate. If it does straddle, scale that slot to [0,1] and recurse. [...] I think it assumed bits were unbiased, but that can be corrected deterministically if the bias is known. (If the bias is unknown, it can be corrected using eg von Neumann's method, but that's non-deterministic and may require discarding an arbitrary number of bits.) Perhaps the method I describe is the same thing? Whatever, its not hard to figure out. -- Paul Hsieh http://www.pobox.com/~qed/ http://bstring.sf.net/ Nov 15 '05 #12

### This discussion thread is closed

Replies have been disabled for this discussion. 