Simple solution to coloured marble problem seems off?

I have a very basic marble problem (paraphrased):

There is a vase with 20 balls, of which 7 are red, 3 are blue, 1 is orange, and 9 are green. What are the odds of drawing a red, a blue, an orange and a green marble from the vase in this specific order ( $P(r b o g)$ )? We draw four marbles without putting them back.

My solution to this problem was straightforward. Using the conventional definition of a probability: $P(A) = \frac{|A|}{|\Omega|}$ if $\Omega$ is the set of all possible outcomes.

I then defined the event of drawing a red, a blue, an orange, a green marble as such. I labeled my marbles 1 through 20 and labelled the red ones 1 through 7, the blue ones 8 through 10, the orange one 11 and the blue ones 12 through 20.

If we call the event of interest $A$, then the set description of $A$ becomes:

$$A = \{(\omega_1, \omega_2, \omega_3, \omega_4), \omega_1 \in \{1, 2, ..., 7\}, \omega_2 \in \{8, 9, 10\}, \omega_3 = 11, \omega_4 \in \{12, 13, ..., 20\}\}$$

If we call the set of all possible outcomes $\Omega$, this set becomes:

$$\Omega = \{(\omega_1, \omega_2, \omega_3, \omega_4), \omega_i\in\{1,2,...,20\},\omega_k\neq\omega_j, k \neq j\}$$

The sizes of the sets are easily calculated as $|A| = 7*3*1*9 = 189$ and $|\Omega| = 20*19*18*17=116280$. With this, the probability of the event becomes $\frac{189}{116280} = 0.001625387$

Subsequently, I tried to check my work by doing a few simulations. I did three:

1. A simulation in Python. I would generate a list with elements 1 through 20. Then I would let Python pick an item from the list and remove it four times. Thus, this would simulate drawing four marbles and not putting them back. I then checked whether I drew a red-blue-orange-green combo.

2. A simulation in C++. I would draw four marbles at random (1-20) and only use the result of the random draw if it made up four different numbers.

3. A simulation in C++. I would shuffle a list with elements 1 through 20 and then use the first four elements.

Simulation 2 confirms my result. Simulations 1 and 3 do not: they give probability 0.00158... when I let them iterate 100,000,000 times or more.

This seems strange, since I'm so sure of my maths. So the question becomes: is my calculation of probability incorrect, or my way of simulation?

• Just go step by step: $\frac 7{20}\times \frac 3{19}\times \frac 1{18}\times \frac 9{17}\approx 0.001625387$ just as you said. All three simulation methods looks sensible, but of course bugs are possible.
– lulu
Feb 11 '18 at 15:25
• As a way to check for bugs, use your various simulations to compute other things. First marble red, that sort of thing. See if you are getting the right answers for other computations.
– lulu
Feb 11 '18 at 15:31
• Thanks for the remark. I decided to indeed do a bug search. What I did is simulate the first four draws using a completely random trial-and-error method to obtain four individual integers. Then I simulated them using the C++ shuffle() function on an array of 20 integers. I collected data on 1,000,000 experiments in a .csv and plotted histograms of the experiments. I expected that all 20 integers would be drawn equally. However, there was an immense bias towards higher numbers when using the shuffle() function! So there was indeed a bug. Feb 11 '18 at 20:00