Monte Carlo Simulation

Question

I'm a student in a Java Programming class. My problem deals with an interpretation of the Monte Carlo Simulation. I'm supposed to find the probability that three quarters or three pennies will be picked out of a purse that has 3 quarters and 3 pennies. Once a coin is picked it is not replaced. The probability should be 0.1XXXXXXX. I keep getting 0 or 1 for my answer. This is what i have so far.

public class CoinPurse {
    public static void main(String[] args) {
        System.out.print("Probability of Drawing 3 coins of the Same Type - ");
        System.out.println(coinPurseSimulation(100));
    }

    /**
     Runs numTrials trials of a Monte Carlo simulation of drawing 
     3 coins out of a purse containing 3 pennies and 3 quarters. 
     Coins are not replaced once drawn.
     @param numTrials - the number of times the method will attempt to draw 3 coins
     @returns a double - the fraction of times 3 coins of the same type were drawn.
     */

    public static double coinPurseSimulation(int numTrials) {
        final int P = 1;
        final int Q = 2;
        int [] purse = {Q, Q, Q, P, P, P};
        int [] drawCoins = new int[3];
        for (int draw = 0; draw < 3; draw ++) {
            int index = (int)(Math.random() * purse.length);
            drawCoins[draw] = purse[index];
            int [] newPurse = new int[purse.length-1];
            int j = 0;
            for (int i =0; i < purse.length; i++) {
                if (i == index) {
                    continue;
                }
                newPurse[j] = purse[i];
                j++;
            }
            purse = newPurse;
        }
        double number = 0.0;
        double result = 0.0;
        for (int i = 0; i < numTrials; i++) {
            result++;
            for (int j = 0; j < numTrials;j++) {
                if(drawCoins[0] == drawCoins [1] && drawCoins[1] == drawCoins[2]) {
                    number++;
                }
            }
        }
        return number/result;
    }
}

Answer 1

The reason you only ever get 0 or 1 is that you only draw (or pick) coins from the purse once, but you then test that draw numTrials * numTrials times. You have two loops (with indices i and j) iterating numTrials time - your logic is a little messed up there.

You can put the first loop (for drawing coins) within a second loop (for running trials) and your code will work. I've put a minimal refactor below (using your code as closely as possible), with two comments afterwards that might help you some more.

public class CoinPurse
{
    public static void main(String[] args)
    {
        System.out.print("Probability of Drawing 3 coins of the Same Type - ");
        System.out.println(coinPurseSimulation(100));
    }

    /**
     * Runs numTrials trials of a Monte Carlo simulation of drawing 3 coins out
     * of a purse containing 3 pennies and 3 quarters. Coins are not replaced
     * once drawn.
     * 
     * @param numTrials
     *            - the number of times the method will attempt to draw 3 coins
     * @returns a double - the fraction of times 3 coins of the same type were
     *          drawn.
     */

    public static double coinPurseSimulation(int numTrials)
    {
        final int P = 1;
        final int Q = 2;

        double number = 0.0;
        double result = 0.0;

        // Changed your loop index to t to avoid conflict with i in your draw
        // loop
        for (int t = 0; t < numTrials; t++)
        {
            result++;

            // Moved your draw without replacement code here
            int[] purse =
            { Q, Q, Q, P, P, P };
            int[] drawCoins = new int[3];

            for (int draw = 0; draw < 3; draw++)
            {
                int index = (int) (Math.random() * purse.length);
                drawCoins[draw] = purse[index];
                int[] newPurse = new int[purse.length - 1];
                int j = 0;
                for (int i = 0; i < purse.length; i++)
                {
                    if (i == index)
                    {
                        continue;
                    }
                    newPurse[j] = purse[i];
                    j++;
                }
                purse = newPurse;
            }

            // Deleted the loop with index j - you don't need to test the same
            // combination numTrials times...
            if (drawCoins[0] == drawCoins[1] && drawCoins[1] == drawCoins[2])
            {
                number++;
            }
        }

        return number / result;
    }
}

Picking coins code

I have some comments on your routing for drawing coins:

1. It works correctly
2. It is rather cumbersome
3. It would have been easier for you to spot the problem if you had broken this bit of code into a separate method.

I'm going to address 3 and then 2.

Break the drawing code out into a method

private static int[] pickCoins(int[] purse, int numPicks)
{
    //A little error check
    if (numPicks > purse.length)
    {
        System.err.println("Can't pick " + numPicks + 
                           " coins from a purse with only " + purse.length + " coins!");
    }

    int[] samples = new int[numPicks];

    // Your sampling code here

    return samples;
}

You can now simply call from within your second loop i.e.

drawCoins = pickCoins(purse, 3);

Sampling algorithm

@pjs's answer recommends using Collections.shuffle() and then taking the first 3 coins in your collection (e.g. an ArrayList). This is a good suggestion, but I'm guessing you haven't been introduced to Collections yet, and may not be 'allowed' to use them. If you are - do use them. If not (as I assume), you might want to think about better ways to randomly draw n items from an r length array without replacement.

One (well accepted) way is the Fisher-Yates shuffle and its derivatives. In effect it involves picking randomly from the unpicked subset of an array.

In Java - an working example could be as follows - it works by moving picked coins to the "end" of the purse and picking only from the first maxInd unpicked coins.

    private static int[] pickCoins(int[] purse, int numCoins)
    {
        int[] samples = new int[numCoins];
        int maxInd = purse.length - 1;

        for (int i = 0; i < numCoins; i++)
        {
            int index = (int) (Math.random() * maxInd);
            int draw = purse[index];
            samples[i] = draw;
            // swap the already drawn sample with the one at maxInd and decrement maxInd
            purse[index] = purse[maxInd];
            purse[maxInd] = draw;
            maxInd -= 1;
        }

        return samples;
    }

Expected results

You say your expected result is 0.1XXXXXXX. As you're learning Monte Carlo simulation - you might need to think about that a little more. The expected result depends on how many trials you do.

First, in this simple example, you can consider the analytic (or in some sense exact) result. Consider the procedure:

1. You draw your first coin - it doesn't matter which one it is
2. Whichever coin it was, there are 2 left in the bag that are the same - the probability of picking one of those is 2 / 5
3. If you picked one of the matching coins in step 2, there is now 1 matching coin left in the bag. The probability of picking that is 1 / 4
4. So, the probability of getting 3 matching coins (of either denomination) is 2 / 5 * 1 / 4 == 2 / 20 == 0.1

Your Monte Carlo programme is trying to estimate that probability. You would expect it to converge on 0.1 given sufficient estimates (i.e. with numTrials high enough). It won't always give a value equal to, or even starting with, 0.1. With sufficient number of trials, it's likely to give something starting 0.09 or 0.1. However, if numTrials == 1, it will give either 0 or 1, because it will draw once and the draw will either match or not. If numTrials == 2, the results can only be 0, 0.5 or 1 and so on.

One of the lessons of doing Monte Carlo simulation to estimate probabilities is to have a sufficiently high sample count to get a good estimate. That in turn depends on the accuracy you want - you can use your code to investigate this once it's working.

Answer 2

You need to move the loop where you generate draws down into the numTrials loop. The way you've written it you're making a single draw, and then checking that one result numTrials times.

I haven't checked the logic for your draw carefully, but that's because I'd recommend a different (and much simpler) approach. Use Collections.shuffle() on your set of quarters and pennies, and check the first three elements after each shuffle.

If done correctly, the answer should be 2 * (3/6) * (2/5) * (1/4), which is 0.1.