function to generate 1 or 0 according to a given probability

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

  1. 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
    #1
    1. Advertising

  2. Mayank Kaushik

    pete Guest

    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
    #2
    1. Advertising

  3. "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
    #3
  4. Mayank Kaushik

    Dingo Guest

    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
    #4
  5. Mayank Kaushik

    pete Guest

    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
    #5
  6. thanks guys, i didnt know about RAND_MAX, i stand enlightened.

    Mayank
    UT-Austin
    Mayank Kaushik, Nov 9, 2005
    #6
  7. Mayank Kaushik

    Tim Rentsch Guest

    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
    #7
  8. Mayank Kaushik

    pete Guest

    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
    #8
  9. "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
    A. Sinan Unur, Nov 9, 2005
    #9
  10. Mayank Kaushik

    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 ?
    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 9, 2005
    #10
  11. 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
    #11
  12. Mayank Kaushik

    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
    #12
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. Replies:
    1
    Views:
    331
    el goog
    Mar 20, 2005
  2. Lord0
    Replies:
    1
    Views:
    562
    Thomas Weidenfeller
    Apr 19, 2006
  3. Farraige
    Replies:
    4
    Views:
    284
    Farraige
    Nov 8, 2006
  4. chiara
    Replies:
    6
    Views:
    465
    Barry Schwarz
    Oct 6, 2005
  5. 2Barter.net
    Replies:
    0
    Views:
    363
    2Barter.net
    Dec 13, 2006
Loading...

Share This Page