Input Modeling Using the KSL

Learning Objectives

• To review the steps of distribution modeling
• To be able to estimate distribution parameters using the KSL
• To be able to fit continuous distributions using the KSL
• To be able to fit discrete distributions using the KSL

Input Distribution Modeling Process

1. Document the process being modeled. Describe the process being modeled and define the random variable to be collected.

2. Develop a plan for collecting the data and then collect the data.

3. Perform graphical and statistical analysis of the data. Using standard statistical analysis tools you should visually and statistically examine your data.

4. Hypothesize distributions: Using what you have learned from steps 1 - 3, you should hypothesize possible distributions for the data.

5. Estimate parameters. Once you have possible distributions in mind you need to estimate the parameters of those distributions.

6. Check the goodness of fit for hypothesized distributions. As part of this step you should perform some sensitivity analysis on your fitted model.

7. Recommend an appropriate distribution

Understand the Difference Between Continuous and Discrete Random Variables

• The key characteristic of a random variable is the set of values that it can have.
  • Variables that take on a countable set of values are discrete.
    
  • Variables that take on an uncountable set of values are said to be continuous.
• For discrete variables we can define a mapping from the set of integers to each possible value.
  • There may be an infinite number of values.
• Continuous variables are represented by the set (or subset) of real numbers.
  • That is, there are an uncountable infinite number of possible values that the variable can take on.
• Keep in mind that measuring devices (e.g. stop watch) are really discrete approximations for reality.
  • However, we need to understand the fundamental values associated with the variable not what our measuring devices limit us to.
• By understanding the modeling situations that produce data, you can hypothesize the appropriate distribution for the distribution fitting process.

Common Modeling Situations for Discrete Distributions

Common Modeling Situations for Continuous Distributions

Basic Statistics to Compute

• mean - Measures the central tendency of the distribution

• variance - Measures the spread about the mean

• coefficient of variation - \(c_v = \frac{\sqrt{Var[X]}}{E[X]}\)

• mode - Measures the most likely value

• skewness - Measures the asymmetry of the distribution about its mean. Helpful in picking an appropriate distribution.

• kurtosis - Measures the degree of peakedness of the distribution

• order statistics - Used when comparing the sample data to the theoretical distribution via P-P plots or Q-Q plots.

• quantiles (1st quartile, median, 3rd quartile) - Summarizes the distribution of the data.

• minimum, maximum, and range - Indicates the range of possible values

The KSL Statistic class and its companion object will estimate all of these quantities except for the mode.

Distribution Parameter Estimation

• There are two main methods for estimating the parameters of distribution functions.
  • method of moments (MOM): The method of moments matches the empirical moments to the theoretical moments of the distribution and attempts to solve the resulting system of equations.
  • method of maximum likelihood (MLE): The method of maximum likelihood attempts to find the parameter values that maximize the joint probability distribution function of the sample.
• The KSL facilitates the development of routines to estimate parameters by providing the ParameterEstimatorIfc interface.
  • Users can perform the estimation process using any algorithm that they deem appropriate.

KSL Classes for Parameter Estimation

The MVBSEstimatorIfc Interface

• MVBSEstimatorIfc - The purpose of this interface is to define a method that takes in data (in an array) and returns the estimated quantities in an array.

/**
 *  Given some data, produce multiple estimated statistics
 *  from the data and store the estimated quantities in
 *  the returned array. It is up to the user to interpret
 *  the array values appropriately.
 */
interface MVBSEstimatorIfc {
    /**
     * The name to associate with each dimension of the
     * array that is returned by estimate(). The names
     * should be unique. The order of the list of names should
     * match the order of elements in the returned array.
     */
    val names: List<String>

    fun estimate(data: DoubleArray): DoubleArray
}

The ParameterEstimatorIfc Interface

/**
 *  This interface defines a general function that uses 1-dimensional data, and estimates the parameters of some 
 * uni-variate distribution via some estimation algorithm. The basic contract is that the returned
 * EstimatedParameters is consistent with the required parameter estimation.
 */
interface ParameterEstimatorIfc {

    /**
     *  Indicates if the estimator requires that the range of the data be checked for a shift
     *  before the estimation process.
     */
    val checkRange: Boolean

    /**
     *  Estimates the parameters associated with some distribution. The returned [EstimationResult]
     *  needs to be consistent with the intent of the desired distribution.
     *  Note the meaning of the fields associated with [EstimationResult]
     */
    fun estimateParameters(data: DoubleArray, statistics: StatisticIfc = Statistic(data)): EstimationResult

}

Estimation Example

• General Steps
  • get the data (i.e. read it or generate it)
  • create the estimator
  • run the estimator on the data
  • review the results

fun main(){
    // define a normal random variable,
    val x = NormalRV(2.0, 5.0)
    // generate some data
    val data = x.sample(100)
    // create the estimator
    val ne = NormalMLEParameterEstimator
    // estimate the parameters
    val result = ne.estimateParameters(data)
    println(result)
}

Estimation Example Results

• estimateParameters(data: DoubleArray, statistics: StatisticIfc): EstimationResult
  • produces an EstimationResult instance
• The data class EstimationResult stores information about the estimation process to be returned as a result of its application on the supplied data.

Estimation Results:
The estimation was a SUCCESS!
Estimation message:
The normal parameters were estimated successfully using a MLE technique

Shift estimation results:
null

Estimated parameters:
RV Type = Normal
Double Parameters {mean=2.4138851463671918, variance=5.571943618235626}
Integer Parameters {}
Double Array Parameters {}

Class EstimationResult

class EstimationResult(
    val originalData: DoubleArray,
    var statistics: StatisticIfc,
    var shiftedData: ShiftedData? = null,
    val parameters: RVParameters? = null,
    var message: String? = null,
    var success: Boolean,
    val estimator: MVBSEstimatorIfc
)

• The success property indicates whether the parameter estimation process was successful.
  • If the success property is true, then the results should be meaningful.
• The statistics property is a statistical summary of the original data.
• The shiftedData property is the data if a shift was required.
• The parameters property is the estimated parameters as an instance of RVParameters.
• The message property is a text (string) representation of the distribution and its parameters.
• The estimator property can be used for bootstrapping results.

Continuous Distribution Recommendation Framework

Providing a Scoring Model

• To implement your own scoring model, you need to sub-class from the abstract base class PDFScoringModel.

• The primary method that needs to be implemented is the abstract score() function.

abstract fun score(data: DoubleArray, cdf: ContinuousDistributionIfc) : Score

• A Score is a data class that indicates whether the score is valid or not and provides the metric that determined the score.
  • You can think of the metric as the units of measure.
• The metric defines the domain (or set of legal values) for the metric and its direction.
  • By direction, we mean whether bigger scores are better than smaller scores or vice versa.
  • This information is used when a set of scores are combined into an overall score.

Default Scoring Models

• P-P Plot Sum of Squared Errors - Based on the concepts found in [@gan-koehler], the sum of squared error related to the P-P plot of the theoretical distribution versus the empirical distribution was developed as a scoring criterion. Let \((x_{(1)}, x_{(2)}, \ldots x_{(n)})\) be the order statistics.
• Let the theoretical distribution be represented with \(\hat{F}(x_{(i)})\) for i= 1, 2, \(\ldots\) n where \(\hat{F}\) is the CDF of the hypothesized distribution. Define the empirical distribution as

\[\tilde{F}_n(x_{(i)}) = \dfrac{i - 0.5}{n}\] Then, the P-P Plot sum of squared error criterion is defined as:

\[ \text{PP Squared Error} = \sum_{i = 1}^n (\tilde{F}_n(x_{(i)}) - \hat{F}(x_{(i)}))^2 \]

Default Scoring Models (continued)

• Anderson Darling criterion - The Anderson-Darling criterion is similar in spirit to the Cramer Von Mises test statistic except that it is designed to detect discrepancies in the tails of the distribution.

\[ A^2 = n\int_{-\infty}^{+\infty}\frac{\Big[ F_n(x) - F(x) \Big]^2}{F(x)(1-F(x))}\,dF(x) \]

• This metric can be computed from data sorted in ascending order (\(x_1, x_2, \cdots, x_n\)), i.e. the order statistics, as: \[ A^2 = -n - \sum_{i=1}^{n}\frac{2i-1}{n}\Big[ \ln(F(x_i)) + \ln(1-F(x_{n+1-i})) \Big] \]

Default Scoring Models (continued)

• Kolmogorov-Smirnov criterion - The Kolmogorov-Smirnov criterion, \(D_{n} = \max \lbrace D^{+}_{n}, D^{-}_{n} \rbrace\), represents the largest vertical distance between the hypothesized distribution and the empirical distribution over the range of the distribution.

\[ \begin{aligned} D^{+}_{n} & = \underset{1 \leq i \leq n}{\max} \Bigl\lbrace \tilde{F}_{n} \left( x_{(i)} \right) - \hat{F}(x_{(i)}) \Bigr\rbrace \\ & = \underset{1 \leq i \leq n}{\max} \Bigl\lbrace \frac{i}{n} - \hat{F}(x_{(i)}) \Bigr\rbrace \end{aligned} \]

\[ \begin{aligned} D^{-}_{n} & = \underset{1 \leq i \leq n}{\max} \Bigl\lbrace \hat{F}(x_{(i)}) - \tilde{F}_{n} \left( x_{(i-1)} \right) \Bigr\rbrace \\ & = \underset{1 \leq i \leq n}{\max} \Bigl\lbrace \hat{F}(x_{(i)}) - \frac{i-1}{n} \Bigr\rbrace \end{aligned} \]

Distribution Recommendation Framework

• The KSL allows the scoring results to be combined into an overall score.
• Suppose distribution \(F_k\) is evaluated by the five scoring models, each resulting in score \(S_i\) for \(i=1,2,\cdots, m\), where \(m\) is the number of scoring models.
• In general, the scores may not have the same scales and the same direction of goodness. The KSL scoring system translates and scales each score \(S_i\) in to a value, \(M_i\) measure, that has a common scale and direction (a bigger score is better).
  
• These value measures are then combined into an overall value (\(V_k\)) for the distribution using weights (\(w_i\)) across the scoring criteria:

\[ V_{k} = \sum_{i=1}^{m}w_i \times M_i \]

• The distribution that has the overall largest value, \(V_k\), is then recommended as the best fitting probability model.

PDFModeler Class Functionality

• Provides functionality to estimate the shift parameter, shift the data to the left, compute confidence intervals for the minimum and maximum of the data via bootstrapping, estimate the range of possible values for the distribution, create distributions from estimated parameters, create the default evaluation model, and work with histograms.

Using the PDFModeler Class

• The PDFModeling class encapsulates the estimation and scoring process.
• The purpose of this class is to serve as the general location for implementing the estimation of distribution parameters across many distributions.
• The general use involves the following:

val d = PDFModeler(data)
val estimationResults: List<EstimationResult> = d.estimateParameters(d.allEstimators)
val scoringResults = d.scoringResults(estimationResults)
val model = d.evaluateScoringResults(scoringResults)
scoringResults.forEach( ::println)

• The scoring results will be updated with the evaluation information and will contain the evaluation scores. The scoring results can be sorted to find the recommended distribution based on the evaluation score. Alternatively, the single function can be used:

val d = PDFModeler(data)
val results  = d.estimateAndEvaluateScores()

• This function returns an instance of PDFModelingResults, which will have the results of the entire fitting process.

Performing Statistical Tests

• The KSL recommendation is not based on P-values or statistical tests; however, the modeler may want to perform statistical tests for the top performing distribution (or others).
• The ContinuousCDFGoodnessOfFit class facilitates the performing of the following goodness of fit tests: Chi-Squared, Kolmogorov-Smirnov, Anderson-Darling, Cramer-Von Mises.

A Note on Histogram Intervals

• The Chi-squared test is based on a histogram summary of the data.
• This involves dividing the range of observations via break points \(b_j\) and tabulating frequencies or proportions of the data falling with the defined intervals.
• Let \(c_{j}\) be the observed count of the \(x\) values contained in the \(j^{th}\) interval \(\left(b_{j-1}, b_{j} \right]\),
• Let \(h_j = c_j/n\) be the relative frequency for the \(j^{th}\) interval, and
• Let \(p_j\) be the theoretical probability associated with the interval.

\[ p_j = P\{b_{j-1} \leq X \leq b_{j}\} = \int\limits_{b_{j-1}}^{b_{j}} f(x) \mathrm{d}x \]

Summarizing the Continuous Distribution Modeling Framework

The KSL continuous distribution modeling framework allows the modeler to:

• Define and select from a set of distributions to evaluate
• Define and select from a set of parameter estimation methods
• Compute bootstrap estimates of the confidence intervals for the estimated parameters
• Define and select from a set of distribution fit quality metrics (scoring models)
• Combine scoring models into an overall measure
• Recommend the best distribution based on weighted preference among evaluation metrics
• Make observations plots, histograms, autocorrelation plots, P-P plots, Q-Q plots, and empirical distribution plots.

Example

One hundred observations of the service time were collected using a portable digital assistant and are shown in here where the first observation is in row 1 column 1, the second observation is in row 2 column 1, the \(21^{st}\) observation is in row 1 column 2, and so forth.

Three Diagnostic Plots

Distribution Fitting

val results  = d.estimateAndEvaluateScores()
println("PDF Estimation Results for each Distribution:")
println("------------------------------------------------------")
results.sortedScoringResults.forEach(::println)
val topResult = results.sortedScoringResults.first()
val scores = results.evaluationModel.alternativeScoresAsDataFrame("Distributions")
println(scores)
println()
val values = results.evaluationModel.alternativeValuesAsDataFrame("Distributions")
println(values)
println()
val distPlot = topResult.distributionFitPlot()
distPlot.showInBrowser("Recommended Distribution ${topResult.name}")
println()
println("** Recommended Distribution** ${topResult.name}")
println()

Distribution Fitting Results

Comments on Distribution Results

• By applying the value scoring model and sorting, we can see the top performing distribution with overall value of 0.9696257711895793.
• The recommended distribution is: 36.83628655364795 + Exponential(mean=145.94151344635202).
• Notice that the shift parameter was automatically estimated for this situation.

Goodness of Fit Analysis

We can use the following code to perform a goodness of fit analysis on the recommended distribution.

val gof = ContinuousCDFGoodnessOfFit(topResult.estimationResult.testData,
  topResult.distribution,
  numEstimatedParameters = topResult.numberOfParameters
)
println(gof)

Goodness of Fit Results

GOF Results for Distribution: Exponential(mean=145.94151344635202)
---------------------------------------------------------

Chi-Squared Test Results:
Bin Label                 P(Bin)       Observed   Expected
  1 [ 0.00, 7.49)         0.050000            3    5.00000
  2 [ 7.49,15.38)         0.050000            4    5.00000
  3 [15.38,23.72)         0.050000            8    5.00000
  4 [23.72,32.57)         0.050000            3    5.00000
  5 [32.57,41.98)         0.050000            8    5.00000
  6 [41.98,52.05)         0.050000            4    5.00000
  7 [52.05,62.87)         0.050000            5    5.00000
  8 [62.87,74.55)         0.050000            8    5.00000
  9 [74.55,87.25)         0.050000            4    5.00000
 10 [87.25,101.16)        0.050000            2    5.00000
 11 [101.16,116.54)       0.050000            9    5.00000
 12 [116.54,133.72)       0.050000            3    5.00000
 13 [133.72,153.21)       0.050000            4    5.00000
 14 [153.21,175.71)       0.050000            6    5.00000
 15 [175.71,202.32)       0.050000            3    5.00000
 16 [202.32,234.88)       0.050000            4    5.00000
 17 [234.88,276.87)       0.050000            7    5.00000
 18 [276.87,336.04)       0.050000            4    5.00000
 19 [336.04,437.20)       0.050000            5    5.00000
 20 [437.20,1596.00)      0.049982            6    4.99822   *** Warning: expected < 5 ***

Goodness of Fit Results

Number of estimated parameters = 1
Number of intervals = 20
Degrees of Freedom = 18
Chi-Squared Test Statistic = 16.000784446070348
P-value = 0.5924925929222188
Hypothesis test at 0.05 level: 
The p-value = 0.5924925929222188 is >= 0.05 : Do not reject hypothesis.

Goodness of Fit Test Results:
K-S test statistic = 0.0399213095618332
K-S test p-value = 0.9954253138552503

Anderson-Darling test statistic = 0.24859556182926212
Anderson-Darling test p-value = 0.9710833730000956

Cramer-Von-Mises test statistic = 0.025283543923658194
Cramer-Von-Mises test p-value = 0.9893094979248238

Review the Diagnostic Plots

Finally, we can review the distribution plots and see that according the the P-P plot, Q-Q plot, histogram overlay, and empirical cumulative distribution overlay plot that the exponential distribution is an excellent probability model for this data.

Making it Super Easy

• All of the previously illustrated analysis can be performed with a few lines of code:

val data = KSLFileUtil.scanToArray(myFile.toPath())
val d = PDFModeler(data)
d.showAllResultsInBrowser()

Modeling Discrete Distributions

• The KSL also has the functionality to fit discrete distributions. That functionality is discussed here.

• The focus is on parameter estimation and statistical testing, not recommending the distribution.

Example Fitting a Negative Binomial

val dist = NegativeBinomial(0.2, theNumSuccesses = 4.0)
val rv = dist.randomVariable
rv.advanceToNextSubStream()
val data = rv.sample(200)
val breakPoints = PMFModeler.makeZeroToInfinityBreakPoints(data.size, dist)
val pf = DiscretePMFGoodnessOfFit(data, dist, breakPoints = breakPoints)
println(pf.chiSquaredTestResults())

• The key challenge is automatically making appropriate bins. Similar to the continuous case the KSL attempts to make bins that have approximately the same expected counts.

Example Results

Chi-Squared Test Results:
Bin Label                 P(Bin)       Observed   Expected
  1 [ 0.00, 3.00)         0.016960            4    3.39200   *** Warning: expected < 5 ***
  2 [ 3.00, 5.00)         0.039322            7    7.86432
  3 [ 5.00, 6.00)         0.029360           10    5.87203
  4 [ 6.00, 7.00)         0.035232            7    7.04643
  5 [ 7.00, 8.00)         0.040265            7    8.05306
  6 [ 8.00, 9.00)         0.044292            5    8.85837
  7 [ 9.00,10.00)         0.047245           10    9.44893
  8 [10.00,11.00)         0.049134            9    9.82689
  9 [11.00,12.00)         0.050028            6    10.0056
 10 [12.00,13.00)         0.050028           11    10.0056
 11 [13.00,14.00)         0.049258            6    9.85162
 12 [14.00,15.00)         0.047851           12    9.57015
 13 [15.00,16.00)         0.045937            7    9.18734
 14 [16.00,17.00)         0.043640           17    8.72798
 15 [17.00,18.00)         0.041073            4    8.21457
 16 [18.00,19.00)         0.038335            7    7.66693
 17 [19.00,20.00)         0.035510            7    7.10200
 18 [20.00,21.00)         0.032669            9    6.53384
 19 [21.00,22.00)         0.029869            5    5.97379
 20 [22.00,23.00)         0.027154            3    5.43072
 21 [23.00,24.00)         0.024556            4    4.91126   *** Warning: expected < 5 ***
 22 [24.00,26.00)         0.041903           13    8.38058
 23 [26.00,28.00)         0.033376            9    6.67520
 24 [28.00,31.00)         0.036998            6    7.39969
 25 [31.00,36.00)         0.036799            7    7.35973
 26 [36.00,Infinity)      0.033207            8    6.64147

Example Results

Number of estimated parameters = 2
Number of intervals = 26
Degrees of Freedom = 23
Chi-Squared Test Statistic = 25.69659399548871
P-value = 0.31539706650371313
Hypothesis test at 0.05 level: 
The p-value = 0.31539706650371313 is >= 0.05 : Do not reject hypothesis.

Plots for Discrete Distributions

• Use the frequency plot

val f = IntegerFrequency(data)
val fp = f.frequencyPlot()
fp.showInBrowser()

• Use the PMF comparison plot

val pmfModeler = PMFModeler(data)
val results = pmfModeler.estimateParameters(setOf(PoissonMLEParameterEstimator))
val e = results.first()
println(e)
val mean = e.parameters!!.doubleParameter("mean")
val pf = PoissonGoodnessOfFit(data.toDoubles(), mean = mean)
println(pf)
val plot = PMFComparisonPlot(data, pf.distribution)
plot.saveToFile("Lab_Count_PMF_Plot")
plot.showInBrowser()

Summary

• The KSL provides a full range of capabilities to fit continuous distributions.
• The KSL will automatically recommend a continuous distribution based on a prescribed scoring model.
• The KSL provides the ability to perform classic statistical tests for continuous distributions.
• The KSL provides the ability to estimate the parameters of discrete distributions and to evaluate their goodness of fit.
• The KSL’s input modeling capabilities are on par with many commercial packages.