Distributed probability random number generator

Question

I want to generate a number based on a distributed probability. For example, just say there are the following occurences of each numbers:

Number| Count           
1    |  150                
2    |  40          
3    |  15          
4    |  3  

with a total of (150+40+15+3) = 208     
then the probability of a 1 is 150/208= 0.72    
and the probability of a 2 is 40/208 = 0.192    

How do I make a random number generator that returns be numbers based on this probability distribution?

I'm happy for this to be based on a static, hardcoded set for now but I eventually want it to derive the probability distribution from a database query.

I've seen similar examples like this one but they are not very generic. Any suggestions?

Answer 1

The general approach is to feed uniformly distributed random numbers from 0..1 interval into the inverse of the cumulative distribution function of your desired distribution.

Thus in your case, just draw a random number x from 0..1 (for example with Random.NextDouble()) and based on its value return

• 1 if 0 <= x < 150/208,
• 2 if 150/208 <= x < 190/208,
• 3 if 190/208 <= x < 205/208 and
• 4 otherwise.

Answer 2

Do this only once:

• Write a function that calculates a cdf array given a pdf array. In your example pdf array is [150,40,15,3], cdf array will be [150,190,205,208].

Do this every time:

• Get a random number in [0,1) , multiply with 208, truncate up (or down: I leave it to you to think about the corner cases) You'll have an integer in 1..208. Name it r.
• Perform a binary search on cdf array for r. Return the index of the cell that contains r.

The running time will be proportional to log of the size of the given pdf array. Which is good. However, if your array size will always be so small (4 in your example) then performing a linear search is easier and also will perform better.

Answer 3

There are many ways to generate a random integer with a custom distribution (also known as a discrete distribution). The choice depends on many things, including the number of integers to choose from, the shape of the distribution, and whether the distribution will change over time.

One of the simplest ways to choose an integer with a custom weight function f(x) is the rejection sampling method. The following assumes that the highest possible value of f is max. The time complexity for rejection sampling is constant on average, but depends greatly on the shape of the distribution and has a worst case of running forever. To choose an integer in [1, k] using rejection sampling:

1. Choose a uniform random integer i in [1, k].
2. With probability f(i)/max, return i. Otherwise, go to step 1.

Other algorithms have an average sampling time that doesn't depend so greatly on the distribution (usually either constant or logarithmic), but often require you to precalculate the weights in a setup step and store them in a data structure. Some of them are also economical in terms of the number of random bits they use on average. These algorithms include the alias method, the Fast Loaded Dice Roller, the Knuth–Yao algorithm, the MVN data structure, and more. See my section "Weighted Choice With Replacement" for a survey.

---

The following C# code implements Michael Vose's version of the alias method, as described in this article; see also this question. I have written this code for your convenience and provide it here.

public class LoadedDie {
    // Initializes a new loaded die.  Probs
    // is an array of numbers indicating the relative
    // probability of each choice relative to all the
    // others.  For example, if probs is [3,4,2], then
    // the chances are 3/9, 4/9, and 2/9, since the probabilities
    // add up to 9.
    public LoadedDie(int probs){
        this.prob=new List<long>();
        this.alias=new List<int>();
        this.total=0;
        this.n=probs;
        this.even=true;
    }
    
    Random random=new Random();
    
    List<long> prob;
    List<int> alias;
    long total;
    int n;
    bool even;

    public LoadedDie(IEnumerable<int> probs){
        // Raise an error if nil
        if(probs==null)throw new ArgumentNullException("probs");
        this.prob=new List<long>();
        this.alias=new List<int>();
        this.total=0;
        this.even=false;
        var small=new List<int>();
        var large=new List<int>();
        var tmpprobs=new List<long>();
        foreach(var p in probs){
            tmpprobs.Add(p);
        }
        this.n=tmpprobs.Count;
        // Get the max and min choice and calculate total
        long mx=-1, mn=-1;
        foreach(var p in tmpprobs){
            if(p<0)throw new ArgumentException("probs contains a negative probability.");
            mx=(mx<0 || p>mx) ? P : mx;
            mn=(mn<0 || p<mn) ? P : mn;
            this.total+=p;
        }
        // We use a shortcut if all probabilities are equal
        if(mx==mn){
            this.even=true;
            return;
        }
        // Clone the probabilities and scale them by
        // the number of probabilities
        for(var i=0;i<tmpprobs.Count;i++){
            tmpprobs[i]*=this.n;
            this.alias.Add(0);
            this.prob.Add(0);
        }
        // Use Michael Vose's alias method
        for(var i=0;i<tmpprobs.Count;i++){
            if(tmpprobs[i]<this.total)
                small.Add(i); // Smaller than probability sum
            else
                large.Add(i); // Probability sum or greater
        }
        // Calculate probabilities and aliases
        while(small.Count>0 && large.Count>0){
            var l=small[small.Count-1];small.RemoveAt(small.Count-1);
            var g=large[large.Count-1];large.RemoveAt(large.Count-1);
            this.prob[l]=tmpprobs[l];
            this.alias[l]=g;
            var newprob=(tmpprobs[g]+tmpprobs[l])-this.total;
            tmpprobs[g]=newprob;
            if(newprob<this.total)
                small.Add(g);
            else
                large.Add(g);
        }
        foreach(var g in large)
            this.prob[g]=this.total;
        foreach(var l in small)
            this.prob[l]=this.total;
    }
    
    // Returns the number of choices.
    public int Count {
        get {
            return this.n;
        }
    }
    // Chooses a choice at random, ranging from 0 to the number of choices
    // minus 1.
    public int NextValue(){
        var i=random.Next(this.n);
        return (this.even || random.Next((int)this.total)<this.prob[i]) ? I : this.alias[i];
    }
}

Example:

 var loadedDie=new LoadedDie(new int[]{150,40,15,3}); // list of probabilities for each number:
                                                      // 0 is 150, 1 is 40, and so on
 int number=loadedDie.nextValue(); // return a number from 0-3 according to given probabilities;
                                   // the number can be an index to another array, if needed

I place this code in the public domain.

Answer 4

I know this is an old post, but I also searched for such a generator and was not satisfied with the solutions I found. So I wrote my own and want to share it to the world.

Just call "Add(...)" some times before you call "NextItem(...)"

/// <summary> A class that will return one of the given items with a specified possibility. </summary>
/// <typeparam name="T"> The type to return. </typeparam>
/// <example> If the generator has only one item, it will always return that item. 
/// If there are two items with possibilities of 0.4 and 0.6 (you could also use 4 and 6 or 2 and 3) 
/// it will return the first item 4 times out of ten, the second item 6 times out of ten. </example>
public class RandomNumberGenerator<T>
{
    private List<Tuple<double, T>> _items = new List<Tuple<double, T>>();
    private Random _random = new Random();

    /// <summary>
    /// All items possibilities sum.
    /// </summary>
    private double _totalPossibility = 0;

    /// <summary>
    /// Adds a new item to return.
    /// </summary>
    /// <param name="possibility"> The possibility to return this item. Is relative to the other possibilites passed in. </param>
    /// <param name="item"> The item to return. </param>
    public void Add(double possibility, T item)
    {
        _items.Add(new Tuple<double, T>(possibility, item));
        _totalPossibility += possibility;
    }

    /// <summary>
    /// Returns a random item from the list with the specified relative possibility.
    /// </summary>
    /// <exception cref="InvalidOperationException"> If there are no items to return from. </exception>
    public T NextItem()
    {
        var rand = _random.NextDouble() * _totalPossibility;
        double value = 0;
        foreach (var item in _items)
        {
            value += item.Item1;
            if (rand <= value)
                return item.Item2;
        }
        return _items.Last().Item2; // Should never happen
    }
}

Answer 5

Thanks for all your solutions guys! Much appreciated!

@Menjaraz I tried implementing your solution as it looks very resource friendly, however had some difficulty with the syntax.

So for now, I just transformed my summary into a flat list of values using LINQ SelectMany() and Enumerable.Repeat().

public class InventoryItemQuantityRandomGenerator
{
    private readonly Random _random;
    private readonly IQueryable<int> _quantities;

    public InventoryItemQuantityRandomGenerator(IRepository database, int max)
    {
        _quantities = database.AsQueryable<ReceiptItem>()
            .Where(x => x.Quantity <= max)
            .GroupBy(x => x.Quantity)
            .Select(x => new
                             {
                                 Quantity = x.Key,
                                 Count = x.Count()
                             })
            .SelectMany(x => Enumerable.Repeat(x.Quantity, x.Count));

        _random = new Random();
    }

    public int Next()
    {
        return _quantities.ElementAt(_random.Next(0, _quantities.Count() - 1));
    }
}

Answer 6

Use my method. It is simple and easy-to-understand. I don't count portion in range 0...1, i just use "Probabilityp Pool" (sounds cool, yeah?)

At circle diagram you can see weight of every element in pool

Here you can see an implementing of accumulative probability for roulette

`

// Some c`lass or struct for represent items you want to roulette
public class Item
{
    public string name; // not only string, any type of data
    public int chance;  // chance of getting this Item
}

public class ProportionalWheelSelection
{
    public static Random rnd = new Random();

    // Static method for using from anywhere. You can make its overload for accepting not only List, but arrays also: 
    // public static Item SelectItem (Item[] items)...
    public static Item SelectItem(List<Item> items)
    {
        // Calculate the summa of all portions.
        int poolSize = 0;
        for (int i = 0; i < items.Count; i++)
        {
            poolSize += items[i].chance;
        }

        // Get a random integer from 0 to PoolSize.
        int randomNumber = rnd.Next(0, poolSize) + 1;

        // Detect the item, which corresponds to current random number.
        int accumulatedProbability = 0;
        for (int i = 0; i < items.Count; i++)
        {
            accumulatedProbability += items[i].chance;
            if (randomNumber <= accumulatedProbability)
                return items[i];
        }
        return null;    // this code will never come while you use this programm right :)
    }
}

// Example of using somewhere in your program:
        static void Main(string[] args)
        {
            List<Item> items = new List<Item>();
            items.Add(new Item() { name = "Anna", chance = 100});
            items.Add(new Item() { name = "Alex", chance = 125});
            items.Add(new Item() { name = "Dog", chance = 50});
            items.Add(new Item() { name = "Cat", chance = 35});

            Item newItem = ProportionalWheelSelection.SelectItem(items);
        }

Answer 7

Here's an implementation using the Inverse distribution function:

using System;
using System.Linq;

    // ...
    private static readonly Random RandomGenerator = new Random();

    private int GetDistributedRandomNumber()
    {
        double totalCount = 208;

        var number1Prob = 150 / totalCount;
        var number2Prob = (150 + 40) / totalCount;
        var number3Prob = (150 + 40 + 15) / totalCount;

        var randomNumber = RandomGenerator.NextDouble();

        int selectedNumber;

        if (randomNumber < number1Prob)
        {
            selectedNumber = 1;
        }
        else if (randomNumber >= number1Prob && randomNumber < number2Prob)
        {
            selectedNumber = 2;
        }
        else if (randomNumber >= number2Prob && randomNumber < number3Prob)
        {
            selectedNumber = 3;
        }
        else
        {
            selectedNumber = 4;
        }

        return selectedNumber;
    }

An example to verify the random distribution:

        int totalNumber1Count = 0;
        int totalNumber2Count = 0;
        int totalNumber3Count = 0;
        int totalNumber4Count = 0;

        int testTotalCount = 100;

        foreach (var unused in Enumerable.Range(1, testTotalCount))
        {
            int selectedNumber = GetDistributedRandomNumber();

            Console.WriteLine($"selected number is {selectedNumber}");

            if (selectedNumber == 1)
            {
                totalNumber1Count += 1;
            }

            if (selectedNumber == 2)
            {
                totalNumber2Count += 1;
            }

            if (selectedNumber == 3)
            {
                totalNumber3Count += 1;
            }

            if (selectedNumber == 4)
            {
                totalNumber4Count += 1;
            }
        }

        Console.WriteLine("");
        Console.WriteLine($"number 1 -> total selected count is {totalNumber1Count} ({100 * (totalNumber1Count / (double) testTotalCount):0.0} %) ");
        Console.WriteLine($"number 2 -> total selected count is {totalNumber2Count} ({100 * (totalNumber2Count / (double) testTotalCount):0.0} %) ");
        Console.WriteLine($"number 3 -> total selected count is {totalNumber3Count} ({100 * (totalNumber3Count / (double) testTotalCount):0.0} %) ");
        Console.WriteLine($"number 4 -> total selected count is {totalNumber4Count} ({100 * (totalNumber4Count / (double) testTotalCount):0.0} %) ");

Example output:

selected number is 1
selected number is 1
selected number is 1
selected number is 1
selected number is 2
selected number is 1
...
selected number is 2
selected number is 3
selected number is 1
selected number is 1
selected number is 1
selected number is 1
selected number is 1

number 1 -> total selected count is 71 (71.0 %) 
number 2 -> total selected count is 20 (20.0 %) 
number 3 -> total selected count is 8 (8.0 %) 
number 4 -> total selected count is 1 (1.0 %)