Monte Carlo Methods Using the KSL

Learning Objectives

• To be able to perform simple Monte Carlo simulation experiments using the KSL
• To review statistical concepts and apply them to the analysis of simple Monte Carlo simulation experiments
• To illustrate generating and collecting statistics for Monte Carlo simulation experiments

Simple Monte Carlo Integration

Suppose we have some function, \(g(x)\), defined over the range \(a \leq x \leq b\) and we want to evaluate the integral: \[ \theta = \int\limits_{a}^{b} g(x) \mathrm{d}x\] Define \(Y\) as follows with \(X \sim U(a,b)\): \[ Y = \left(b-a\right)g(X) \] By definition,

\[ E_{X}\lbrack g(X) \rbrack = \int\limits_{a}^{b} g(x)f_{X}(x)\mathrm{d}x\]

Monte Carlo Estimator for \(\hat{\theta}\)

A good estimator for \(E\lbrack Y \rbrack\) is the sample average of observations of \(Y\). Let \(Y_{1}, Y_{2},...Y_{n}\) be a random sample of observations of \(Y\). Let \(X_{i}\) be the \(i^{th}\) observation of \(X\). Substituting each \(X_{i}\) yields the \(i^{th}\) observation of \(Y\),

\[Y_{i} = \left(b-a\right)g(X_{i})\]

Then, the sample average of is:

\[\begin{aligned} \bar{Y}(n) & = \frac{1}{n}\sum\limits_{i=1}^{n} Y_{i} = \left(b-a\right)\frac{1}{n}\sum\limits_{i=1}^{n}\left(b-a\right)g(X_{i})\\ & = \left(b-a\right)\frac{1}{n}\sum\limits_{i=1}^{n}g(X_{i}) = \hat{\theta} \\\end{aligned}\]

where \(X_{i} \sim U(a,b)\).

Example

Suppose that we want to estimate the area under \(f(x) = x^{\frac{1}{2}}\) over the range from \(1\) to \(4\). That is, we want to evaluate the integral:

\[\theta = \int\limits_{1}^{4} x^{\frac{1}{2}}\mathrm{d}x = \dfrac{14}{3}=4.6\bar{6}\] According to the previously presented theory, we need to generate \(X_i \sim U(1,4)\) and then compute \(\bar{Y}\), where \(Y_i = (4-1)\sqrt{X{_i}}= 3\sqrt{X{_i}}\).

KSL Code for Example

fun main() {
    val a = 1.0
    val b = 4.0
    val ucdf = UniformRV(a, b, streamNum = 3)
    val stat = Statistic("Area Estimator")
    val n = 100 // sample size
    for (i in 1..n) {
        val x = ucdf.value
        val gx = sqrt(x)
        val y = (b - a) * gx
        stat.collect(y)
    }
    System.out.printf("True Area = %10.3f %n", 14.0 / 3.0)
    System.out.printf("Area estimate = %10.3f %n", stat.average)
    println("Confidence Interval")
    println(stat.confidenceInterval)
}

KSL Results for Example

True Area =      4.667
Area estimate =      4.781
Confidence Interval
[4.608646560421988, 4.952515649272401]

Notice that the estimate is not equal to the true area. This difference is due to sampling error. We can determine the sample size need to bound the error.

Determining the Sample Size

The confidence interval for a point estimator can serve as the basis for determining how many observations to have in the sample. The quantity:

\[ h = t_{1-(\alpha/2), n-1} \dfrac{s}{\sqrt{n}} \]

is called the half-width of the confidence interval. You can place a bound, \(\epsilon\), on the half-width by picking a sample size that satisfies:

\[ h = t_{1-(\alpha/2), n-1} \dfrac{s}{\sqrt{n}} \leq \epsilon \]

We call \(\epsilon\) the margin of error for the bound and \(t_{p, \nu}\) is the \(100p\) percentage point of the Student t-distribution with \(\nu\) degrees of freedom.

Determining the Sample Size Via Normal Approximation

The required sample size can be approximated using the normal distribution. Solving for \(n\) yields:

\[n \geq \left(\dfrac{t_{1-(\alpha/2), n-1} \; s}{\epsilon}\right)^2\]

As \(n\) gets large, \(t_{1-(\alpha/2), n-1}\) converges to the \(100(1-(\alpha/2))\) percentage point of the standard normal distribution \(z_{1-(\alpha/2)}\). This yields the following approximation:

\[ n \geq \left(\dfrac{z_{1-(\alpha/2)} s}{\epsilon}\right)^2 \]

Using Statistic Class to Compute the Sample Size

The normal approximation is implemented in the function, estimateSampleSize of the Statisitic class’s companion object.

/** Estimate the sample size based on a normal approximation
* @param desiredHW the desired half-width (must be bigger than 0)
* @param stdDev    the standard deviation (must be bigger than or equal to 0)
* @param level     the confidence level (must be between 0 and 1)
* @return the estimated sample size
*/
fun estimateSampleSize(desiredHW: Double, stdDev: Double, level: Double): Long {
   require(desiredHW > 0.0) { "The desired half-width must be > 0" }
   require(stdDev >= 0.0) { "The desired std. dev. must be >= 0" }
   require(!(level <= 0.0 || level >= 1.0)) { "Confidence Level must be (0,1)" }
   val a = 1.0 - level
   val a2 = a / 2.0
   val z = Normal.stdNormalInvCDF(1.0 - a2)
   val m = z * stdDev / desiredHW * (z * stdDev / desiredHW)
   return (m + .5).roundToLong()
}

Example

Suppose we are required to estimate the output from a simulation so that we are 99% confidence that we are within \(\pm 0.1\) of the true population mean. After taking a pilot sample of size \(n=10\), we have estimated \(s=6\). What is the required sample size?

For a 99% confidence interval, we have \(\alpha = 0.01\) and \(\alpha/2 = 0.005\). Thus, \(z_{1-0.005} = z_{0.995} = 2.5758293064439264\). Because the margin of error, \(\epsilon\) is \(0.1\), we have that,

\[ n \geq \left(\dfrac{z_{1-(\alpha/2)}\; s}{\epsilon}\right)^2 = \left(\dfrac{2.5758293064439264 \times 6}{0.1}\right)^2 = 23885.63 \approx 23886 \]

KSL Code to Compute Sample Size

The following KSL code will easily compute the sample size.

fun main() {
    val desiredHW = 0.1
    val s0 = 6.0
    val level = 0.99
    val n = Statistic.estimateSampleSize(
        desiredHW = desiredHW,
        stdDev = s0,
        level = level
    )
    println("Sample Size Determination")
    println("desiredHW = $desiredHW")
    println("stdDev = $s0")
    println("Level = $level")
    println("recommended sample size = $n")
}

Sample Size Determination
desiredHW = 0.1
stdDev = 6.0
Level = 0.99
recommended sample size = 23886

Sample Size for a Proportion

If the quantity of interest is a proportion, then a different method can be used. In particular, a \(100\times(1 - \alpha)\%\) large sample two sided confidence interval for a proportion, \(p\), has the following form:

\[ \hat{p} \pm z_{1-(\alpha/2)} \sqrt{\dfrac{\hat{p}(1 - \hat{p})}{n}} \]

where \(\hat{p}\) is the estimate for \(p\). From this, you can determine the sample size via the following equation:

\[ n = \left(\dfrac{z_{1-(\alpha/2)}}{\epsilon}\right)^{2} \hat{p}(1 - \hat{p}) \]

Example

Suppose that we want to simulate a normally distributed random variable, \(X\), with \(E[X] = 10\) and \(Var[X] = 16\). From the simulation, we want to estimate the true population mean with 99 percent confidence for a half-width \(\pm 0.50\) margin of error. In addition to estimating the population mean, we want to estimate the probability that the random variable exceeds 8. That is, we want to estimate, \(p= P(X>8)\). We want to be 95% confident that our estimate of \(P(X>8)\) has a margin of error of \(\pm 0.05\). What are the sample size requirements needed to meet the desired margins of error?

Indicator Variables

Let \(X\) represent the unknown random variable of interest. Then, we are interested in estimating \(E[X]=\mu\) and \(p = P[X > 8]\). To estimate a probability it is useful to define an indicator variable. An indicator variable has the value 1 if the condition associated with the probability is true and has the value 0 if it is false. To estimate \(p\), define the following indicator variable:

\[ Y_i = \begin{cases} 1 & X_i > 8\\ 0 & X_i \leq 8 \\ \end{cases} \] This definition of the indicator variable \(Y_i\) allows \(\bar{Y}\) to be used to estimate \(p\).

\[ \hat{p} = \bar{Y} = \frac{1}{n}\sum_{i=1}^{n}Y_i = \frac{\#\lbrace X_i > 8\rbrace}{n} \]

Example KSL Solution

fun main() {
    val rv = NormalRV(10.0, 16.0, streamNum = 3)
    val estimateX = Statistic("Estimated X")
    val estOfProb = Statistic("Pr(X>8)")
    val r = StatisticReporter(mutableListOf(estOfProb, estimateX))
    val n0 = 20 // sample size
    for (i in 1..n0) {
        val x = rv.value
        estimateX.collect(x)
        estOfProb.collect(x > 8)
    }
    println(r.halfWidthSummaryReport())
}

Example Sample Size Determination

To determine the sample size so that we get a half-width that is less than \(\epsilon = 0.1\) with 99% confidence, we can use the previous equations. Since the table reports a half-width for a 95% confidence interval, we need to solve for \(s\) in terms of \(h\). Rearranging the half-width equation yields,

\[ s = \dfrac{h\sqrt{n}}{t_{1-(\alpha/2), n-1}} \]

Using the results and \(\alpha = 0.05\) yields,

\[ s_0 = \dfrac{h\sqrt{n_0}}{t_{1-(\alpha/2), n_0-1}} = \dfrac{2.2248\sqrt{20}}{t_{0.975, 19}}= \dfrac{2.2248\sqrt{20}}{2.093024}=4.7536998 \]

Example Sample Size Determination

For a 99% confidence interval, we have \(\alpha = 0.01\) and \(\alpha/2 = 0.005\). Thus, \(z_{0.995} = 2.576\). Because the margin of error, \(\epsilon\) is \(0.5\), we have that,

\[ n \geq \left(\dfrac{z_{1-(\alpha/2)}\; s}{\epsilon}\right)^2 = \left(\dfrac{2.5758293064439264 \times 4.7536998}{0.5}\right)^2 = 599.73 \approx 600 \]

Example Sample Size Determination for Proportion

To determine the sample size for estimating \(p=P(X>8)\) with 95% confidence to \(\pm 0.1\), we can use Equation @ref(eq:pSampleSize)

\[ n = \left(\dfrac{z_{1-(\alpha/2)}}{\epsilon}\right)^{2} \hat{p}(1 - \hat{p})=\left(\dfrac{z_{0.975}}{0.05}\right)^{2} (0.80)(1 - 0.80)=\left(\dfrac{1.96}{0.05}\right)^{2} (0.80)(0.20) = 245.86 \approx 246 \]

By using the maximum of \(\max{(246, 600)=600}\), we can re-run the simulation for this number of replications. Doing so, yields,

The half-width criteria are met.

Sample Size Determination via the KSL

fun main() {
    val rv = NormalRV(10.0, 16.0, streamNum = 3)
    val estimateX = Statistic("Estimated X")
    val estOfProb = Statistic("Pr(X>8)")
    val r = StatisticReporter(mutableListOf(estOfProb, estimateX))
    val n0 = 20 // sample size
    for (i in 1..n0) {
        val x = rv.value
        estimateX.collect(x)
        estOfProb.collect(x > 8)
    }
    println(r.halfWidthSummaryReport())
    val desiredHW = 0.5
    val s0 = estimateX.standardDeviation
    val level = 0.99
    val n = Statistic.estimateSampleSize(desiredHW = desiredHW,stdDev = s0,level = level)
    println("Sample Size Determination")
    println("desiredHW = $desiredHW")
    println("stdDev = $s0")
    println("Level = $level")
    println("recommended sample size = $n")
    val m = Statistic.estimateProportionSampleSize(desiredHW = 0.05,pEst = estOfProb.average,level = 0.95)
    println("Recommended sample size for proportion = $m")
}

Sample Size Determination via the KSL

The results match the values computed from the equations.

Sample Size Determination
desiredHW = 0.5
stdDev = 4.753707020199372
Level = 0.99
recommended sample size = 600
Recommended sample size for proportion = 246

Simulating the Game of Craps

Consider the game of “craps” as played in Las Vegas. The basic rules of the game are as follows: one player, the “shooter”, rolls a pair of dice. If the outcome of that roll is a 2, 3, or 12, the shooter immediately loses; if it is a 7 or an 11, the shooter wins. In all other cases, the number the shooter rolls on the first toss becomes the “point”, which the shooter must try to duplicate on subsequent rolls. If the shooter manages to roll the point before rolling a 7, the shooter wins; otherwise the shooter loses. It may take several rolls to determine whether the shooter wins or loses. After the first roll, only a 7 or the point have any significance until the win or loss is decided. Using the KSL random and statistic packages give answers and corresponding estimates to the following questions. Be sure to report your estimates in the form of confidence intervals.

1. Before the first roll of the dice, what is the probability the shooter will ultimately win?
2. What is the expected number of rolls required to decide the win or loss?

Game of Craps Solution

fun main() {
    val d1 = DUniformRV(1, 6, streamNum = 3)
    val d2 = DUniformRV(1, 6, streamNum = 4)
    val probOfWinning = Statistic("Prob of winning")
    val numTosses = Statistic("Number of Toss Statistics")
    val numGames = 5000
    for (k in 1..numGames) {
        var winner = false
        val point = d1.value.toInt() + d2.value.toInt()
        var numberoftoss = 1
        if (point == 7 || point == 11) {
            // automatic winner
            winner = true
        } else if (point == 2 || point == 3 || point == 12) {
            // automatic loser
            winner = false
        } else {// now must roll to get point...}

Game of Craps Solution

      else { // now must roll to get point
            var continueRolling = true
            while (continueRolling) {
                // increment number of tosses
                numberoftoss++
                // make next roll
                val nextRoll = d1.value.toInt() + d2.value.toInt()
                if (nextRoll == point) {// hit the point, stop rolling
                    winner = true
                    continueRolling = false
                } else if (nextRoll == 7) {// crapped out, stop rolling
                    winner = false
                    continueRolling = false
                }
            }
        }
        probOfWinning.collect(winner)
        numTosses.collect(numberoftoss.toDouble())
    }

Game of Craps Results

As we can see from the results of the statistics reporter, the probability of winning is just a little less than 50 percent.

Half-Width Statistical Summary Report - Confidence Level (95.000)% 

Name                            Count         Average      Half-Width 
-------------------------------------------------------------------------------- 
Prob of winning                   5000         0.4960          0.0139 
Number of Toss Statistics         5000         3.4342          0.0856 
-------------------------------------------------------------------------------- 

News Vendor Model

Sly’s convenience store sells BBQ wings for 25 cents a piece. They cost 15 cents a piece to make. The BBQ wings that are not sold on a given day are purchased by a local food pantry for 2 cents each. Assuming that Sly decides to make 30 wings a day, what is the expected revenue for the wings provided that the demand distribution is as show in in this table:

• What is the possible amount sold? If \(D\) is the demand for the day and we start with \(q\) items, then the most that can be sold is \(\min(D, q)\).
• What is the possible amount left over (available for salvage)? The amount left over at the end of the day is \(\max(0, q-D)\).
• Since the news vendor buys \(q\) items for \(c\), the cost of the items for the day is \(c \times q\)

News Vendor Model Solution

• The profit has the following form: \(G(D,q) = s \times \min(D, q) + u \times \max(0, q-D) - c \times q\)

fun main() {
    val q = 30.0 // order qty
    val s = 0.25 //sales price
    val c = 0.15 // unit cost
    val u = 0.02 //salvage value
    val values = doubleArrayOf(5.0, 10.0, 40.0, 45.0, 50.0, 55.0, 60.0)
    val cdf = doubleArrayOf(0.1, 0.3, 0.6, 0.8, 0.9, 0.95, 1.0)
    val dCDF = DEmpiricalRV(values, cdf, streamNum = 3)
    val stat = Statistic("Profit")
    val n = 100.0 // sample size
    var i = 1
    while (i <= n) {
        val d = dCDF.value
        val amtSold = minOf(d, q)
        val amtLeft = maxOf(0.0, q - d)
        val g = s * amtSold + u * amtLeft - c * q
        stat.collect(g)
        i++
    }
}

News Vendor Model Results

Running the model for 100 days results in the following output.

Count =      100.000000 
Average =    1.677500 
Std. Dev. =      2.201324 
Half-width =     0.436790 
95.0% CI = [1.2407095787617952, 2.1142904212382048]

Simple Inspection Process

A manufacturing process has a simple inspection process. The process manufactures items and places them in boxes of with four items per box. The manufacturing process has a defect rate of 15%. That is, the probability of an individual item being defective is 0.15. Assume that the chance of an item being defective is independent of any other item being defective. An inspector tests each box by randomly sampling one item from the box. If the selected item is good, the box is passed; if the selected item is defective, the box is rejected. We are interested in estimating the expected number of boxes that the inspector will look at before the first box is rejected.

Simple Inspection Process KSL Model

The probability that an item is defective can be modeled with a Bernoulli random variable with probability of success equal 0.15. In the following code, we model the inspection process with a function that counts the number of boxes that are inspected until the first rejected box.

fun main() {
    val itemRV = BernoulliRV(probOfSuccess = 0.15, streamNum = 3)
    val itemsPerBox = 4
    val stat = Statistic("Num Until Rejection")
    val sampleSize = 100
    for (i in 1..sampleSize) {
        val countBoxes = numUntilFirstRejection(itemRV, itemsPerBox)
        stat.collect(countBoxes)
    }
    print(String.format("%s \t %f %n", "Count = ", stat.count))
    print(String.format("%s \t %f %n", "Average = ", stat.average))
    print(String.format("%s \t %f %n", "Std. Dev. = ", stat.standardDeviation))
    print(String.format("%s \t %f %n", "Half-width = ", stat.halfWidth))
    println((stat.confidenceLevel * 100).toString() + "% CI = " + stat.confidenceInterval)
}

Simple Inspection Process KSL Model

/**
 *  This function counts the number of boxes inspected until the first
 *  box is found that has a randomly selected item that is defective.
 */
fun numUntilFirstRejection(itemRV: BernoulliRV, itemsPerBox: Int): Double {
    require(itemsPerBox >= 1) { "There must be at least 1 item per box" }
    var count = 0.0
    do {
        count++
        // randomly generate the box
        val box = itemRV.sample(itemsPerBox)
        // randomly sample from the box
        val inspection = KSLRandom.randomlySelect(box)
        // stops if bad = 1.0 is found, continues if good = 0.0 is found
    } while (inspection != 1.0)
    return count
}

Stochastic Activity Network

Consider a project consisting of nine jobs (A, B, \(\cdots\), I) with the precedence relationships and activity times. The table provides the minimum, mode, and maximum for each activity. We will assume that the activity times can be modeled using a triangular distribution with the appropriate parameters.

Stochastic Activity Network

Project Network Diagram

• Develop a simulation model that will estimate the expected time to completion for the project. In addition, estimate the probability that the project will be completed within 50 days. Finally, estimate the probability that the project will be completed between 42 and 48 days.

Stochastic Activity Network KSL Solution

• Data structures to represent the network.

    val activityRVs = mapOf<String, TriangularRV>(
        "A" to TriangularRV(2.0, 5.0, 8.0, streamNum = 2),
        "B" to TriangularRV(6.0, 9.0, 12.0, streamNum = 3),
        "C" to TriangularRV(6.0, 7.0, 8.0, streamNum = 4),
        "D" to TriangularRV(1.0, 4.0, 7.0, streamNum = 5),
        "E" to TriangularRV(7.0, 8.0, 9.0, streamNum = 6),
        "F" to TriangularRV(5.0, 14.0, 17.0, streamNum = 7),
        "G" to TriangularRV(3.0, 12.0, 21.0, streamNum = 8),
        "H" to TriangularRV(3.0, 6.0, 9.0, streamNum = 9),
        "I" to TriangularRV(5.0, 8.0, 11.0, streamNum = 10),
    )

    val paths = mutableListOf<List<String>>(
        listOf("A", "B", "D", "F", "H", "I"),
        listOf("A", "E", "F", "H", "I"),
        listOf("A", "C", "D", "F", "H", "I"),
        listOf("A", "B", "C", "G", "H", "I"),
    )

Stochastic Activity Network KSL Solution

• Function to generate activity times

fun generateActivityTimes(activities: Map<String, TriangularRV>) : Map<String, Double> {
    val map = mutableMapOf<String, Double>()
    for ((name, rv) in activities){
        map[name] = rv.value
    }
    return map
}

Stochastic Activity Network KSL Solution

• Function to compute the time to complete the network.

fun timeToCompletion(activityTimes: Map<String, Double>, paths: List<List<String>>): Double {
    val times = mutableListOf<Double>()
    for (path in paths) {
        var pathTime = 0.0
        for (activity in path) {
            pathTime = pathTime + activityTimes[activity]!!
        }
        times.add(pathTime)
    }
    return times.max()
}

Stochastic Activity Network KSL Solution

• Collecting the statistics on the generated networks.

    val timeToCompleteStat = Statistic("Time to Completion")
    val probLT50Days = Statistic("P(T<=50)")
    val probWith42and48Days = Statistic("P(42<=T<=48)")
    val interval = Interval(42, 48)
    val sampleSize = 1000
    for (i in 1..sampleSize) {
        val a = generateActivityTimes(activityRVs)
        val t = timeToCompletion(a, paths)
        timeToCompleteStat.collect(t)
        probLT50Days.collect(t <= 50)
        probWith42and48Days.collect(interval.contains(t))
    }
    val sr = StatisticReporter(mutableListOf(timeToCompleteStat, probLT50Days, probWith42and48Days))
    println(sr.halfWidthSummaryReport())

Stochastic Activity Network KSL Results

• The results indicate that on average it takes about 47.8 days to complete the network.

Statistical Summary Report

News Vendor Model as a Monte-Carlo Experiment

Because running Monte-Carlo experiments is very common, the KSL has provided some classes to support simple experimentation on Monte-Carlo problems within the ksl.utilities.mcintegration package.

Implement the MCReplicationIfc

To use the MCExperiment class the user can either subclass MCExperiment or provide a function that represents a replication of an individual observation. This can be done by using the MCReplicationIfc functional interface.

fun interface MCReplicationIfc {
    /**
     * @param j the current replication number. Could be used to implement more advanced
     * sampling that needs the current replication number. For example, antithetic sampling.
     */
    fun replication(j: Int): Double
}

News Vendor KSL Implementation

class NewsVendor(var demand: RVariableIfc) : MCReplicationIfc {
    var orderQty = 30.0 // order qty

    var salesPrice = 0.25 //sales price

    var unitCost = 0.15 // unit cost

    var salvageValue = 0.02 //salvage value

    override fun replication(j: Int): Double {
        val d = demand.value
        val amtSold = minOf(d, orderQty)
        val amtLeft = maxOf(0.0, orderQty - d)
        return salesPrice * amtSold + salvageValue * amtLeft - unitCost * orderQty
    }

}

News Vendor KSL Experiment

This code creates an instance of the news vendor problem and supplies the same demand distribution used in the previous section. The instance of NewsVendor, which implements the MCReplicationIfc interface is supplied to the MCExperiment as a constructor parameter. Then the desired half-width bound is set and the simulation executed.

fun main() {
    val values = doubleArrayOf(5.0, 10.0, 40.0, 45.0, 50.0, 55.0, 60.0)
    val cdf = doubleArrayOf(0.1, 0.3, 0.6, 0.8, 0.9, 0.95, 1.0)
    val dCDF = DEmpiricalRV(values, cdf, streamNum = 3)
    val nv = NewsVendor(dCDF)
    val exp = MCExperiment(nv)
    exp.desiredHWErrorBound = 0.01
    exp.runSimulation()
    println(exp)
}

News Vendor KSL Experiment Results

Monte Carlo Simulation Results
initial Sample Size = 30
max Sample Size = 10000
reset Stream OptionOn = false
Estimated sample size needed to meet criteria = 2045.0
desired half-width error bound = 0.01
actual half-width = 0.00999980532316225
error gap (hw - bound) = -1.946768377510122E-7
-----------------------------------------------------------
The half-width criteria was met!
-----------------------------------------------------------

News Vendor KSL Experiment Results

**** Sampling results ****
Number of macro replications executed = 2047.0
Number of micro replications per macro replication = 100
Total number of observations = 204700.0
ID 2
Name Statistic_1
Number 2047.0
Average 1.5049719101123593
Standard Deviation 0.23069885711869198
Standard Error 0.005099017725643279
Half-width 0.00999980532316225
Confidence Level 0.95
Confidence Interval [1.494972104789197, 1.5149717154355216]

MCExperiment Key Concepts

• The sampling is performed in the form of micro and macro replications. A micro replication observes the desired quantity for a pre-determined sample size. Then, the micro replications are repeated for a pre-determined overall number of macro replications.
• Let \(r\) be the number of micro replications and let \(n\) be the number of macro replications. Let \(Y_{ij}\) be the \(i^{th}\) observation of the \(r\) micro replications and let \(j\) be the \(j^{th}\) observation of the \(n\) macro replications.
• Let \(\bar{Y}_{\cdot j} = \frac{1}{r}\sum_{i=1}^{r}Y_{ij}\) be the sample average for the \(j^{th}\) macro replication. Let \(k\) be the sample size of the initial set of macro replications. Thus, they form a random sample of size \(k\), as follows:

\[ (\bar{Y}_{\cdot 1},\bar{Y}_{\cdot 2}, \dots, \bar{Y}_{\cdot k} ) \]

• The simulation is performed in two loops: an outer loop called the macro replications and an inner loop called the micro replications.

MCExperiment Key Concepts

• The user specifies a desired (half-width) error bound (\(\epsilon\)), an initial sample size (\(k\)), and a maximum sample size limit (\(M\)) for the macro replications.
  
• The initial sample size is used to generate a pilot sample from which an estimate of the number of samples needed to meet the half-width criteria. Let’s call the estimated sample size, \(m\). If \(m > k\), then an additional \((m-k)\) samples will be taken or until the error criteria is met or the maximum number of samples \(M\) is reached.
• Thus, if \(m > M\), and the error criterion is not met during the macro replications a total of \(M\) observations will be observed. Thus, the total number of macro replications will not exceed \(M\).
  
• If the error criteria is met before \(M\) is reached, the number of macro replications \(n\) will be somewhere between \(k\) and \(M\). For each of the \(n\), macro replications, a set of micro replications will be executed.
• Let \(r\) be the number of micro replications. The micro replications represent the evaluation of \(r\) observations of the Monte-Carlo evaluation. Thus, the total number of observations will be \(n \times r\).

Summary of Key Concepts

• Using the KSL’s random variate generation methods and the functionality for collecting and reporting statistics Monte Carlo experiments can be easily setup and executed.

• The quantities to be estimated from Monte Carlo experiments are random variables. We use the Monte Carlo experiment to generate a random sample from the experiment.

  • The random sample from the experiment provides estimates such as the sample average.
    • To estimate a probability, we can use an indicator variable.
  • We should report the sampling error or a confidence interval on the quantities estimated from the Monte Carlo experiment.

• Because we are performing an experiment, we should plan the sampling by understanding the sampling error and if possible, pre-plan the required sample size.

• Besides the random variate generation and statistical classes within the KSL, we can use the KSL’s MCExperiment class to execute Monte Carlo experiments.