Introductory Note

This is a write-up based on a research paper my group mates and I wrote for our Stat 142 (Introduction to Computational Statistics) course, under Engr. Xavier Bilon, when we were completing our Bachelor’s degree in Statistics at University of the Philippines - Diliman.

The paper, entitled “Simulated Normality (SiNo) Test”, was authored by:

Special thanks to my friend Elle Ordonio for helping me proofread this write-up.

Click on the hyperlinks above and connect with us on LinkedIn!

Description

The Simulated Normality (or SiNo) Test is a Lilliefors-like Kernel Density Estimate Test for Normality. It uses the integral of the L1-distance over the reals as the test statistic. The null distribution is estimated using Monte Carlo simulations. The realized sample is first standardized (z-score) and then, the kernel density estimate is obtained based on the standardized sample. This kernel density estimate is then compared with the standard normal density.

Test Formulation

Suppose that \(X_{1}\) and \(X_{2}\) are continuous random variables with probability density functions \(f_{X_{1}}\left(\cdot\right)\) and \(f_{X_{2}}\left(\cdot\right)\), respectively. Consider the following integral:

\[\int_{\mathbb{R}}{\left|f_{X_{1}}\left(x\right)-f_{X_{2}}\left(x\right) \right|}dx\]

If \(X_{1}\) and \(X_{2}\) were identically distributed, then

\[\int_{\mathbb{R}}{\left|f_{X_{1}}\left(x\right)-f_{X_{2}}\left(x\right) \right|}dx=0\]

because \(f_{X_{1}}\left(x\right)=f_{X_{2}}\left(x\right)\) for all \(x\in\mathbb{R}\). However, if \(X_{1}\) and \(X_{2}\) came from different distributions, then

\[\int_{\mathbb{R}}{\left|f_{X_{1}}\left(x\right)-f_{X_{2}}\left(x\right) \right|}dx\gt 0\]

From these, a natural interpretation of the integral arises. Larger integral values correspond to “more differences” between the two distributions (i.e., between the distributions of the random variables \(X_{1}\) and \(X_{2}\). To illustrate, consider the graph below.

Figure 1. N\left(0,1\right) (red) and N\left(2, 0.25\right) (blue) densities
Figure 1. \(N\left(0,1\right)\) (red) and \(N\left(2, 0.25\right)\) (blue) densities


The green region represents the difference between the two distributions. The area of this region is simply the corresponding form of the integral shown earlier. Thus, in this case, it would be

\[\int_{\mathbb{R}}{\left|\frac{1}{\sqrt{2\pi}}e^{-\frac{x^{2}}{2}} -\frac{1}{\frac{1}{2}\sqrt{2\pi}}e^{-\frac{1}{2}\left(\frac{x-2}{\frac{1}{2}} \right)^{2}}\right|}dx\gt 0\]

The absolute value is taken so that deviations in different directions do not “cancel out” and since we weigh deviations in either direction equally. Of course, another way to quantify the difference between two distributions is by utilizing the cumulative distribution functions instead of the density functions.

In many existing tests of normality, the empirical distribution function is used as the “observed” cumulative distribution function and compare this with the cumulative distribution function of the specified normal distribution (or the standard normal distribution). However, one can also utilize the availability of computational resources and proceed in a different direction.

Suppose \(\mathcal{F}\) is the parametric family of normal distributions. Moreover, let \(\left(X_{1},X_{2},\dots,X_{n}\right)\) be a random sample of size \(n\), where \(n\in\mathbb{Z}^{+}\), from a population with probability density function \(f_{X}(\cdot)\). Suppose \(f_{X}^{*}(\cdot)\) is the kernel density estimate, with certain parameters, for the true density \(f_{X}(\cdot)\) based on the realized random sample \(\left(x_{1},x_{2},\dots,x_{n}\right)\). The null and alternative hypotheses are:

\[H_{0}:f_{X}(\cdot)\in\mathcal{F}\] \[H_{1}:f_{X}(\cdot)\notin\mathcal{F}\]

Then, from Cao and Lugosi (2005), a test statistic \(T\) for testing the hypothesis that \(f_{X}(\cdot)\in\mathcal{F}\) is

\[\inf_{g\in\mathcal{F}}{\int_{\mathbb{R}}{\left|f_{X}^{*}(x)-g\right|}dx}\]

Here, we explore a parameter-free hypothesis test for normality inspired by Cao and Lugosi (2005) and by Lilliefors (1967). The null hypothesis is that \(f_{X}(\cdot)\in\mathcal{F}\).

From the original realized random sample \(\left(x_{1},x_{2},\dots,x_{n}\right)\), generate \(\left(y_{1},y_{2},\dots,y_{n}\right)\) such that for each \(i\in\{1,2,\dots,n\}\),

\[y_{i}=\frac{x_{i}-\overline{X}}{s_{X}}\]

, where \(\overline{X}\) and \(s_{X}\) are the sample mean and sample standard deviation of the original random sample, respectively. Now, suppose that \(f_{Y}^{*}(\cdot)\) is the kernel density estimate based on \(\left(y_{1},y_{2},\dots,y_{n}\right)\). Then, to test the hypothesis that \(f_{X}(\cdot)\in\mathcal{F}\), get the test statistic

\[T=\int_{\mathbb{R}}{\left|\phi(x)-f_{Y}^{*}(x)\right|}dx\]

, where \(\phi\left(\cdot\right)\) is the density function of the standard normal distribution. Instead of using the estimates as the parameters for the reference distribution similar to what Pavia (2015) outlined, the estimates are used to transform the sample before performing kernel density estimation and the reference distribution will always just be the standard normal density.

To aid in interpreting the test statistic, note that \(\overline{X}\) and \(s_{X}\) are both consistent and unbiased estimators of \(E(X)\) and \(V(X)\), respectively, where \(X\sim f_{X}(\cdot)\). Thus, as \(n\to+\infty\), if \(X\sim N\left(\mu,\sigma^{2}\right)\), then

\[\frac{X-\overline{X}}{s_{X}}\approx N(0,1)\]

, where \(u\in\mathbb{R}\) and \(\sigma^{2}>0\) are the parameters of the underlying normal distribution. If \(f_{X}(\cdot)\sim N\left(\mu,\sigma^{2} \right)\), it is expected that \(T\approx 0\).

Clearly, the minimum value of \(T\) is \(0\) and the set of possible values of \(T\) is a subset of \(\mathbb{R}^{+}\). For a given sample size, a higher test statistic implies greater evidence that the true population distribution from which the random sample was taken does not belong to the normal family of distributions.

Deriving the Null Distributions

The null distributions will be estimated using Monte Carlo simulations. However, the ancillarity of \(T\) under the null hypothesis must be established first.

Define \(K\left(\cdot ;z_{1},z_{2},\dots,z_{n}\right)\) as the kernel density estimate (with certain parameters) based on \(\left(z_{1},z_{2}, \dots,z_{n}\right)\). Then, re-expressing \(T\),

\[T\left(X_{1},X_{2},\dots,X_{n}\right)=\int_{\mathbb{R}}{\left|\phi(x)- K\left(x;\frac{X_{1}-\overline{X}}{s_{X}},\frac{X_{2}-\overline{X}}{s_{X}}, \dots,\frac{X_{n}-\overline{X}}{s_{X}}\right)\right|}dx\]

Let \(c_{1},c_{2}\in\mathbb{R}\). Note that \(\overline{X}\) is location-scale equivariant and that \(s_{X}\) is scale equivariant and location invariant. Then,

\[T\left(c_{1}X_{1}+c_{2},c_{1}X_{2}+c_{2},\dots,c_{1}X_{n}+c_{2}\right)\] \[=\int_{\mathbb{R}}{\left|\phi(x)-K\left(x;\frac{c_{1}X_{1}+c_{2}-c_{1} \overline{X}-c_{2}}{c_{1}s_{X}},\frac{c_{1}X_{2}+c_{2}-c_{1}\overline{X}- c_{2}}{c_{1}s_{X}},\dots,\frac{c_{1}X_{n}+c_{2}-c_{1}\overline{X}-c_{2}}{ c_{1}s_{X}}\right)\right|}dx\] \[=\int_{\mathbb{R}}{\left|\phi(x)-K\left(x;\frac{X_{1}-\overline{X}}{s_{X}} ,\frac{X_{2}-\overline{X}}{s_{X}},\dots,\frac{X_{n}-\overline{X}}{s_{X}}\right) \right|}dx=T\left(X_{1},X_{2},\dots,X_{n}\right)\]

Thus, \(T\left(c_{1}X_{1}+c_{2},c_{1}X_{2}+c_{2},\dots,c_{1}X_{n}+c_{2}\right) =T\left(X_{1},X_{2},\dots,X_{n}\right)\). In other words, \(T\) is location-scale invariant.

Now, assume that the null hypothesis is true. That is, \(f_{X}(\cdot)\in\mathcal{F}\). Note that \(F\), the family of normal distributions, is a location-scale family of densities. Therefore, since \(F\) is a location-scale family and \(T\) is location-scale invariant, it follows that \(T\) is an ancillary statistic.

Consequently, the distribution of \(T\) does not depend on the parameters of \(f_{X}(\cdot)\). That is, the distribution of \(T\) is invariant with respect to \(\mu\) and \(\sigma^{2}\) (the parameters that determine \(f_{X}(\cdot)\) or \(N\left(\mu,\sigma^{2}\right)\)).

Because \(T\) is ancillary, the distribution of \(T\) (under the null hypothesis) can be estimated by completing Monte Carlo simulations with the standard normal distribution as the population distribution (and the distribution of \(T\) will be the same, regardless of what “kind” of normal distribution the population distribution really is).

A total of \(M\in\mathbb{Z}^{+}\) random samples from \(N(0,1)\) will be taken. Suppose that \(T_{i}\) is the test statistic for the \(i\)th random sample. The estimated distribution of \(T\) can then be obtained from \(\{T_{1},T_{2}, \dots,T_{M}\}\). This is then the estimated null distribution.

The estimated null distributions for select sample sizes are shown below.

Figure 2. Estimated Null Distributions
Figure 2. Estimated Null Distributions


Implementation

R Script

The test is implemented in R. Four existing functions is used in the core mechanics of the program or function. Particularly, the following functions are used: stats::density, cubature::cubintegrate, stats::approxfun, and spatstat::CDF. To perform kernel density estimation, stats::density is used with mostly default parameters (although, some precision-related parameters are made variable in the function implementing the SiNo test). Some parameters of interest are bw and kernel, which are the bandwidth used for the estimation and the kernel smoothing function used for density estimation, respectively. At least as a first step, the default values are used, which are bw = bw.nrd0 and kernel = “gaussian”, which means that the Gaussian kernel is used and the bandwidth that will be selected will be based on the rule-of-thumb for bandwidth selection of the Gaussian kernel smoothing function.

Moreover, to approximate the integral in the test statistic, cubature::cubintegrate is used. Default values are used, with the exception of method, where instead of using the default method = “hcubature”, method = “pcubature” is used. This is because the number of dimensions is small (precisely 1) and it is reasonable to expect that most density estimates that will be encountered are smooth. By definition, the density estimate is similar to a smoothened version of the histogram and many common parametric continuous distributions (like the normal distribution - the one of most interest in this case) are smooth.

Lastly, to perform calculations like obtaining the p-value, the functions stats::approxfun and spatstat::CDF are used. This is because stats::density returns an object with class density. However, it is better to obtain a formula (or a function) for the estimated density function (afterall, the density is a function) that can be easily used for calculating probabilities. Thus, stats::approxfun is used on the output of stats::density to obtain an approximate function of the kernel density estimate (or simply said, to obtain the estimated density function). After that, spatstat::CDF is used to generate the cumulative distribution function, which are then used to calculate p-values (the cumulative distribution function is doubly useful here since SiNo Test is always right-tailed).

Ultimately, the current R Implementation of the SiNo Test is shown below.

sino_test <- function(x, simulations = 1000, precision = 1024, def_lim = 10,
                      ref_dist = NULL, sig_level = 0.05,
                      show_message = TRUE) {

  ## Function for calculating test statistic
  sino_stat <- function(dat, prec, dlim) {
    dat <- (dat - mean(dat)) / sd(dat)
    kde <- density(dat, n = prec, cut = dlim)
    kde_pdf <- approxfun(kde, yleft = 0, yright = 0)
    abs_diff <- function(v) {
      return(abs(kde_pdf(v) - dnorm(v)))
    }
    return(cubintegrate(abs_diff, lower = -Inf, upper = Inf,
                        method = "pcubature")$integral)
  }

  ## Function for estimating the null distribution
  sino_null <- function(size, prec, dlim, sims) {
    i <- 1
    test_stats <- rep(0, sims)
    while (i <= sims) {
      set.seed(i)
      rand_samp <- rnorm(size)
      test_stats[i] <- sino_stat(rand_samp, prec, dlim)
      i <- i + 1
    }
    kde <- density(test_stats, n = prec, cut = dlim)
    return(kde)
  }

  ## Perform the test
  if (is.null(ref_dist)) {
    if ((length(x) >= 5) && (length(x) <= 2000)) {
      null_dist <- SinoTest::null_densities[[length(x) - 4]]
    } else {
      null_dist <- sino_null(length(x), precision, def_lim, simulations)
    }
  } else {
    null_dist <- ref_dist
  }
  test_stat <- sino_stat(x, precision, def_lim)
  if (class(null_dist) == "density") {
    cdf_null <- CDF(null_dist)
  } else if (class(null_dist) == "function") {
    cdf_null <- null_dist
  }
  p_val <- 1 - cdf_null(test_stat)
  if (p_val <= sig_level) {
    interp <- "Reject the null hypothesis"
  } else {
    interp <- "Fail to reject the null hypothesis"
  }
  p_message <- ifelse(p_val == 0, "< 0.00001", p_val)
  if (show_message) {
    message(paste("\nSimulated Normality Test\n",
                  "Null Hypothesis: The population distribution is a ",
                  "normal distirbution.\n",
                  "Data: ", deparse(substitute(x)), "\n",
                  "Test Statistic: ", test_stat, "\n",
                  "p-value: ", p_message, "\n",
                  "Result: ", interp,
                  " at ", sig_level, " level of significance\n", sep = ""))
  }
  results <- list(data_name = deparse(substitute(x)), data = x,
                  statistic = test_stat, p_value = p_val,
                  simulations = simulations, precision = precision,
                  def_lim = def_lim,
                  ref_dist = ifelse(is.null(ref_dist), cdf_null,
                                    deparse(substitute(ref_dist))))

}

R Package and Installation

An R package, SinoTest, for the Simulated Normality Test was created, which includes the following:

  • sino_test - function implementing the test
  • sino_null - function for estimating the null distirbution
  • null_densities - estimates for the null distributions for \(n=5,6,7,\dots,2000\)

It is currently hosted on GitHub: https://github.com/jptuazon/SinoTest. One can install and load the package with the following script:

library(devtools)
install_github("jptuazon/SinoTest", dependencies = TRUE)
library(SinoTest)

The package includes documentation for the functions. If you encounter issues installing the package, this may help. Moreover, pre-installing (or loading before installing SinoTest) the packages stat, spatstat.explore, and cubature may resolve the issue.

Please note that this is still currently being developed. The scripts and documentations may still change and improve.

Examples

Below are some examples implementing the SiNo Test:

sino_test(rnorm(500))
## 
## Simulated Normality Test
## Null Hypothesis: The population distribution is a normal distirbution.
## Data: rnorm(500)
## Test Statistic: 0.0400743654581433
## p-value: 0.976213810495743
## Result: Fail to reject the null hypothesis at 0.05 level of significance
sino_test(rexp(500))
## 
## Simulated Normality Test
## Null Hypothesis: The population distribution is a normal distirbution.
## Data: rexp(500)
## Test Statistic: 0.432928927889442
## p-value: < 0.00001
## Result: Reject the null hypothesis at 0.05 level of significance
sino_test(runif(100), sig_level = 0.1)
## 
## Simulated Normality Test
## Null Hypothesis: The population distribution is a normal distirbution.
## Data: runif(100)
## Test Statistic: 0.289930511891535
## p-value: 0.000577103409977742
## Result: Reject the null hypothesis at 0.1 level of significance

Power Analysis

Here, the power of the SiNo test is estimated using 10000 Monte Carlo simulations for various distributions. Based on Wijekularathna, Manage, and Scariano (2019), several long-tailed symmetric distributions, short-tailed symmetric distributions, and asymmetric distributions are considered:

  • Long-tailed Symmetric
    • Laplace(0, 1)
    • \(t(5)\)
    • \(t(20)\)
  • Short-tailed Symmetric
    • Logistic(0, 1)
    • U(0, 1)
    • Beta(5)
  • Asymmetric
    • \(\chi^{2}(15)\)
    • Exp(1)
    • Lognormal(0, 1)

Only sample sizes \(n=10,20,30,50,75,100,200,500,1000, 1500, 2000\) are considered. The power of the test for the specified distributions and sample sizes are summarized below.

Figure 3. Power of SiNoTest
Figure 3. Power of SiNoTest


Based on Figure 3, the power of the test seems to be generally largest for asymmetric distributions (colored red). Relatively, the power seems to be moderate in general for short-tailed (colored green) and long-tailed (colored blue) symmetric distributions, except for \(t(20)\). The power of the test is low when the true distribution is \(t(20)\) - this is not surprising since the \(t\)-distribution approaches the standard normal distribution as the degrees of freedom approach \(+\infty\).

Complete Paper

The complete paper includes everything here, as well as additional details and analyses (such as the test’s robustness to non-sampling errors). It will be available to everyone for download soon.

References

  • Cao, R., & Lugosi, G. (2005, December). Goodness-of-Fit Tests Based on the Kernel Density Estimator. Scandinavian Journal of Statistics, 32(4), 599-616. https://www.jstor.org/stable/4616900

  • Lilliefors, H. W. (1967, June). On the Kolmogorov-Smirnov Test for Normality with Mean and Variance Unknown. Journal of the American Statistical Association, 62(318), 399-402. https://www.jstor.org/stable/2283970

  • Pavia, J. M. (2015). Testing Goodness-of-Fit with the Kernel Density Estimator: GoFKernel. Journal of Statistical Software, 66(1), 1-27. https://doi.org/10.18637/jss.v066.c01

  • Wijekularathna, D. K., Manage, A. B., & Scariano, S. M. (2019). Power analysis of several normality tests: A Monte Carlo Simulation Study. Communications in Statistics - Simulation and Computation, 51(3), 757-773. https://doi.org/10.1080/03610918.2019.1658780