3.8 Monte-Carlo Experiments

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. Figure 3.7 illustrates the classes within the package.

Classes in ksl.utilities.mcintegration

Figure 3.7: Classes in ksl.utilities.mcintegration

A Monte-Carlo experiment is conceptualized as a two-step process: 1) an initial pilot sample, and 2) the sampling to meet the desired half-width criteria. The initial pilot sample determines based on the equations of Section 3.3.2 how many samples are needed to meet the criteria. Then, the remaining samples are executed until a maximum number of samples is reached or the half-width criteria is met.

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’s use some notation to represent these concepts.

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:

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

From this sample, the number of macro replications needed to meet the criteria is determined using Equation (3.8). The reason that this sample is formed from the \(r\) micro replications is to help ensure that the observations \(\bar{Y}_{\cdot j}\) are approximately normally distributed. According to the Central Limit Theorem, sample averages tend to be normally distributed as the number of observations within the sample increase towards infinity.

The simulation is performed in two loops: an outer loop called the macro replications and an inner loop called the micro replications. 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\).

By default, the number of macro replications should be relatively small and the number of micro replications large. Specific settings will be problem dependent. The default initial sample size, \(k\) is 30, with a maximum number of macro replications of \(M = 10000\). The default half-width error bound is \(0.001\). The default setting of the number of micro replications, \(r\), is 100. Again, these are all adjustable by the user via the properties of the MCExperiment class. The initial pilot sample is executed automatically as part of the overall sampling; however, using the runInitialSample() method, the user can just execute the initial sample and review the results before setting up the overall sampling. In the following code, we re-implement the news vendor problem using the MCExperiment class.

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
}

A Kotlin functional interface allows the function to be supplied easily as a parameter of a function. An interface with only one abstract method is called a functional interface, or a Single Abstract Method (SAM) interface. Using the functional programming constructs of Kotlin the user can supply a lambda expression as the function, implement the interface, or reference a class or object method that has the same functional signature. In the example presented here, for simplicity and clarity, we implement the interface.

class NewsVendor(var demand: RVariableIfc) : MCReplicationIfc {
    var orderQty = 30.0 // order qty
        set(value) {
            require(value > 0)
            field = value
        }
    var salesPrice = 0.25 //sales price
        set(value) {
            require(value > 0)
            field = value
        }
    var unitCost = 0.15 // unit cost
        set(value) {
            require(value > 0)
            field = value
        }
    var salvageValue = 0.02 //salvage value
        set(value) {
            require(value > 0)
            field = 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
    }

}

Notice the function, replication which implements computation of the profit equation for the news vendor problem. The key input parameters of the problem have also been implemented as properties of the class so that different problem instances can be easily created. The user of the class must provide a random variable to represent the demand. The running of the experiment is implemented in the following code.

Example 3.13 (News Vendor Via MCExperiment) This example illustrates how to use the MCExperiment class to run a Monte-Carlo experiment for the news vendor problem.

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)
    val nv = NewsVendor(dCDF)
    val exp = MCExperiment(nv)
    exp.desiredHWErrorBound = 0.01
    exp.runSimulation()
    println(exp)
}

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.

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!
-----------------------------------------------------------
**** 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]

As we can see from the results, the desired half-width criteria of 0.01 was met based on 2047 macro replications.

The KSL also provides a class called MC1DIntegration that is a subclass of MCExperiment which facilitates the evaluation of 1-dimensional integrals as previously illustrated in Section 3.2. For that situation a function and a sampling distribution can be provided.

  • Let \(f(x)\) be the probability distribution for the random variable supplied by the sampler.
  • Let \(g(x)\) be the function that needs to be integrated.
  • Let \(h(x)\) be a factorization of \(g(x)\) such that \(g(x) = h(x)*f(x)\), that is \(h(x) = g(x)/f(x)\)

The following code illustrates the use of the MC1DIntegration class to evaluate the following integral:

\[\theta = \int\limits_{0}^{\pi} \sin (x) \mathrm{d}x\] With \(g(x) = \sin(x)\), we need to specify a sampling distribution, \(f(x)\) over the range of \([0, \pi]\). By choosing, \(f(x)\) as the uniform distribution over \([0, \pi]\), we have that

\[f_{X}(x) = \begin{cases} \frac{1}{\pi} & 0 \leq x \leq \pi\\ 0 & \text{otherwise} \end{cases}\]

Thus,

\[h(x) = \begin{cases} \pi \sin(x) & 0 \leq x \leq \pi\\ 0 & \text{otherwise} \end{cases}\]

In the following code, notice that we multiply by \(\pi\) in the functional evaluation to account for the range of sampling.

Example 3.14 (1-D Integration via MCExperiment) This example illustrates how to use the MCExperiment class to run a Monte-Carlo experiment on a one-dimensional integration problem.

fun main() {
    class SinFunc : FunctionIfc {
        override fun f(x: Double): Double {
            return Math.PI * sin(x)
        }
    }

    val f = SinFunc()
    val mc = MC1DIntegration(f, UniformRV(0.0, Math.PI))
    println()
    mc.runSimulation()
    println(mc)
}

The results based on the default settings do not meet the desired half-width criteria. We leave it as an exercise for the reader to further explore this problem.

Monte Carlo Simulation Results
initial Sample Size = 30
max Sample Size = 10000
reset Stream OptionOn = false
Estimated sample size needed to meet criteria = 63129.0
desired half-width error bound = 0.001
actual half-width = 0.0025130156865938546
error gap (hw - bound) = 0.0015130156865938546
-----------------------------------------------------------
The half-width criteria was not met!
The user should consider one of the following:
1. increase the desired error bound using desiredHWErrorBound
2. increase the number of macro replications using maxSampleSize
2. increase the number of micro replications using microRepSampleSize
-----------------------------------------------------------
**** Sampling results ****
Number of macro replications executed = 10000.0
Number of micro replications per macro replication = 100
Total number of observations = 1000000.0
ID 2
Name Statistic_1
Number 10000.0
Average 1.9998788102129432
Standard Deviation 0.09754281939519215
Standard Error 9.754281939519214E-4
Half-width 0.0025130156865938546
Confidence Level 0.99
Confidence Interval [1.9973657945263494, 2.002391825899537]