Generating Random Variates Using the KSL

Learning Objectives

• To be able to generate random variates using the KSL
• To understand how to use the KSL for basic probability computations

Basic Random Variate Generation Methods

Appendix A describes four basic methods for generating random variates:

• inverse transform
• convolution
• acceptance/rejection
• mixture distributions

Exponential Random Variate Via Inverse Transform

Consider a random variable, \(X\), that represents the time until failure for a machine tool. Suppose \(X\) is exponentially distributed with an expected value of \(1.\overline{33}\). Generate a random variate for the time until the first failure using a uniformly distributed value of \(u = 0.7\).

fun main() {
    val u = 0.7
    val mean = 1.333333
    val x = rExpo(mean, u)
    println("Generated X = $x")
}

fun rExpo(mean: Double, u: Double) : Double {
    return -mean*ln(1.0 - u)
}

Negative Binomial Via Convolution

Use the following pseudo-random numbers \(u_{1} = 0.35\), \(u_{2} = 0.64\), \(u_{3} = 0.14\), generate a random variate from a shifted negative binomial distribution having parameters \(r=3\) and \(p= 0.3\).

fun main() {
    val u1 = 0.35
    val u2 = 0.64
    val u3 = 0.14
    val p = 0.3
    val x1 = rGeom(p, u1)
    val x2 = rGeom(p, u2)
    val x3 = rGeom(p, u3)
    val x = x1 + x2 + x3
    println("Generated X = $x")
}

fun rGeom(p: Double, u: Double): Double {
    val n = ln(1.0 - u)
    val d = ln(1.0 - p)
    return 1.0 + floor(n / d)
}

Random Variate Generation

• The KSL has the capability to generate random variates from both discrete and continuous distributions. The ksl.utilities.random.rvariable package supports this functionality.
• If d is a reference to an instance of a sub-class of type RVariable, then d.value generates a random value via the GetValueIfc interface.
• The RVariableIfc interface implements the RandomIfc interface, which provides stream control of the underlying random number stream.
• The SampleIfc interface provides a set of methods for sampling 1 or more values.

Random Variable Interfaces

The KSL can be used to model many different types of discrete random variables.

• BernoulliRV(probOfSuccess)
• BinomialRV(pSucces, numTrials)
• ConstantRV(constVal) a degenerate probability mass on a single value that cannot be changed
• DEmpiricalRV(values, cdf) values is an array of numbers and cdf is an array representing the cumulative distribution function over the values
• DUniformRV(min, max)
• GeometricRV(probOfSucces) with range is 0 to infinity
• NegativeBinomialRV(probOfSuccess, numSuccess) with range is 0 to infinity. The number of failures before the \(r^{th}\) success.
• PoissonRV(mean)
• ShiftedGeometricRV(probOfSucces) range is 1 to infinity
• VConstantRV(constVal) a degenerate probability mass on a single value that can be changed

The KSL can be used to model many different types of continuous random variables.

• BetaRV(alpha1, alpha2)
• ChiSquaredRV(degreesOfFreedom)
• ExponentialRV(mean)
• GammaRV(shape, scale)
• GeneralizedBetaRV(alpha1, alpha2, min, max)
• JohnsonBRV(alpha1, alpha2, min, max)
• LaplaceRV(mean scale)
• LogLogisticRV(shape, scale)

The KSL can be used to model many different types of continuous random variables.

• LognormalRV(mean, variance)
• NormalRV(mean, variance)
• PearsonType5RV(shape, scale)
• PearsonType6RV(alpha1, alpha2, beta)
• StudentTRV(degreesOfFreedom)
• TriangularRV(min, mode, max)
• UniformRV(min, max)
• WeibullRV(shape, scale)

To implement your own random variable classes

• Extend the RVarable class, which implements the RVariableIfc interface.
• The key method to implement is the protected generate() method, which should return the generated random value according to some algorithm.
• The KSL generally use the inverse transform method to generate random variates.
• Templates are available for generating via the acceptance rejection method via the AcceptanceRejectionRV and ratio of uniforms RatioOfUniformsRV classes.
• Random variables that are governed by a mixture of distributions can be formulated by using the MixtureRv class, which requires the random variables being mixed and the mixing distribution.
• Shifted random variables can be generated with the help of the ShiftedRV class.
• Generation from truncated distributions is facilitated via the TruncatedRV class.

Creating and using random variables

• Make the instance of the random variable
• Call the value property or value() function to get new values of the variate

fun main() {
    // create a normal mean = 20.0, variance = 4.0 random variable, 
    // using stream number 1
    val n = NormalRV(20.0, 4.0, streamNum = 1)
    print(String.format("%3s %15s %n", "n", "Values"))
    // generate some values
    for (i in 1..5) {
        // the value property returns a generated value
        val x = n.value
        print(String.format("%3d %15f %n", i, x))
    }
}

  n          Values 
  1       17.718732 
  2       19.056360 
  3       19.003682 
  4       21.875759 
  5       18.466600 

Using the SampleIfc Interface

• sample() returns a new value
• sample(sampleSize: Int) returns an array of values
• sampleInto(values: DoubleArray) fills an array with values
• sampleAsRows(sampleSize: Int, nRows: Int = 1) makes a 2-D array with samples in the rows
• sampleAsColumns(sampleSize: Int, nRows: Int = 1) makes a 2-D array with samples in the columns
• sampleInto(matrix: Array<DoubleArray>) fills a 2-D array with samples

fun main() {
    // create a triangular random variable with min = 2.0, mode = 5.0, max = 10.0
    val t = TriangularRV(2.0, 5.0, 10.0, streamNum = 1)
    // sample 5 values
    val sample = t.sample(5)
    print(String.format("%3s %15s %n", "n", "Values"))
    for (i in sample.indices) {
        print(String.format("%3d %15f %n", i + 1, sample[i]))
    }
}

Construct Random Variables with Particular Streams

• All random variables implement the stream control interface.
• This allows resetting, advancing, antithetic, common random numbers, etc. as previously discussed.

fun main() {
    // create a normal mean = 20.0, variance = 4.0, with the stream
    val n = NormalRV(20.0, 4.0, streamNum = 3)
    print(String.format("%3s %15s %n", "n", "Values"))
    for (i in 1..5) {
        // value property returns generated values
        val x = n.value
        print(String.format("%3d %15f %n", i, x))
    }
}

Creating and Using a Discrete Empirical Random Variable

The values are provided via an array and the probabilities must be specified in the form of the cumulative probabilities, such that the last element is 1.0.

fun main() {
    // values is the set of possible values
    val values = doubleArrayOf(1.0, 2.0, 3.0, 4.0)
    // cdf is the cumulative distribution function over the values
    val cdf = doubleArrayOf(1.0 / 6.0, 3.0 / 6.0, 5.0 / 6.0, 1.0)
    //create a discrete empirical random variable
    val n1 = DEmpiricalRV(values, cdf, streamNum = 1)
    println(n1)
    print(String.format("%3s %15s %n", "n", "Values"))
    for (i in 1..5) {
        print(String.format("%3d %15f %n", i, n1.value))
    }
    println()
    println(n1.probabilityPoints)
}

Creating and Using a Mixture Type Random Variable

Suppose the time that it takes to pay with a credit card, \(X_{1}\), is exponentially distributed with a mean of \(1.5\) minutes and the time that it takes to pay with cash, \(X_{2}\), is exponentially distributed with a mean of \(1.1\) minutes. In addition, suppose that the chance that a person pays with credit is 70%. Then, the overall distribution representing the payment service time, \(X\), has an hyper-exponential distribution with parameters \(\omega_{1} = 0.7\), \(\omega_{2} = 0.3\), \(\lambda_{1} = 1/1.5\), and \(\lambda_{2} = 1/1.1\). Generate 5 random variates from this distribution.

fun main() {
    // rvs is the list of random variables for the mixture
    val rvs = listOf(ExponentialRV(1.5), ExponentialRV(1.1))
    // cdf is the cumulative distribution function over the random variables
    val cdf = doubleArrayOf(0.7, 1.0)
    //create a mixture random variable
    val he = MixtureRV(rvs, cdf, streamNum = 1)
    print(String.format("%3s %15s %n", "n", "Values"))
    for (i in 1..5) {
        print(String.format("%3d %15f %n", i, he.value))
    }
}

Creating and Using a Random Variable from a Truncated Distribution

Suppose \(X\) represents the distance between two cracks in highway. Suppose that \(X\) has an exponential distribution with a mean of 10 meters. Generate 5 random distance values restricted between 3 and 6 meters.

fun main() {
    val cdf = Exponential(mean = 10.0)
    val rv = TruncatedRV(cdf, Interval(0.0, Double.POSITIVE_INFINITY),
        Interval(3.0, 6.0), streamNum = 1)
    print(String.format("%3s %15s %n", "n", "Values"))
    for (i in 1..5) {
        print(String.format("%3d %15f %n", i, rv.value))
    }
}

Creating and Using a Shifted Random Variable

Suppose \(X\) represents the time to setup a machine for production. From past time studies, we know that it cannot take any less than 5.5 minutes to prepare for the setup and that the time after the 5.5 minutes is random with a Weibull distribution with shape parameter \(\alpha = 3\) and scale parameter \(\beta = 5\). Generate 5 random observations of the setup time.

fun main() {
    val w = WeibullRV(shape = 3.0, scale = 5.0)
    val rv = ShiftedRV(5.0, w, streamNum = 1)
    print(String.format("%3s %15s %n", "n", "Values"))
    for (i in 1..5) {
        print(String.format("%3d %15f %n", i, rv.value))
    }
}

Using the Acceptance-Rejection Method

Consider the following PDF over the range \(\left[-1,1\right]\). Generate 1000 random variates from this distribution. \[ f(x) = \begin{cases} \frac{3}{4}\left( 1 - x^2 \right) & -1 \leq x \leq 1\\ 0 & \text{otherwise}\\ \end{cases} \]

Using the Acceptance-Rejection Method (continued)

• We can use \(g(x) = 3/4\) as the majorizing function, which results in \(c=3/2\), and \(w(x)\)

\[ w(x) = \begin{cases} \frac{1}{2} & -1 \leq x \leq 1\\ 0 & \text{otherwise}\\ \end{cases} \]

• Notice that, \(w(x)\) results in random variables, \(W\), where \(W \sim U(-1,1)\).
  
• Thus, the proposal distribution is simply a uniform distribution over the range from -1 to 1.

Using the Acceptance-Rejection Method (continued)

fun main() {
    // proposal distribution
    val wx = Uniform(-1.0, 1.0)
    // majorizing constant, if g(x) is majorizing function, then g(x) = w(x)*c
    val c = 3.0 / 2.0
    val rv = AcceptanceRejectionRV(wx, c, fOfx, streamNum = 1)
    val h = Histogram.create(-1.0, 1.0, 100)
    for (i in 1..1000) {
        h.collect(rv.value)
    }
    println(h)
    val hp = h.histogramPlot()
    hp.showInBrowser()
    val dp = DensityPlot(h) { x -> fOfx.pdf(x) }
    dp.showInBrowser()
}

object fOfx : PDFIfc {
    override fun pdf(x: Double): Double {
        if ((x < -1.0) || (x > 1))
            return 0.0
        return (0.75 * (1.0 - x * x))
    }

    override fun domain(): Interval {
        return Interval(-1.0, 1.0)
    }

}

Using KSLRandom to Generate Random Variates

• Although not recommended the KSLRandom object has functions defined for many of the distributions.
• Each method is named with an "r" in front of the distribution name.
• By using an import of KSLRandom functions, users can more conveniently call these methods.

fun main() {
    val v = rUniform(10.0, 15.0) // generate a U(10, 15) value
    val x = rNormal(5.0, 2.0) // generate a Normal(mu=5.0, var= 2.0) value
    val n = rPoisson(4.0).toDouble() //generate from a Poisson(mu=4.0) value
    print(String.format("v = %f, x = %f, n = %f %n", v, x, n))
}

Generating from Arrays and Lists

• KSLRandom has functions for sampling from arrays and lists without replacement.
• Extension functions for arrays and lists have also been defined.

fun main() {
    // create a list
    val strings = listOf("A", "B", "C", "D")
    // randomly pick from the list, with equal probability
    for (i in 1..5) {
        println(KSLRandom.randomlySelect(strings))
    }
    println()
    for (i in 1..5) {
        println(strings.randomlySelect())
    }
}

A
B
B
D
A

Permuting and Sampling from Populations

fun main() {
    // create an array to hold a population of values
    val y = DoubleArray(10)
    for (i in 0..9) {
        y[i] = (i + 1).toDouble()
    }
    // create the population
    val p = DPopulation(y, streamNum = 1)
    println(p.contentToString())
    println("Print the permuted population")
    // permute the population
    p.permute()
    println(p.contentToString())

[1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
Print the permuted population
[8.0, 10.0, 2.0, 4.0, 3.0, 5.0, 1.0, 6.0, 9.0, 7.0]

Permuting and Sampling from Populations

    // directly permute the array using KSLRandom
    println("Permuting y")
    KSLRandom.permute(y)
    println(y.contentToString())
    // sample from the population
    val x = p.sample(5)
    println("Sampling 5 from the population")
    println(x.contentToString())

Permuting y
[2.0, 4.0, 5.0, 9.0, 6.0, 8.0, 3.0, 1.0, 7.0, 10.0]
Sampling 5 from the population
[2.0, 4.0, 4.0, 9.0, 5.0]

Permuting and Sampling from Populations

    // create a string list and permute it
    val strList: MutableList<String> = ArrayList()
    strList.add("a")
    strList.add("b")
    strList.add("c")
    strList.add("d")
    println("The mutable list")
    println(strList)
    KSLRandom.permute(strList)
    println("The permuted list")
    println(strList)
    println("Permute using extension function")
    strList.permute()
    println(strList)
}

The mutable list
[a, b, c, d]
The permuted list
[c, a, b, d]
Permute using extension function
[b, d, c, a]

Creating and Using Random Variable Functionals

• The KSL also contains an algebra for working with random variables.
• Let \(f(\cdot)\) be an arbitrary function and let \(X\) be a random variable. Then, the quantity \(Y = f(X)\) is also a random variable.
• We can defined two random variables and then a third that is the sum of the first two random variables.
• Example: Generate Erlang distributed random variables as a convolution of exponential random variables.

fun main(){
    var erlang: RVariableIfc = ExponentialRV(10.0, streamNum = 1)
    for(i in 1..4) {
        erlang = erlang + ExponentialRV(10.0)
    }
    println("stream number = ${erlang.streamNumber}")
    val sample = erlang.sample(1000)
    val stats = Statistic(sample)
    println()
    print(stats)
    sample.writeToFile("erlang.txt")
}

Creating and Using Random Variable Functionals

• The interface RVariableIfc and base class RVariable provides the ability to construct new random variables that are functions of other random variables by overriding the \((+, -, \times, \div)\) operators.
• New random variables can be defined as functions of other random variables.
• If \(X \sim N(\mu, \sigma^2)\) then the random variable \(Y=e^X\) will be lognormally distributed \(LN(\mu_l,\sigma_{l}^{2})\).

fun main(){
    // define a lognormal random variable, y, in terms of a normal random variable
    val x = NormalRV(2.0, 5.0, streamNum = 1)
    val y = exp(x)
    // generate from y
    val ySample = y.sample(1000)
    println(ySample.statistics())
}

Creating and Using Random Variable Functionals

• If \(Y_1 \sim Gamma(\alpha_1, 1)\) and \(Y_2 \sim Gamma(\alpha_2, 1)\), then \(X = Y_1/(Y_1 + Y_2)\) has a Beta(\(\alpha_1, \alpha_2\)) distribution.

fun main(){
    // define a beta random variable in terms of gamma
    val alpha1 = 2.0
    val alpha2 = 5.0
    val y1 = GammaRV(alpha1, 1.0, streamNum = 2)
    val y2 = GammaRV(alpha2, 1.0)
    val betaRV = y1/(y1+y2)
    val betaSample = betaRV.sample(500)
    println("stream number = ${betaRV.streamNumber}")
    println(betaSample.statistics())
}

Using Probability Models

• The ksl.utilities.distributions package facilitates computing from distributions.

• The DistributionIfc interface requires the following functions.

• cdf(b: Double) - computes the cumulative probability, \(F(b) = P(X \leq b)\)

• cdf(a: Double, b: Double) - computes the cumulative probability, \(P( a \leq X \leq b)\)

• complementaryCDF(b: Double) - computes the cumulative probability, \(1-F(b) = P(X > b)\)

• mean() - returns the expected value (mean) of the distribution

• variance() - returns the variance of the distribution

• standardDeviation() - returns the standard deviation of the distribution

• invCDF(p: Double) - returns the inverse of the cumulative distribution function \(F^{-1}(p)\). This is performed by numerical search if necessary

Key Points About Probability Models

• Discrete and continuous distributions are defined by the interfaces they implement.
• The discrete distributions such as the geometric, binomial, etc. implement the DiscreteDistributionIfc and PMFIfc interfaces.
• Instance of distributions are mutable. Their parameters can change. Instances of random variables are immutable. Their parameters cannot change.
• All distributions have functions to create random variables.

Illustrative Calculations

• Create an instance of the Binomial distribution and use it to compute some probability values.

fun main() {
    // make and use a Binomial(p, n) distribution
    val n = 10
    val p = 0.8
    println("n = $n")
    println("p = $p")
    val bnDF = Binomial(p, n)
    println("mean = " + bnDF.mean())
    println("variance = " + bnDF.variance())
    // compute some values
    print(String.format("%3s %15s %15s %n", "k", "p(k)", "cdf(k)"))
    for (i in 0..10) {
        print(String.format("%3d %15.10f %15.10f %n", i, bnDF.pmf(i), bnDF.cdf(i)))
    }
    println()

Illustrative Calculations Output

n = 10
p = 0.8
mean = 8.0
variance = 1.5999999999999996
  k            p(k)          cdf(k) 
  0    0.0000001024    0.0000001024 
  1    0.0000040960    0.0000041984 
  2    0.0000737280    0.0000779264 
  3    0.0007864320    0.0008643584 
  4    0.0055050240    0.0063693824 
  5    0.0264241152    0.0327934976 
  6    0.0880803840    0.1208738816 
  7    0.2013265920    0.3222004736 
  8    0.3019898880    0.6241903616 
  9    0.2684354560    0.8926258176 
 10    0.1073741824    1.0000000000 

Mutating the Distribution

• Instances of RVariable cannot be mutated. However, instances of distributions, for example Binomial can have their parameters mutated.
• Distributions can also create random variables.

    // change the probability and number of trials
    bnDF.probOfSuccess = 0.5
    bnDF.numTrials = 20
    println("mean = " + bnDF.mean())
    println("variance = " + bnDF.variance())
    // make random variables based on the distributions
    val brv = bnDF.randomVariable(1)
    print(String.format("%3s %15s %n", "n", "Values"))
    // generate some values
    for (i in 1..5) {
        // value property returns generated values
        val x = brv.value.toInt()
        print(String.format("%3d %15d %n", i, x))
    }
}

Some Useful Binomial Distribution Companion Object Functions

• binomialPMF(j: Int, n: Int, p: Double) - directly computes the probability for the value \(j\)
• binomialCDF(j: Int, n: Int, p: Double) - directly computes the cumulative distribution function for the value \(j\)
• binomialCCDF(j: Int, n: Int, p: Double)- directly computes the complementary cumulative distribution function for the value of \(j\)
• binomialInvCDF(x: Double, n: Int, p: Double) - directly computes the inverse cumulative distribution function

Some Useful Normal Distribution Companion Object Functions

• stdNormalCDF(z: Double) - the cumulative probability for a \(Z ~ N(0,1)\) random variable, i.e. \(F(z) = P(Z \leq z)\)
• stdNormalComplementaryCDF(z: Double) - returns \(1-P(Z \leq z)\)
• stdNormalInvCDF(p: Double) - returns \(z = F^{-1}(p)\) the inverse of the cumulative distribution function

Some Useful Student-T Distribution Companion Object Functions

• cdf(dof: Double, x: Double) - computes the cumulative distribution function for \(x\) given the degrees of freedom
• invCDF(dof: Double, p: Double) - computes the inverse cumulative distribution function or t-value for the supplied probability given the degrees of freedom.

Additional Functionality

• RArrays - defines extension functions for randomly sampling from arrays and some lists
• MVSampleIfc - defines the interface for generating random arrays of data
• MVRVariableIfc - the multi-variate analog for the RVariableIfc
• MVRVariable - an abstract base class for defining multi-variate random variables
• MVIndependentMarginals - a concrete implementation for generating independent vectors that have specified random variates for each coordinate of the vector.