An Urn Experiment: Part 2

In the last post we developed some code to simulate the experiment of choosing marbles from an urn containing r red and g green marbles. We choose one marble at a time and without replacement until we see a red marble. In this post, we will build a function called simUrn that will repeat this experiment many times and return all the values of the random variable. Recall that the random variable is the number of marbles we had to draw before seeing a red marble.

Here is the code we have so far:

rv = 0;
(*build urn*)
reds = Table["R", {3}];
greens = Table["G", {4}];
urn = Join[reds, greens];
(*draw*)
draw = RandomChoice[urn];
rv++;
While[draw === "G",
(*build*)
reds = Table["R", {3}];
greens = Table["G", {4 – rv}];
urn = Join[reds, greens];
(*draw*)
draw = RandomChoice[urn];
(*draw*)
rv++

];
rv

We will be using a Module so that we can localize all the variables. We also need to introduce a new variable called rvList that will store the value of the random variable after each simulation. Finally, we will use a counter variable called counter that will repeat the experiment until the desired number of simulations is achieved.

Here are two templates for Module:

ModuleFunctionTemplates

Both template require two arguments. The first argument is a list containing the local variables and the second argument contains the main body of the function. In the second template, we can initialize the values of the local variables. The code we developed last time (with some minor adjustments) will be the main body of the Module. We need to add a While loop to repeat the experiment until a specified number of times and after each simulation we want to store the value of the random variable.

WhileFunctionTemplate.png

Let’s say we want to simulate the experiment 1000 times, our “test” will be counter < 1000. Our body will be mostly the code above that actually does the simulation.

Here is a high-level view of the Module code:

Module[{counter = 0, rv, reds, greens, urn, draw, rvList={}},]

While[counter < 1000,

(*simulate the experiment*)

(*store the value of the random variable*)

];

(*return the list containing the random variable values*)
rvList
]

Notice that two of the local variables are initialized (counter and rvList) while the others aren’t. It is quite often that the two templates of Module are used together.

Let’s now talk about how to build functions. There are two ways to extend the Wolfram Language. One way is to Set a symbol with an expression or a number: x=RandomReal[]. Another way is to use SetDelayed: x:=RandomReal[]. The former method evaluates the right hand side of the equal sign and assigns the result to x. The latter version only evaluates the right hand side when the symbol x is called. We want the latter version.

Finally our function simUrn needs to be customizable. We want to be able to specify a number of parameters (arguments). Number of reds, number of greens, and the number of simulations. Let’s call these values r, g, counterMax respectively. We will define our function as follows:

simUrn[r_,g_,counterMax_]:=Module[…]

Ok, so let’s put all the code together.

simUrn[r_,g_,counterMax_]:=Module[{counter = 0, rv, reds, greens, urn, draw, rvList={}},

While[counter < 1000,

(*simulate the experiment*)
rv = 0;
(*build urn*)
reds = Table["R", {3}];
greens = Table["G", {4}];
urn = Join[reds, greens];
(*draw*)
draw = RandomChoice[urn];
rv++;
While[draw === "G",
(*build*)
reds = Table["R", {3}];
greens = Table["G", {4 – rv}];
urn = Join[reds, greens];
(*draw*)
draw = RandomChoice[urn];
(*draw*)
rv++
]
(*store the value of the random variable*)
AppendTo[rvList,rv]
];
(*return the list containing the random variable values*)
rvList
]

We are almost done. Now we need to make sure we pass on the arguments (number of reds, number of greens, etc) to the code. Currently, the code will run with 3 red and 4 green marbles and will simulate the experiment 1000 times. Here is the adjusted code.

simUrn[r_, g_, counterMax_] :=
Module[{counter = 0, rv, reds, greens, urn, draw, rvList = {}},
While[counter < counterMax,(*simulate the experiment*)rv = 0;
(*build urn*)reds = Table["R", {r}];
greens = Table["G", {g}];
urn = Join[reds, greens];
(*draw*)draw = RandomChoice[urn];
rv++;
While[draw === "G",
(*build*)reds = Table["R", {r}];
greens = Table["G", {g – rv}];
urn = Join[reds, greens];
(*draw*)
draw = RandomChoice[urn];
rv++]
(*store the value of the random variable*)
AppendTo[rvList, rv];
counter++;
];
(*return the list containing the random variable values*)
rvList]

Let’s test this code by running the function simUrn[3, 4, 100]. My output is shown below. You output will contain different numbers.

{2,2,2,3,3,2,4,2,3,1,1,3,1,1,1,4,1,4,2,1,4,3,1,3,2,4,2,1,2,4,1,2,2,2,3,1,2,2,4,3,2,3,3,1,3,3,1,3,2,1,3,2,1,1,2,1,1,3,2,3,1,1,1,2,4,3,1,1,3,2,3,3,4,4,4,1,5,4,5,5,1,3,1,2,1,1,1,1,1,4,3,3,1,1,2,4,1,4,4,3}

We will now wrap the code inside Histogram to a get visual of the distribution. Here are some templates on how to use this built-in function.
HistogramFunctionTemplates.png

We will use the last version, use Automatic as the bin specification “bspec” and “Probability” as the height specification “hSpec”. Instead of copy-pasting the list that we got, we can just put the code used to generate the list of values. Here is the code that will simulate the experiment 1000 times and show a histogram showing the probabilities associated with each value of the random variable.

Histogram[simUrn[3, 4, 1000],Automatic, "Probability"]

Histogram1.png

Leave a Reply

Your email address will not be published. Required fields are marked *