function to generate 1 or 0 according to a given probability

Discussion in 'C Programming' started by Mayank Kaushik, Nov 8, 2005.

1. Mayank KaushikGuest

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

Mayank Kaushik, Nov 8, 2005

2. peteGuest

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

pete, Nov 8, 2005

3. A. Sinan UnurGuest

"Mayank Kaushik" <> wrote in
news::

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

return 0;
}

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

A. Sinan Unur, Nov 8, 2005
4. DingoGuest

A. Sinan Unur wrote:
> "Mayank Kaushik" <> wrote in
> news::
>
> > 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?

Dingo, Nov 9, 2005
5. peteGuest

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

pete, Nov 9, 2005
6. Mayank KaushikGuest

thanks guys, i didnt know about RAND_MAX, i stand enlightened.

Mayank
UT-Austin

Mayank Kaushik, Nov 9, 2005
7. Tim RentschGuest

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

Tim Rentsch, Nov 9, 2005
8. peteGuest

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

pete, Nov 9, 2005
9. A. Sinan UnurGuest

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

A. Sinan Unur, Nov 9, 2005
10. Guest

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/

, Nov 9, 2005
11. Michael WojcikGuest

In article <>, 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

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

Michael Wojcik, Nov 10, 2005
12. Guest

Michael Wojcik wrote:
> 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 13, 2005