0

I'm neither a statistician nor a mathematician so apologies if this question is asked and answered elsewhere or if I've misdescribed my problem. I am trying to figure out in a thread in stackoverflow how to compute all possible permutations and their probabilities when $b$ indistinguishable balls are randomly inserted into $B$ indistinguishable buckets where each bucket has a maximum size $z$.

I'm mostly looking for algorithms so I can code this up but first I need to make sure I understand probability correctly. To check the correctness of the code, I have made the following probability trees for $b=3$, $B=3$, $z=3$ and for $b=3$, $B=3$, $z=2$:

bins and buckets http://bents.us/Pictures/bins_balls.png

Where each node is the current possible state of the system and each edge is labeled with the probability that a newly inserted ball will move the system from one state to another. Then the probability of being in any one state I compute by summing the probabilities of each path into that state.

Is this correct please? I don't trust myself at all when it comes to stats ever since that annoying Monty Hall thing.

Community
  • 1
John Bent
  • 11
  • 3
  • It might be worthwhile noticing that your diagrams implicitly describe a probability distribution over the set of all possible states, one in which the balls are filled sequentially. That is *not* the same thing as the uniform distribution over all allowable distributions of the balls in the bins! Thus, the very first thing to decide is *what distribution are you trying to create?* – whuber Jun 24 '16 at 00:36
  • Thanks @whuber! Unfortunately I either don't understand your reply or I've poorly written the question. I'm not filling the buckets sequentially. I'm just merging states which seem discrepant but are actually identical, just the balls and buckets are indistinguishable. For example, the initial insert makes actually three different states: [1,0,0] or [0,1,0] or [0,0,1]. But since balls and buckets are indistinguishable, I combine all three of these and use the [1,0,0] notation. Is that the misunderstanding? I don't know what distribution I'm trying to create. – John Bent Jun 24 '16 at 03:33
  • Check http://stats.stackexchange.com/questions/204254/how-many-ways-there-are-to-add-up-7-numbers-to-sum-of-180/204258#204258 – Tim Jun 24 '16 at 13:07
  • Compare your question to the one at http://stats.stackexchange.com/questions/184348 . – whuber Jun 24 '16 at 13:43

2 Answers2

1

The easiest way to simulate the system is directly -- draw a bin for each ball. If the bin is full, discard the ball.

If there is no bin constraint, this reduces to a uniform discrete distribution.

Now, you don't particularly care if 2 balls are in bin 0 or in bin 2 -- it seems that you are trying to compute the distribution of the bin counts

E.g. if, for parameters (3,3,2), you get [2,0,1], your bin counts are 0=1 1=1 2=1

My approach was to simply use a monte carlo type simulation and to average the bin count distribution over 10000 runs, as shown below

#####################
# Constrained N-nomial simulator
# Copyright 2016 Gregory R. Bronner

import numpy as np
RUNS=10000
SEED=0
import random


def compute_n_nomial_probabilities(b, B,z):
    if b>B*z:
        raise Exception("No place to put the extra balls")


    #probability distribution...
    results=np.zeros(z+1, dtype=int)


    for r in xrange(RUNS):
        arr=compute_one_bin(b,B,z)
        #now do the counts by bin...
        results+=np.bincount(arr,minlength=z+1)

    #now divide...
    results=results.astype(np.double)/float(RUNS*B)
    return results



def compute_one_bin(b,B,z):
    arr=np.zeros(B, dtype=np.int)
    #compute a bunch of random integers
    rnums=np.random.randint(0,B,10*b)
    succ=0
    for val in rnums:
        if arr[val]<z:
            succ+=1
            arr[val]+=1
            if succ==b:
                return arr

    ### on the off chance that we didn't use up all of our balls, keep simulating...

    while(True):
        val=np.random.randomint(0,B)
        succ += 1
        arr[val] += 1
        if succ == b:
            return arr


if __name__=="__main__":
    random.seed(SEED)
    print compute_n_nomial_probabilities(3,3,3)
    print compute_n_nomial_probabilities(3,3,2)
gbronner
  • 111
  • 3
1

To answer my own question, yes, both probability trees as shown are correct assuming we trust several hundred thousand repetitions of the following random simulator:

#! /usr/bin/env python

from __future__ import division
import random

def simulate(num_bucks,items,bsize,iterations=50000):
  perms = dict()
  for n in range(iterations):
    buckets = [0] * num_bucks 
    for i in range(items):
      while True:
        b = random.randint(0,num_bucks-1)
        if buckets[b] < bsize:
          break  # kludge, loop until we find an unfilled bucket 
      buckets[b] +=1
    buckets.sort()
    buckets = tuple(reversed(buckets))
    try:
      perms[buckets]['count'] += 1
    except KeyError:
      perms[buckets] = {'perm' : buckets, 'count' : 1}

  for perm in perms.values(): 
    perm['prob'] = perm['count'] / iterations
  return perms


def main():
  perms =    simulate(num_bucks=3,items=3,bsize=2)
  for perm in perms.values():
    print perm

if __name__ == "__main__":
    main()

Has output like:

(1, 1, 1) 0.22394
(2, 1, 0) 0.77606
John Bent
  • 11
  • 3