Collecting Statistics Using the KSL

Learning Objectives

• To be able to collect statistics using classes within the KSL
• To understand the basics of statistical computations supported by the KSL

Collecting Statistics

• The ksl.utilities.statistics package implements statistical collection.
  
• The methods of the CollectorIfc interface are illustrated in the figure on the next slide.
• collect methods observe and collect various quantities based on the presented values.
• The collect method has been overloaded to facilitate collection of double values, arrays of double values, and boolean values.
• A collector can be reset. Resetting a collector should set the state of the collector as if no values had been observed.

The CollectorIfc Interface

Major Classes and Interface of the Statistics Package

• There are two major kinds of statistics, one of which assumes that the values presented must be weighted, the WeightedStatisticIfc interface and the WeightedStatistic class.
• The other branch of classes, derived from AbstractStatistic do not necessarily have to be weighted.
• The main classes to be discussed are Statistic, and Histogram.
• There are classes to tabulate frequencies, IntegerFrequency, and state related frequencies, StateFrequency, including single-step transition probabilities.
• The BatchStatistic class provides for computing statistics across batches of observations.

Major Classes and Interface of the Statistics Package

Creating and Using a Statistic

• The first statistic directly collects the generated values. The second statistic is designed to collect \(P(X\geq 20.0)\) by observing whether the generated value meets this criteria as defined by the boolean expression x >= 20.0.

fun main() {
    // create a normal mean = 20.0, variance = 4.0 random variable
    val n = NormalRV(20.0, 4.0, streamNum = 3)
    // create a Statistic to observe the values
    val stat = Statistic("Normal Stats")
    val pGT20 = Statistic("P(X>=20")
    // generate 100 values
    for (i in 1..100) {
        // getValue() method returns generated values
        val x = n.value
        stat.collect(x)
        pGT20.collect(x >= 20.0)
    }
    println(stat)
    println(pGT20)
    val reporter = StatisticReporter(mutableListOf(stat, pGT20))
    println(reporter.halfWidthSummaryReport())
}

toString() Output for Statistic Class

ID 1
Name Normal Stats
Number 100.0
Average 20.370190128861807
Standard Deviation 2.111292233346322
Standard Error 0.2111292233346322
Half-width 0.4189261806189412
Confidence Level 0.95
Confidence Interval [19.951263948242865, 20.78911630948075]
Minimum 15.020744984423821
Maximum 25.33588436770212
Sum 2037.0190128861807
Variance 4.457554894588499
Weighted Average 20.370190128861797
Weighted Sum 2037.0190128861796
Sum of Weights 100.0
Weighted Sum of Squares 41935.76252316213
Deviation Sum of Squares 441.2979345642614
Last value collected 21.110736402119805
Last weighted collected 1.0
Kurtosis -0.534855387072145
Skewness 0.20030433873223502
Lag 1 Covariance -0.973414579833684
Lag 1 Correlation -0.22057990840016864
Von Neumann Lag 1 Test Statistic -2.2136062395518343
Number of missing observations 0.0
Lead-Digit Rule(1) -1

Using a Statistical Reporter

• Any of the methods or properties of the Statistic class can be used to request specific statistical quantities.
• The StatisticReporter class, uses a list of objects that implement the StatisticIfc interface and facilitates the writing and printing of various statistical summary reports.

    val reporter = StatisticReporter(mutableListOf(stat, pGT20))
    println(reporter.halfWidthSummaryReport())

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

Name                           Count          Average      Half-Width 
--------------------------------------------------------------------------- 
Normal Stats                    100           20.3702          0.4189 
P(X>=20                         100            0.5100          0.0997 
--------------------------------------------------------------------------- 

Statistic Class Companion Object Methods

• estimateSampleSize(desiredHW: Double, stdDev: Double, level: Double) : Long - returns the approximate sample size necessary to reach the desired half-width at the specified confidence level given the estimate of the sample standard deviation.
• median(data: DoubleArray) : Double - the statistical median of the supplied array
• percentile(data: DoubleArray, p: Double) : Double - the \(p^{th}\) percentile of the data
• quantile(data: DoubleArray, q: Double) : Double - the \(q^{th}\) percentile of the data based on definition 7 of quantiles as per the R definitions.

Histograms

• A histogram tabulates counts and frequencies of observed data over a set of contiguous intervals.
  
• Let \(b_{0}, b_{1}, \cdots, b_{k}\) be the breakpoints (end points) of the class intervals such that \(\left(b_{0}, b_{1} \right], \left(b_{1}, b_{2} \right], \cdots, \left(b_{k-1}, b_{k} \right]\) form \(k\) disjoint and adjacent intervals. The intervals do not have to be of equal width.
  
• Also, \(b_{0}\) can be equal to \(-\infty\) resulting in interval \(\left(-\infty, b_{1} \right]\) and \(b_{k}\) can be equal to \(+\infty\) resulting in interval \(\left(b_{k-1}, +\infty \right)\).
• Define \(\Delta b_j = b_{j} - b_{j-1}\) and if all the intervals have the same width (except perhaps for the end intervals), \(\Delta b = \Delta b_j\).
• The key parameters of a histogram are the break points associated with the bins.

Histogram Class

Histogram Companion Object

The Histogram class’s companion object provides a number of methods for creating break points and recommending break points based on sample data. The function, createBreakPoints has input parameters:

• The first bin lower limit (\(b_{0}\)): This is the starting point of the range over which the data will be tabulated.
• The number of bins (\(k\)))
• The width of the bins, (\(\Delta b\))

The function recommendBreakPoints computes a reasonable lower limit, number of bins, and bin width based on the statistics.

Using the Histogram Class

This example illustrates how to configure the breakpoints for a histogram and to use the histogram to observe observations from a random sample.

fun main() {
    val d = ExponentialRV(2.0, streamNum = 1)
    val data = d.sample(1000)
    var bp = Histogram.recommendBreakPoints(data)
    bp = Histogram.addPositiveInfinity(bp)
    val h = Histogram(breakPoints = bp)
    for (x in data) {
        h.collect(x)
    }
    println(h)
    val plot = h.histogramPlot()
    plot.showInBrowser("Exponentially Distributed Data")
}

Output from the Histogram Class

-------------------------------------
Number of bins = 25
First bin starts at = 0.0
Last bin ends at = Infinity
Under flow count = 0.0
Over flow count = 0.0
Total bin count = 1000.0
Total count = 1000.0
-------------------------------------
Bin Range        Count CumTot Frac  CumFrac
  1 [ 0.00, 0.75)  = 333 333.0 0.333000 0.333000 
  2 [ 0.75, 1.50)  = 198 531.0 0.198000 0.531000 
  3 [ 1.50, 2.25)  = 150 681.0 0.150000 0.681000 
  4 [ 2.25, 3.00)  = 112 793.0 0.112000 0.793000 
  5 [ 3.00, 3.75)  = 75 868.0 0.075000 0.868000 
  6 [ 3.75, 4.50)  = 48 916.0 0.048000 0.916000 
  7 [ 4.50, 5.25)  = 28 944.0 0.028000 0.944000 
  8 [ 5.25, 6.00)  = 19 963.0 0.019000 0.963000 
  9 [ 6.00, 6.75)  = 16 979.0 0.016000 0.979000 
 10 [ 6.75, 7.50)  = 1 980.0 0.001000 0.980000 
etc.

Using the CachedHistogram Class

• The KSL also provides the CachedHistogram class, which stores the observed data within a data cache (array) until sufficient data has been observed to provide reasonable breakpoints.
  • The default size of the cache is 512 observations.

  fun main() {
    val d = ExponentialRV(2.0, streamNum = 1)
    val data = d.sample(1000)
    val ch = CachedHistogram()
    for (x in data) {
        ch.collect(x)
    }
    println(ch)
    val plot = ch.histogramPlot()
    plot.showInBrowser("Exponentially Distributed Data")
  }

Frequency Tabulation

• The KSL will also tabulate count frequencies when the values are only integers. This is accomplished with the IntegerFrequency class

Frequency Tabulation Example

In this example, an instance of the IntegerFrequency class is used to tabulate observations from a binomial random variable. A bar chart representation of the frequency is also shown within a browser window.

fun main() {
    val f = IntegerFrequency(name = "Frequency Demo")
    val bn = BinomialRV(0.5, 100, streamNum = 2)
    val sample = bn.sample(10000)
    f.collect(sample)
    println(f)
    val plot = f.frequencyPlot()
    plot.showInBrowser("Frequency Demo Plot")
}

Frequency Tabulation Example Output

Frequency Tabulation Frequency Demo
----------------------------------------
Number of cells = 37
Lower limit = -2147483648
Upper limit = 2147483647
Under flow count = 0
Over flow count = 0
Total count = 10000.0
Minimum value observed = 30
Maximum value observed = 69
----------------------------------------
Value    Count   Proportion
30   1.0     1.0E-4
33   3.0     3.0E-4
34   7.0     7.0E-4
35   10.0    0.001
36   17.0    0.0017
37   23.0    0.0023
38   48.0    0.0048
39   86.0    0.0086
40   101.0       0.0101
etc.
----------------------------------------

State Frequency Tabulation

The KSL provides the ability to define labeled states and to tabulate frequencies and proportions related to the visitation and transition between the states. This functionality is available in the StateFrequency class.

fun main() {
    val sf = StateFrequency(6)
    val states = sf.states
    for (i in 1..10000) {
        val state = KSLRandom.randomlySelect(states)
        sf.collect(state)
    }
    println(sf)
    val plot = sf.frequencyPlot(proportions = true)
    plot.showInBrowser("State Frequency Demo Plot")
}

State Frequency Output

State Frequency Tabulation for: Identity#1
State Labels
State{id=1, number=0, name='State:0'}
State{id=2, number=1, name='State:1'}
State{id=3, number=2, name='State:2'}
State{id=4, number=3, name='State:3'}
State{id=5, number=4, name='State:4'}
State{id=6, number=5, name='State:5'}
State transition counts
[288, 272, 264, 282, 265, 286]
[283, 278, 283, 286, 296, 266]
[286, 298, 263, 264, 247, 282]
[271, 263, 275, 279, 280, 294]
[274, 305, 273, 281, 296, 268]
[254, 277, 282, 270, 313, 255]

State Frequency Output

State transition proportions

[0.17380808690404345,0.16415208207604104,0.15932407966203982,0.17018708509354255,0.15992757996378998,0.17260108630054316]
[0.16725768321513002,0.16430260047281323,0.16725768321513002,0.1690307328605201,0.17494089834515367,0.15721040189125296]
[0.174390243902439,0.18170731707317073,0.1603658536585366,0.16097560975609757,0.15060975609756097,0.1719512195121951]
[0.16305655836341756,0.1582430806257521,0.1654632972322503,0.16787003610108303,0.1684717208182912,0.17689530685920576]
[0.1614614024749558,0.17972893341190335,0.16087212728344136,0.16558632881555688,0.17442545668827342,0.15792575132586917]
[0.15384615384615385,0.16777710478497881,0.17080557238037553,0.16353725015142337,0.18958207147183526,0.15445184736523318]

State Frequency Output

Frequency Tabulation Identity#1
----------------------------------------
Number of cells = 6
Lower limit = 0
Upper limit = 5
Under flow count = 0
Over flow count = 0
Total count = 10000
----------------------------------------
Value    Count   Proportion
0    1657    0.1657
1    1693    0.1693
2    1640    0.164
3    1662    0.1662
4    1697    0.1697
5    1651    0.1651
----------------------------------------

Batch Means Statistics

The method of batch means is based on observations \((X_{1}, X_{2}, X_{3}, \dots, X_{n})\). The idea is to group the output into batches of size, \(b\), such that the averages of the data within a batch are more nearly independent and possibly normally distributed.

\[\begin{multline*} \underbrace{X_1, X_2, \ldots, X_b}_{batch 1} \cdots \underbrace{X_{b+1}, X_{b+2}, \ldots, X_{2b}}_{batch 2} \cdots \\ \underbrace{X_{(j-1)b+1}, X_{(j-1)b+2}, \ldots, X_{jb}}_{batch j} \cdots \underbrace{X_{(k-1)b+1}, X_{(k-1)b+2}, \ldots, X_{kb}}_{batch k} \end{multline*}\]

Let \(k\) be the number of batches each of size \(b\), where, \(b = \lfloor \frac{n}{k}\rfloor\). Define the \(j^{th}\) batch mean (average) as:

\[ \bar{X}_j(b) = \dfrac{1}{b} \sum_{i=1}^b X_{(j-1)b+i} \]

Confidence Intervals on Batch Means

To form a \((1 - \alpha)\)% confidence interval, we simply treat this new series like a random sample and compute approximate confidence intervals using the sample average and sample variance of the batch means series:

\[ \bar{Y}(k) = \dfrac{1}{k} \sum_{j=1}^k Y_j \] The sample variance of the batch process is based on the \(k\) batches: \[ S_b^2 (k) = \dfrac{1}{k - 1} \sum_{j=1}^k (Y_j - \bar{Y}^2) \]

Confidence Intervals on Batch Means

Finally, if the batch process can be considered independent and identically distributed the \(1-\alpha\) level confidence interval can be written as follows: \[ \bar{Y}(k) \pm t_{\alpha/2, k-1} \dfrac{S_b (k)}{\sqrt{k}} \] This represents an interval estimate for the batch means estimator.

The BatchStatistic Class

The BatchStatistic class works with data as it is presented to its collect method. Three properties of the BatchStatistic class that are important are:

• minNumBatches – This represents the minimum number of batches required. The default value for this attribute is determined by BatchStatistic. MIN_NUM_BATCHES, which is set to 20.
• minBatchSize – This represents the minimum size for forming initial batches. The default value for this attribute is determined by BatchStatistic. MIN_NUM_OBS_PER_BATCH, which is set to 16.
• maxNumBatchesMultiple – This represents a multiple of minimum number of batches which is used to determine the upper limit (maximum) number of batches. For example, if maxNumBatchesMultiple = 2 and the minNumBatches = 20, then the maximum number of batches we can have is 40 (2*20).

Example BatchStatistic

fun main() {
    val d = ExponentialRV(2.0)
    // number of observations
    val n = 1000
    // minimum number of batches permitted
    // there will not be less than this number of batches
    val minNumBatches = 40
    // minimum batch size permitted
    // the batch size can be no smaller than this amount
    val minBatchSize = 25
    // maximum number of batch multiple
    //  The multiple of the minimum number of batches
    //  that determines the maximum number of batches
    val maxNBMultiple = 2
    // In this example, since 40*25 = 1000, the batch multiple does not matter
    val bm = BatchStatistic(minNumBatches, minBatchSize, maxNBMultiple)
    for (i in 1..n) {
        bm.collect(d.value)
    }

Example BatchStatistic

    val bma = bm.batchMeans
    var i = 0
    for (x in bma) {
        println("bm($i) = $x")
        i++
    }

    // this re-batches the 40 down to 10
    val reformed = bm.reformBatches(10)
    println(Statistic(reformed))
}