brisela.py

#!/usr/bin/env python

from random import shuffle
import numpy as np

nb = 2 # Number of Bruna's in deck
ng = 3 # Number of Gisela's in deck
ncards = 60 # Total number of cards in deck

# Construct the deck
deck = ['b' for i in range(nb)]
deck.extend(['g' for i in range(ng)])
deck.extend(['l' for i in range(ncards - nb - ng)])

ndraws = np.zeros(ncards, dtype=int)

# Shuffle it
for n in range(10000):
    found_b = False
    found_g = False
    shuffle(deck)
    for i, c in enumerate(deck):
        if c == 'b':
            found_b = True
        elif c == 'g':
            found_g = True
        if found_b and found_g:
            ndraws[i] += 1
            break

prob = np.array(ndraws, dtype=float) / float(np.sum(ndraws))
print(np.cumsum(prob))

Output for nb = 2 and ng = 3:

[ 0.      0.0031  0.0098  0.0181  0.0322  0.0465  0.0618  0.08    0.0985
  0.118   0.1439  0.1683  0.196   0.2223  0.2489  0.2781  0.3053  0.3344
  0.3624  0.3919  0.4221  0.4508  0.4774  0.5043  0.5286  0.5543  0.5807
  0.6091  0.6368  0.6588  0.6813  0.7045  0.7287  0.7512  0.773   0.7907
  0.8091  0.8303  0.8492  0.8656  0.8808  0.8951  0.9091  0.9202  0.9316
  0.9418  0.9508  0.958   0.9654  0.9716  0.9779  0.9835  0.9875  0.9917
  0.9947  0.9969  0.9978  0.9994  1.      1.    ]

Output for nb = 4 and ng = 4:

[ 0.      0.0086  0.027   0.0483  0.0761  0.1084  0.1445  0.1819  0.2235
  0.2662  0.306   0.3489  0.386   0.4295  0.4704  0.5122  0.5484  0.5837
  0.6155  0.654   0.6845  0.7172  0.747   0.7704  0.7961  0.8164  0.8363
  0.8546  0.8721  0.8857  0.9024  0.9155  0.9275  0.9381  0.9471  0.9546
  0.9606  0.9674  0.9728  0.9777  0.9825  0.9865  0.9893  0.9922  0.9946
  0.9958  0.9969  0.9981  0.9987  0.9991  0.9994  0.9995  0.9999  1.      1.
  1.      1.      1.      1.      1.    ]