Chapter 2: urn probabilities #1

To program the simulation, we first define the number of repetitions, the number of balls in the urn, and the list of balls (from 1 to the number of balls) in the urn:

REPS=10000; %number of repetitions in the simulation

Nballs=118; %number of balls

URN=1:Nballs; %define the urn [first 54 are red]

We next make a matrix of NaN (‘not a number’) values whose size is equal to the number of simulation repetitions x the number of observations per repetition:

D=nan(REPS,2); %‘empty’ observation matrix

Finally, we write a for-loop that executes the simulation (fills the matrix with simulated observations, replacing nan’s), and find the simulation runs that meet the criterion that both observations are red (i.e., less than or equal to 54). The length of the list (rows in D) found to match this criterion, divided by the total number of simulation runs (REPS), is our simulation-based approximation to the desired probability calculation.

for n=1:REPS, Unow=URN(randperm(Nballs)); D(n,:)=Unow(1:2); end

sum(and(D(:,1)<=54,D(:,2)<=54))/REPS

Or, we could write a version of the loop that does not explicitly simulate the mixing of the balls using the RANDPERM command, instead using the RANDI command:

for n=1:REPS, D(n,1)=randi(Nballs); unow=exclude(URN,D(n,1));

D(n,2)=unow(randi(Nballs-1)); end

sum(and(D(:,1)<=54,D(:,2)<=54))/REPS