Lets say we have an urn with 100 balls (20 red, 30 blue, 50, green). Whenever a ball is dawn there is a chance it will be kept and not put it back. This chance with which the ball is kept is different for each colour (red = 0.5, blue = 0.3, green = 0.2). So for example the first ball randomly drawn is blue and with a 0.3 chance it is kept otherwise it is put back into the urn and another random ball is drawn. This is repeated until 50 balls have been kept. What is the most likely distribution of coloured balls/expected number of balls of each colour of the 50 kept?

The application I need it for is because I am simulating animals eating a certain amount of total food (the 50 balls) made up of different plant types (different colours) that are available differently in the landscape (number of balls per colour and number of total balls) and that are eaten with a certain preference (probability to keep a ball).

I have calculated this by simulating this draw 10000 times in the programming language R and calculating the mean number of balls of each colour. Code is below in case you are interested but it is not needed to answer my question. For 50 balls taken the most likely distribution is roughly:

  • Red: 14.3367 balls
  • Blue: 15.8905 balls
  • Green: 19.7728 balls

However I do not want to simulate it but just be able to calculate it without simulating this draw as it takes much longer especially if I have more than 3 colours (or plant types). I would very much appreciate your answer or if you can tell me where to find one or how this type of problem is called, as I was unsucessful in googling this problem. Thank you in advance!

#The find_prob is number of balls in the urn divided by total balls and the keep_prob is the chance the ball will be kept. 

blue_find_prob<-red_find_prob + total_blues/total_balls




for (i in 1:10000) { # repeated 1000 times to get closer to the real number

while ((red+blue+green < needed)) {
  blue_find_prob<-red_find_prob + total_blues/total_balls
  colour_drawn<-runif(1,min = 0, max = 1)
  keep_colour<-runif(1,min = 0, max = 1)
  if (colour_drawn < red_find_prob){
    if (keep_colour < red_keep_prob) {
      red = red + 1
      total_balls = total_balls - 1
      total_reds = total_reds -1
    }else {}
  }else if (colour_drawn < blue_find_prob) {
    if (keep_colour < blue_keep_prob) {
      blue = blue + 1
      total_balls = total_balls - 1
      total_blues = total_blues - 1
    }else {}
  }else {
    if (keep_colour < green_keep_prob) {
      green = green + 1
      total_balls = total_balls - 1
    }else {}

Just to make sure it is not misunderstood: I do not need help with the code I need a mathematical formula to calculate the most likely distribution (of course if you find a griveous mistake in the code please point it out).

EDIT: I have finally found a similar question not sure if it exactly what I need, I'll have a closer look once I am back home: Drawing balls from a bin with a color-specific probabilistic discard step

  • 1
    $\begingroup$ When you say "the most likely distribution of coloured balls of the 50 drawn", do you mean "the most likely distribution of coloured balls of the 50 balls that are kept"? $\endgroup$
    – Arthur
    Nov 5 '20 at 0:11
  • $\begingroup$ You surely don't mean "the most likely distribution of colors", since the distribution you give is impossible; there can't be $14.3367$ red balls, for example. Do you perhaps mean the expected number of balls of each color? $\endgroup$
    – saulspatz
    Nov 5 '20 at 0:58
  • $\begingroup$ Yes, thank you I have edited the post to make it clearer. $\endgroup$
    – Alex
    Nov 5 '20 at 1:21
  • $\begingroup$ +1 fascinating question, which I am clueless on. $\endgroup$ Nov 5 '20 at 1:24
  • $\begingroup$ Are you familiar with Markov chains? You don't need to use them to actually compute the answer, but it's easier to discuss the problem if you're familiar with them. Otherwise, I'll tell you what you need to know in my answer. $\endgroup$
    – saulspatz
    Nov 5 '20 at 13:08

I don't think there's a closed-form formula that will solve this problem in any reasonable way. There is a procedure for getting the exact solution, but it's harder to program than simulation.

I will talk about computing the expected number of balls of each color remaining in the urn once $50$ balls have been kept. If we know that, we can at once compute the expected number of balls of each color that are outside the urn.

We can describe the situation at any point by and ordered triple $(r,g,b)$, where $r$ is the number of red balls in the urn, $b$ is the number of blue balls, and $g$ the number of green balls. We call this the "state" of the system. If we are in the state $(r,g,b)$ and we draw a ball, the next state will be one of $$(r,g,b)\\(r-1,g,b)\\(r,g-1,b)\\(r,g,b-1)$$ and the probability of transitioning to each of these states is easily computed. Of course, if $r+g+b=50$, then the process is over, and finishes in this state. The states with $r+g+b=50$ are called "absorbing". All other states are called "transient".

The salient point about this process is that the probability of transitioning from one state to another depends only on the current state, not how we got there. We say that the process has "no memory." These facts mean that we are dealing with a finite-state, absorbing Markov chain. I'll refer to the wiki article to explain what I'm saying, but it's not necessary for you to master all the theory to follow what I'm saying.

Let $A$ be the set of absorbing states. For $(r,g,b)\in A$ let $P(r,g,b)$ be the probability that the process stops in state $(r,g,b)$. The expected number of red balls remaining in the urn is $$\sum_{(r,g,b)\in A}rP(r,g,b)$$ and similar formulas hold for the other colors.

We first to get an idea of the scale of this problem. Brute force calculation shows that there are $651$ absorbing states and $16,275$ transient states. The matrix $Q$ in the wiki article, which encodes, the transition probabilities among the transient states, is $16275\times16275$ as in the fundamental matrix $N=(I-Q)^{-1}$ We don't want to store $Q$ in memory, or to compute its inverse. We can avoid this, because of a special property of this particular Markov chain.

Once the chain has left a particular state, it can never return to it. Let us number the states in such a way that if state $s_1$ has a higher number than state $s_2$ then it is not possible to go from $s_1$ to $s_2$. We can accomplish this by sorting a list of all states in order of decreasing values of $r+g+b$. (It doesn't matter how balls with the same number of balls are sorted among themselves. We can use dictionary order for example.) Then we number each state with its index in the sorted list. What this means is that the matrix $Q$ will be upper triangular, which greatly simplifies the calculations. We will need the matrix $I-Q$ . Let $U=I-Q$, where $U$ stands for "upper."

Now on the wiki under "Absorbing Probabilities" it says

the probability of being absorbed in the absorbing state $j$ when starting from transient state $i$, which is the $(i,j)$-entry of the matrix $$B=NR.$$ $N$ is defined as $U^{-1}$ and since we don't want to compute inverses, we rewrite this as $$UB=R\tag1$$ We only care about the probability of ending in state $j$ when $i$ is the initial state, so we only care about the first row of $B$. So far as I can see, unfortunately, there's no way to compute the first row of $B$ without computing the other rows, but we don't have to save the values.

In equation $(1)$, $U$ is $16,275\times16,275$, and $B$ and $R$ are $16,625\times651$. We can solve for $B$ column-by-column, meaning that we have $651$ systems of linear equations, each of which consists of $16,275$ equations in $16,275$ unknowns.

Since we have arranged for $U$ to be upper-triangular, each such system is easily solved by back-substitution. Once we have gotten the probability $P(r,g,b)$ we can increment the expectations: $$E_r=E_r+rP(r,g,b)\\E_b=E_b+bP(r,g,b)\\E_g=E_g+gP(r,g,b)$$ and we don't need the solution to that system any more, so we can reuse the space.

Also, notice that $U$ is a very sparse matrix. At most $4$ entries in any row or column are non-zero. That should make it possible to program the back-substitution to be very efficient.

This would be all very well, if you just had just this one specific problem, but your question suggests that you may have many of them, with differing numbers of colors. I'm not sure how easy it would be do automate this solution, so it would handle any problem of this type automatically, (subject to size constraints, of course.)

I'll have to think about that.


I wrote a python program to perform the calculations described above. (Sorry, I don't know R).

from itertools import product 

balls = (20,30,50)
probs = (.5,.3,.2)  # probability that ball will be kept
kept = 50
N = sum(balls)-kept
M= len(balls)
expect = list(balls)

absorbing = []
transient = []

for state in product(*[range(b+1) for b in balls]):
    if sum(state)== N:
    elif sum(state)>N:
def successors(state):
    answer = []
    for idx, c in enumerate(state):
        if c >0:
    return answer

def from2(s,t):
    # probability of moving from s to t in one move
    # Pre: t is a successor of s
    i = {i for i in range(M) if s[i] !=t[i]}.pop()
    return probs[i]*s[i]/sum(s)
# sort the transient states in decreasing order
# of the number of balls.  Make an inverted list
# from state to index.

transient.sort(key=sum, reverse=True)
tran2idx = {s:i for i,s in enumerate(transient)}

# Q is the transition matrix for transient states to transients states.
# U is I-Q
# R is is the transition matrix for transient states to absorbing states
# In the equation UB = R, B[i,j] is the probability of reaching absorbing 
# state j from transient state i.  Sorting the transient states  assures that
# U is upper triangular, so we can solve by back substitution.
# We do this column-by-column

b = len(transient)*[0]  # current column of B

for j, a in enumerate(absorbing):   
    for t in reversed(transient):
        # p is probability of moving out of state t
        p =sum(s*probs[i] for i,s in enumerate(t))/sum(t)
        r = 0
        for s in successors(t):
            if s == a:
                r += from2(t,s)
            elif sum(s) > N:
                k = tran2idx[s]
                r += b[k]*from2(t,s) 
        i = tran2idx[t]
        b[i] = r/p
    for i in range(M):
        expect[i] -= a[i] * b[0]
for i in range(M):
    print(balls[i], probs[i], expect[i])

This produced

20 0.5 14.325163888283411
30 0.3 15.908709486720031
50 0.2 19.76612662499662 

in close agreement with your simulation results.

This took $3$ minutes to run on my old $2012$ Mac mini. If you want to do it for more balls and more colors, the number of states will explode. You could program it in C, and the problem is a natural for parallelization, since we deal with each column independently, but I think that simulation is the way to go.

The fact that your simulation produced results so close to the theoretical answers, should encourage you to use simulation, I would think.

  • $\begingroup$ If the answer solves your problem, please accept it by clicking the check mark. $\endgroup$
    – saulspatz
    Nov 9 '20 at 14:01
  • $\begingroup$ Thank you very much for your answer by the way, it is very satisfying that I have the theoretical background to my approach and of course that they are the same : D $\endgroup$
    – Alex
    Jun 9 at 22:04

This is a more knowledgeable you from the future (wow!) Here is a piece of code that is way faster (at least in R) because it uses vector calculation. It is a deterministic approach and just an approximation (although quite good)

The basis of how it works is that it just takes the find probability times the keep_probabilty. This number is taken times an "amount_removed". The closer you set that amount removed to 0, the more accurate it will be. Then that amount is removed and the find probabilities are updated. The theory behind it is that normally you couldnt just take find probability times keep probability times total needed, because it does not take into account that once you remove something you lower the probability for that to be found. However, if the amount you remove goes close to 0 then it does take it into account the next time you remove a very small amount as you have updated the find probability.

Results in comparison to the approach in the question and saulspatz's answer:

Amount keep_chance question approach saulspatz's answer this answer
20 0.5 14.3367 14.325163888283411 14.32446
30 0.3 15.8905 15.908709486720031 15.90252
50 0.2 19.7728 19.76612662499662 19.77302

Code in R:

#Your input

#thinks for the algorithm
max_amount_removed<-needed/100 #the lower you set max_amount_removed the more accurate but also longer it will take

while(needed > max_amount_removed) {
  #this is to go with smaller steps towards the end to get a lot closer to the total needed
  if (needed < 1.5 * max_amount_removed && counter_max_amount_removed < 50) { 
    max_amount_removed<-max_amount_removed / 2
    counter_max_amount_removed <- counter_max_amount_removed + 1
  taken<-max_amount_removed*(find_prob*preference) #this is the main step
  counter<-counter + 1
print(counter) #shows how many iterations it took

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.