function to generate 1 or 0 according to a given probability

[Mayank Kaushik]

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]

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]

"Mayank Kaushik" <pr***************@yahoo.com> 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 <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 <1u**@llenroc.ude.invalid>
(reverse each component and remove .invalid for email address)

[Dingo]

A. Sinan Unur wrote:

"Mayank Kaushik" <pr***************@yahoo.com> 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 <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]

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]

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

Mayank
UT-Austin

[Tim Rentsch]

pete <pf*****@mindspring.com> 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]

Tim Rentsch wrote:

pete <pf*****@mindspring.com> 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]

"Dingo" <di*************@aol.com> 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

[websnarf]

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/

[Michael Wojcik]

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

[websnarf]

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/