How to pick a number based on probability?

Question

I want to select a random number from 0,1,2,3...n, however I want to make it that the chance of selecting k|0<k<n will be lower by multiplication of x from selecting k - 1 so x = (k - 1) / k. As bigger the number as smaller the chances to pick it up.

As an answer I want to see the implementation of the next method:

int pickANumber(n,x)

This is for a game that I am developing, I saw those questions as related but not exactly that same:

• How to pick an item by its probability
• C Function for picking from a list where each element has a distinct probabili

Answer 1

p1 + p2 + ... + pn = 1
p1 = p2 * x
p2 = p3 * x
...
p_n-1 = pn * x

Solving this gives you:

p1 + p2 + ... + pn = 1
(p2 * x) + (p3 * x) + ... + (pn * x) + pn = 1
((p3*x) * x) + ((p4*x) * x) + ... + ((p_n-1*x) * x) + pn = 1
....
pn* (x^(n-1) + x^(n-2) + ... +x^1 + x^0) = 1
pn*(1-x^n)/(1-x) = 1
pn = (1-x)/(1-x^n)

This gives you the probability you need to set to pn, and from it you can calculate the probabilities for all other p1,p2,...p_n-1

Now, you can use a "black box" RNG that chooses a number with a distribution, like those in the threads you mentioned.

A simple approach to do it is to set an auxillary array:

aux[i] = p1 + p2 + ... + pi

Now, draw a random number with uniform distribution between 0 to aux[n], and using binary search (aux array is sorted), get the first value, which matching value in aux is greater than the random uniform number you got

---

Original answer, for substraction (before question was editted):

For n items, you need to solve the equation:

p1 + p2 + ... + pn = 1
p1 = p2 + x
p2 = p3 + x
...
p_n-1 = pn + x

Solving this gives you:

p1 + p2 + ... + pn = 1
(p2 + x) + (p3 + x) + ... + (pn + x) + pn = 1
((p3+x) + x) + ((p4+x) + x) + ... + ((p_n-1+x) + x) + pn = 1
....
pn* ((n-1)x + (n-2)x + ... +x + 0) = 1
pn* x = n(n-1)/2
pn = n(n-1)/(2x)

This gives you the probability you need to set to pn, and from it you can calculate the probabilities for all other p1,p2,...p_n-1

Now, you can use a "black box" RNG that chooses a number with a distribution, like those in the threads you mentioned.

---

Be advised, this is not guaranteed you will have a solution such that 0<p_i<1 for all i, but you cannot guarantee one given from your requirements, and it is going to depend on values of n and x to fit.

Answer 2

Edit This answer was for the OPs original question, which was different in that each probability was supposed to be lower by a fixed amount than the previous one.

Well, let's see what the constraints say. You want to have P(k) = P(k - 1) - x. So we have:

P(0)

P(1) = P(0) - x

P(2) = P(0) - 2x ...

In addition, Sum_k P(k) = 1. Summing, we get:

1 = (n + 1)P(0) -x * n / 2 (n + 1),

This gives you an easy constraint between x and P(0). Solve for one in terms of the other.

Answer 3

For this I would use the Mersenne Twister algorithm for a uniform distribution which Boost provides, then have a mapping function to map the results of that random distribution to the actual number select.

Here's a quick example of a potential implementation, although I left out the quadtratic equation implementation since it is well known:

int f_of_xib(int x, int i, int b)
{
    return x * i * i / 2 + b * i;
}

int b_of_x(int i, int x)
{
    return (r - ( r ) / 2 );
}


int pickANumber(mt19937 gen, int n, int x)
{
    // First, determine the range r required where the probability equals i * x
    // since probability of each increasing integer is x higher of occuring.
    // Let f(i) = r and given f'(i) = x * i then r = ( x * i ^2 ) / 2 + b * i
    // where b = ( r - ( x * i ^ 2 ) / 2 ) / i . Since r = x when i = 1 from problem
    // definition, this reduces down to b = r - r / 2. therefore to find r_max simply
    // plugin x to find b, then plugin n for i, x, and b to get r_max since r_max occurs
    // when n == i.

    // Find b when 
    int b = b_of_x(x);
    int r_max = f_of_xib(x, n, b);

    boost::uniform_int<> range(0, r_max);
    boost::variate_generator<boost::mt19937&, boost::uniform_int<> > next(gen, range);

    // Now to map random number to desired number, just find the positive value for i
    // when r is the return random number which boils down to finding the non-zero root
    // when 0 = ( x * i ^ 2 ) / 2 + b * i - r
    int random_number = next();

    return quadtratic_equation_for_positive_value(1, b, r);
}



int main(int argc, char** argv)
{
    mt19937 gen;
    gen.seed(time(0));

    pickANumber(gen, 10, 1);

    system("pause");
}