Velocity Reviews > function to generate 1 or 0 according to a given probability

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

Mayank Kaushik
Guest
Posts: n/a

 11-08-2005
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

pete
Guest
Posts: n/a

 11-08-2005
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 <stdio.h>
#include <stdlib.h>

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

A. Sinan Unur
Guest
Posts: n/a

 11-08-2005
"Mayank Kaushik" <(E-Mail Removed)> wrote in
news:(E-Mail Removed) 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 <stdio.h>
#include <stdlib.h>
#include <time.h>

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.3 ? "success" : "failure");
}

return 0;
}

Sinan
--
A. Sinan Unur <(E-Mail Removed)>
(reverse each component and remove .invalid for email address)

Dingo
Guest
Posts: n/a

 11-09-2005

A. Sinan Unur wrote:
> "Mayank Kaushik" <(E-Mail Removed)> wrote in
> news:(E-Mail Removed) 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 <stdio.h>
> #include <stdlib.h>
> #include <time.h>
>
> int success(double prob) {
> return prob >= ((double) rand())/RAND_MAX;
> }
>

What if prob == 0.0 and the call to rand() returns 0?

pete
Guest
Posts: n/a

 11-09-2005
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

Mayank Kaushik
Guest
Posts: n/a

 11-09-2005
thanks guys, i didnt know about RAND_MAX, i stand enlightened.

Mayank
UT-Austin

Tim Rentsch
Guest
Posts: n/a

 11-09-2005
pete <(E-Mail Removed)> 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.)

pete
Guest
Posts: n/a

 11-09-2005
Tim Rentsch wrote:
>
> pete <(E-Mail Removed)> 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

A. Sinan Unur
Guest
Posts: n/a

 11-09-2005
"Dingo" <(E-Mail Removed)> wrote in news:1131495313.121949.181570

>
> 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

websnarf@gmail.com
Guest
Posts: n/a

 11-09-2005
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 ?

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/