Justine Lequesne Centre

Henri Becquerel

Philippe Regnault

vsgoftest: An R Package for Goodness-of-Fit Testing Based on Kullback-Leibler Divergence

Keywords: R statistical computing environment, goodness-of-fit tests, Shannon entropy, Kullback-Leibler divergence, sample spacing based estimation

), of various classical families of distributions. The theoretical framework of the so-called Vasicek-Song (VS) tests is summarized and followed by a detailed description of the different features of the package. The power and computational time performances of VS tests are studied through their comparison with other GOF tests. Application to real datasets illustrates the easy-to-use functionalities of the vsgoftest package.

Introduction

Goodness-of-fit (GOF) tests constitute a classical tool in deciding of the compatibility of data with a theoretical (probability) distribution. The present work proposes a package for the R statistical computing environment [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF] performing GOF tests based on Shannon entropy and Kullback-Leibler divergence, together with a methodological guide and applications.

Precisely, we consider fitting numeric (real valued) data either to a unique distribution, the so-called simple null hypothesis test H 0 : P = P 0 (θ) against H 1 : P ̸ = P 0 (θ), (1) or to a parametric family, the so-called composite null hypothesis test H 0 : P ∈ P 0 (Θ) against H 1 : P / ∈ P 0 (Θ).

(2)

The set P 0 (Θ) = {P ∈ D : P = P 0 (θ), θ ∈ Θ}, with Θ ⊂ R d , is a parametric subfamily of the set D of all probability distributions absolutely continuous with respect to Lebesgue measure on R, i.e., probability distributions with a density function. Decision is to be taken from the observation x n 1 = (x 1 , . . . , x n ) of a sample X n 1 = (X 1 , . . . , X n ) of size n of independent and identically distributed random variables drawn from P ∈ D.

Classically, a GOF test procedure is derived by computing some distance-like functional between the observations and the null distribution, or family of distributions, the null hypothesis being rejected when the distance is larger than a critical value. Kolmogorov-Smirnov, Cramérvon Mises and Anderson-Darling tests constitute some of the most commonly used GOF tests.

Their test statistics measure discrepancy between the empirical cumulative distribution function of the sample and the cumulative distribution function of the null distribution; these tests are refered to as EDF tests in the following; see [START_REF] Stephens | EDF statistics for goodness of fit and some comparisons[END_REF]. The tests in this paper are based on the Kullback-Leibler (KL) divergence of the density of the sample with respect to the null density.

GOF tests based on KL divergence have been introduced by [START_REF] Vasicek | A test for normality based on sample entropy[END_REF] for testing normality. Vasicek normality test relies on the maximum entropy property satisfied by the normal distribution: amoung all distributions with density and finite variance, Shannon entropy is maximized by the normal distribution. Vasicek test statistic is a monotone function of the entropy difference between the null normal distribution and the observed one. This has been subsequently extended to GOF tests of uncategorical data for numerous families of distributions satisfying a maximum entropy property -say maximum entropy (ME) distributions; see Section 2 for references. GOF tests based on entropy differences are known to have higher power than classical GOF tests in numerous cases; see Section 4.1. [START_REF] Song | Goodness-of-fit tests based on Kullback-Leibler discrimination information[END_REF] considers GOF tests based on KL divergence, for a large class of distributions including all classical distribution families. The test statistic is an estimate of the KL divergence between the sample and the null distributions. It is asymptotically normally distributed. When applied to ME distributions, it is equal to the difference between the entropies of the null distribution and the sample one, yielding the same decision rule as Vasicek test. This paper presents the implementation of Vasicek and Song tests (VS tests) for various families of distributions: uniform, normal, log-normal, exponential, gamma, Weibull, Pareto, Fisher, Laplace and Beta distributions. For further details on the theoretical aspects of VS tests, see [START_REF] Girardin | Entropy-based goodness-of-fit tests -a unifying framework. Application to DNA replication[END_REF] in which a unifying framework for tests based on entropy difference and KL divergence is provided.

Numerous R packages perform GOF tests for various families of distributions. The functions chisq.test, ks.test and shapiro.test of the stats package perform respectively the chi-squared test of adequacy to a discrete distribution, the Kolmogorov-Smirnov GOF test for any theoretical continuous distribution and the Shapiro-Wilk normality test. The packages goftest developed in [START_REF] Faraway | goftest: Classical Goodness-of-Fit Tests for Univariate Distributions[END_REF] and goft in Gonzalez-Estrada and Villasenor-Alva (2016) perform respectively Cramér-von Mises and Anderson-Darling GOF tests, and tests based on the ratios of variance and other moment estimators. KScorrect in Novack-Gottshall and Wang (2016) performs the Lilliefors-corrected Kolmogorov-Smirnov GOF test. Numerous GOF tests of the exponential or two-parameter Weibull distributions are available in EWGoF [START_REF] Krit | EWGoF: Goodness-of-Fit Tests for the Exponential and Two-Parameter Weibull Distributions[END_REF] while nortest and normtests are dedicated to testing normality. The dbEmpLikeGOF package developed in [START_REF] Miecznikowski | dbEmpLikeGOF: An R Package for Nonparametric Likelihood Ratio Tests for Goodness-of-Fit and Two-Sample Comparisons Based on Sample Entropy[END_REF] proposes GOF normality and uniformity tests based on empirical likelihood ratio. These tests are closely related with VS tests; similarities and differences between them are highlighted in the following sections.

The test procedure implemented in vsgoftest uses either the asymptotic distribution of the test statistic or Monte-Carlo simulation, depending on sample size or user's choice. Optional arguments are included for handling particular situations such as samples with numerous ties. They also contribute to make the procedure flexible and fully parameterizable. Besides these practical aspects, the paper presents a comprehensive review of the literature dealing with power properties of VS tests. Monte Carlo simulations are conducted to illustrate their performance when applied to discriminate between close distributions.

The paper is organized as follows. The theoretical framework of GOF tests based on Shannon entropy difference and KL divergence is briefly presented in Section 2. The functionalities of vsgoftest are presented in Section 3. The tests performed by vsgoftest are compared to other GOF tests in Section 4. More precisely, power comparisons to classical GOF tests are presented in Section 4.1; Section 4.2 focuses on the comparison of vsgoftest and dbEmpLikeGOF test procedures, which rely on very close theoretical frameworks but significantly differ in some of their features. Finally, applications to real data in Section 5 illustrate the usage of the proposed functionalities.

Entropy difference and KL divergence based GOF tests

The Shannon entropy of a distribution P with density function p on R has been defined in [START_REF] Shannon | A mathematical theory of communication[END_REF] as

S(P ) := - ∫ R p(x) log p(x)dx.
(3)

Entropy measures the uncertainty or variability of a distribution. The maximum entropy principle under moment constraints, or ME method, favours distributions with highest entropy for their highest degree of uncertainty; see [START_REF] Shannon | A mathematical theory of communication[END_REF] and [START_REF] Jaynes | Information theory and statistical mechanics[END_REF]. Among all distributions supported by a given finite length interval I in R, entropy is maximum and equals log |I| for the uniform distribution, where |I| denotes the length of I. Hence, the entropy difference log |I| -S(P ) can be thought as a distance-like measure between P and the uniform distribution.

Similarly, among all continuous distributions supported within R with mean µ and variance σ 2 , Shannon entropy is maximum for the normal distribution N (µ, σ 2 ) and equals

S(N (µ, σ 2 )) = ln(σ √ 2πe). ( 4
)
The entropy difference S(N (µ, σ 2 )) -S(P ) is nonnegative and thus defines a distance-like measure between any distribution with mean µ and variance σ 2 and the N (µ, σ 2 ) distribution.

Based on this property, [START_REF] Vasicek | A test for normality based on sample entropy[END_REF] derives a normality test, with a test statistic expressed in terms of entropy differences, defined as follows

K mn := n 2mS { n ∏ i=1 (X (i+m) -X (i-m) ) } 1/n = √ 2πe exp ( V mn -S ( N ( X, S 2 ))) ,
where

X := 1 n n ∑ i=1 X i and S 2 := 1 n n ∑ i=1 (X i -X) 2
are the empirical estimators of respectively the mean and the variance of the sample

X n 1 := (X 1 , . . . , X n ), X (1) ≤ • • • ≤ X (n)
denotes the order statistics associated to X n 1 and

V mn := 1 n n ∑ i=1 log ( n 2m [ X (i+m) -X (i-m) ] ) (5)
is the non-parametric Vasicek estimator of S(P ) based on spacings, with

X (i) = X (1) if i < 1 and X (i) = X (n) if i > n; the window size m ∈ N * is smaller than n/2.
The test statistic has been extended to various families of ME distributions under moment constraints; see [START_REF] Dudewicz | Entropy-based tests of uniformity[END_REF], [START_REF] Ebrahimi | Testing exponentiality based on Kullback-Leibler information[END_REF], [START_REF] Choi | Testing goodness-of-fit for laplace distribution based on maximum entropy[END_REF], [START_REF] Mergel | Test of goodness-of-fit for the inverse-gaussian distribution[END_REF], [START_REF] Mudholkar | An entropy characterization of the inverse Gaussian distribution and related goodness-of-fit test[END_REF], among many others.

A unifying framework for any exponential family of distributions is proposed in [START_REF] Girardin | Entropy-based goodness-of-fit tests -a unifying framework. Application to DNA replication[END_REF], with asymptotic properties, consistency and application to biology; see also [START_REF] Lequesne | Entropy-based goodness-of-fit test: Application to the Pareto distribution[END_REF], [START_REF] Lequesne | Tests statistiques basés sur la théorie de l'information[END_REF] and Lequesne (2015a) for power efficiencies, GOF tests of Pareto distributions and extension to generalized entropies.

The Kullback-Leibler (KL) divergence of a distribution P with respect to another one Q, is defined as

K(P |Q) := ∫ R p(x) log p(x) q(x) dx, ( 6
)
if P is absolutely continuous with respect to Q, with respective densities p and q, and as +∞ if not; see [START_REF] Kullback | On information and sufficiency[END_REF]. The KL divergence is linked to Shannon entropy through the relation

K(P |Q) = -S(P ) - ∫ R p(x) log q(x)dx. ( 7
)
The KL divergence is not a mathematical distance because of lack of both symmetry and triangular inequality, but it satisfies K(P |Q) ≥ 0, with K(P |Q) = 0 if and only if P = Q, and thus constitutes a natural measure of discrepancy for GOF tests. [START_REF] Song | Goodness-of-fit tests based on Kullback-Leibler discrimination information[END_REF] proposes GOF tests based on KL divergence for either simple (1) or composite (2) null hypothesis. Precisely, thanks to (7), the test statistic I mn is the estimator of K(P |P 0 (θ)), defined by

I mn := -V mn - 1 n n ∑ i=1 log p 0 (X i , θ n ), (8)
where V mn is the Vasicek estimator (5) of S(P ), and θ n is either the maximum likelihood estimator (MLE) of θ satisfying

1 n n ∑ i=1 log p 0 (X i , θ n ) = max θ∈Θ 1 n n ∑ i=1 log p 0 (X i , θ),
or θ itself in case of a simple null hypothesis (1).

The KL divergence K(P |Q) for a maximum entropy distribution Q under moment constraints reduces to the entropy difference S(Q)-S(P ) for all P satisfying the same moment constraints; see [START_REF] Csiszár | I-divergence geometry of probability distributions and minimization problems[END_REF]. This Pythagorean equality allowed [START_REF] Girardin | Entropy-based goodness-of-fit tests -a unifying framework. Application to DNA replication[END_REF]

K mn = √ 2πe exp(-I mn ),
yielding identical decision rules.

Based on the asymptotic properties of V mn proven by [START_REF] Dudewicz | Entropy-based tests of uniformity[END_REF] for testing uniformity, [START_REF] Song | Goodness-of-fit tests based on Kullback-Leibler discrimination information[END_REF] establishes the asymptotic behavior of I mn , independently of the null hypothesis: I mn is consistent and asymptotically normally distributed provided the null distribution belongs to the class

F = { P ∈ D : sup x: 0<F (x)<1 |p ′ (x)| p 2 (x) F (x)[1 -F (x)] < γ, γ > 0 } , ( 9
)
where F is the cumulative distribution function of P and p its density with derivative p ′ (almost everywhere). The class F contains the most classical distributions such as uniform (γ = 0), normal, exponential and gamma (γ = 1), Fisher (γ = (2 + ν 2 )/ν 2 where ν 2 is the second degree of freedom), Pareto (γ = (µ + 1)/µ, where µ is the shape parameter), etc. For

P 0 (Θ) ⊂ F , if m/ log n ---→ n→∞ 0 and m(log n) 2/3 /n 1/3 ---→ n→∞ 0, ( 10
) then √ 6mn[I mn -log(2m) + ψ(2m)] D -→ N (0, 1), ( 11
)
where ψ(m) is the digamma function. The asymptotic bias log(2m) -ψ(2m) of I mn is that of -V mn . [START_REF] Song | Goodness-of-fit tests based on Kullback-Leibler discrimination information[END_REF] suggests a bias correction in the asymptotic distribution ( 11) for moderate sample sizes:

√ 6mn [I mn -b mn ] D -→ N (0, 1), ( 12
)
where 12), an asymptotic p-value for the related VS test is given by

b mn := log(2m) -log(n) -ψ(2m) + ψ(n + 1) + 2m n R 2m-1 - 2 n m ∑ i=1 R i+m-2 , with R m := ∑ m j=1 1/j. From (
p = 1 -Φ -1 ( √ 6mn [I mn (x n 1 ) -b mn ] ) , ( 13
)
where I mn (x n 1 ) denotes the value of the statistic I mn for the observations x n 1 = (x 1 , . . . , x n ), and Φ denotes the cumulative distribution function of the normal distribution. According to [START_REF] Song | Goodness-of-fit tests based on Kullback-Leibler discrimination information[END_REF], the asymptotic p-value (13) provides accurate results for sample sizes n larger than 80. For small sample sizes, Monte Carlo simulations should be preferred. A large number N of replications of X n 1 drawn from P 0 ( θ n ) (or P 0 (θ) in case of simple null hypothesis) are generated. The test statistic I i mn is computed for each replication i, 1 ≤ i ≤ N. The p-value is then given by the empirical mean ( [START_REF] Song | Goodness-of-fit tests based on Kullback-Leibler discrimination information[END_REF] proposes to minimize I mn -that is to maximize V mn , with respect to m, yielding the most conservative test. The KL divergence K(P |P 0 (θ)) being nonnegative, values of m for which I mn is negative are excluded, leading to choose m subject to the constraint

∑ N i=1 1 {I i mn >Imn(x n 1 )} )/N. For choosing m,
V mn ≤ - 1 n n ∑ i=1 log p 0 ( .; θ n ) . ( 14
)
Finally, the window size proposed by [START_REF] Song | Goodness-of-fit tests based on Kullback-Leibler discrimination information[END_REF] -say the optimal window size, is

m := min { m * ∈ argmax m∈N * { V mn : V mn ≤ - 1 n n ∑ i=1 log p 0 ( X i , θ n ) } : 1 ≤ m * < ⌊n 1/3-δ ⌋ } , ( 15
) for some δ < 1/3 and the VS test statistic is then

I mn = -V mn - 1 n n ∑ i=1 log p 0 ( X i ; θ n ) . ( 16
)
The upper bound n 1/3-δ for the window size m is chosen so that conditions (10) are fulfilled and hence that asymptotic normality (11) holds. No optimal choice of δ exists; it depends on the family of distributions of the null hypothesis; see Section 3 for details.

The package vsgoftest presented below performs VS GOF tests for several parametric families of ME distributions: uniform, normal, log-normal, exponential, Pareto, Laplace, Weibull, Fisher, gamma and beta distributions. These families of distribution, all included in the class F given by ( 9), have been chosen so that the package covers a large variety of applications. Note that the package dbEmpLikeGOF performs uniformity and normality VS tests, with an alternative choice for the window size. Precisely, the test statistic is nI mn + 1/2, and the window m is chosen, between 1 and n 1/2 , minimizing nI mn . The constraint ( 14) is not considered. The asymptotic distribution of nI mn is not used, p-values being computed from a pre-calculated table for small sample sizes or via Monte-Carlo simulation; see [START_REF] Miecznikowski | dbEmpLikeGOF: An R Package for Nonparametric Likelihood Ratio Tests for Goodness-of-Fit and Two-Sample Comparisons Based on Sample Entropy[END_REF] and [START_REF] Vexler | Empirical likelihood ratios applied to goodness-of-fit tests based on sample entropy[END_REF]. This alternative methodological approach leads to different decisions that may be less reliable, particularly when applied to heavy tailed samples. Other differences in the coding structure make vsgoftest faster than dbEmpLikeGOF, especially when Monte-Carlo simulation is performed. These points will be detailled in Section 4.2.

The package vsgoftest

The vsgoftest package provides functions for estimating Shannon entropy of absolutely continuous distributions and testing the goodness-of-fit of some theoretical family of distributions to a vector of real numbers. It also provides functions for computing the density, cumulative density and quantile functions of Pareto and Laplace distributions, as well as for generating samples from these distributions.

The vsgoftest package is available on CRAN mirrors and can be installed by executing the command install.packages('vsgoftest')

Alternatively, the latest (under development) version of the vsgoftest package is also available and can be installed in R from the github repository of the project as follows:

#Package devtools must be installed devtools::install_github(repo = 'pregnault/vsgoftest')

The package is structured around two functions, entropy.estimate and vs.test; the first one computes the spacing based estimator (5) from a numeric sample, the second one performs Vasicek-Song GOF test for usual parametric families of distributions based on the test statistic (16). A comprehensive presentation of their usage is proposed in Sections 3.1 and 3.2, with numerous examples. Section 3.3 provides further technical information about the structure of the package.

Function entropy.estimate for estimating Shannon entropy

The function entropy.estimate computes the spacing based estimate (5) of Shannon entropy (3) from a numeric sample. Two arguments have to be provided:

• x: the numeric sample;

• window: an integer between 1 and half of the sample size, specifying the window size of the spacing-based estimator (5).

It returns the estimate of Shannon entropy of the sample. Here is an example for a sample drawn from a normal distribution with parameters µ = 0 and σ 2 = 1. The estimate returned by entropy.estimate obviously depends on the window selected by the user, as illustrated by the following chunck.

sapply(1:10, function(w) entropy.estimate(x = samp, window =w))

[1] 1.205018 1.346352 1.378732 1.387337 1.391691 1.393512 1.394428 [8] 1.394728 1.394486 1.392669

One may select the window size that maximizes the entropy estimate, as follows.

n <-100 #sample size V <-sapply(1:(n/2 -1), function(w) entropy.estimate(x = samp, window =w)) which.max(V) #Choose window that maximizes entropy

Function vs.test for testing GOF to a specified model

The function vs.test performs the VS test, as described in Section 2; setting two non-optional arguments is required:

• x: the numeric sample;

• densfun: a character string specifying the theoretical family of distributions of the null hypothesis. Available families of distributions are: uniform, normal, log-normal, exponential, gamma, Weibull, Pareto, Fisher and Laplace distributions. They are referred to by the symbolic name in R of their density function. For example, set densfun = 'dnorm' to test GOF of the family of normal distributions; see Table 1 for details.

It returns an object of class htest, i.e., a list whose main components are:

• statistic: the value of VS test statistic ( 16) for the sample, with optimal window size defined by (15);

• parameter: the optimal window size;

• estimate: the maximum likelihood estimate of the parameters of the null distribution (for the test (2) with composite null hypothesis);

• p.value: the p-value associated to the sample.

By default, vs.test performs the composite VS test of the family of distributions densfun for the sample x. The p-value is estimated by means of Monte-Carlo simulation if the sample size is smaller than 80, or through the asymptotic distribution (11) of the VS test statistic otherwise.

In the following example, a normally distributed sample is simulated. VS test rejects the null hypothesis that this sample is drawn from a Laplace distribution, but does not reject the normality hypothesis (for a significant level set to 0.05).

set.seed(5) samp <-rnorm(50,2,3) vs.test(x = samp, densfun = 'dlaplace')

Distribution Call (densfun) Parameters Density Default δ Uniform "dunif" a < b 1 b -a 1 1 [a,b] (x) 1/12 Normal "dnorm" µ ∈ R, σ ∈ R + 1 √ 2πσ exp ( - 1 2σ 2 (x -µ) 2 ) 1/12 Log-normal "dlnorm" µ ∈ R, σ ∈ R + 1 x √ 2πσ exp ( -(log x-µ) 2 2σ 2 ) 1 1 R+ (x) 1/12 Exponential "dexp" λ ∈ R λ exp(-λx)1 1 R+ (x) 1/12 Pareto "dpareto" c ∈ R, µ ∈ R + µc µ 1 x µ+1 1 1 [c,∞[ (x) 1/12 Laplace "dlaplace" µ ∈ R, σ ∈ R + 1 2σ exp ( - |x -µ| σ ) 1/12 Weibull "dweibull" a, b ∈ R + a b a x a-1 exp [ - ( x b ) a ] 1 1 R+ (x) 2/15 Fisher "df" ν 1 , ν 2 ∈ R + ( d1x d1x+d2 ) d1/2 ( 1 -d1x d1x+d2 ) d2/2 xB(d 1 /2, d 2 /2) 1 1 R+ (x) 2/15 Gamma "dgamma" α, β ∈ R + β α Γ(α) x α-1 exp(-βx)1 1 R+ (x) 2/15 Beta "dbeta" α, β ∈ R x α-1 (1 -x) β-1 B(α, β) 1 1 R+ (x) 2/15
Table 1: Families of distributions supported as the null by the vsgoftest package. The column denoted "Default δ" corresponds to the default setting for the parameter δ; see (15 For performing a simple null hypothesis GOF test, the additional argument param has to be set to a numeric vector, consistent with the parameter requirements for the null distribution. In such case, the MLE of the parameter(s) of the null distribution has not to be computed and hence the component estimate in results is not available. Vasicek's estimates V mn are computed for all m from 1 to n 1/3-δ , where δ < 1/3; the test statistic is I mn for m the optimal window size, as defined in (15). The choice of δ depends on the family of distributions of the null hypothesis. Precisely, for Weibull, Pareto, Fisher, Laplace and Beta, δ is set by default to 2/15, while for uniform, normal, log-normal, exponential and gamma, it is set to 1/12. These default settings result from numerous experimentations. Still, the user can choose another value through the optional argument delta. Note that upper-bounding the window size by n 1/3-δ is only required when the asymptotic normality of I mn is used to compute asymptotic p-values from (11). When the p-value are computed by means of Monte-Carlo simulation, this upper-bound can be extended to n/2 by adding extend = TRUE, which may lead to a more reliable test, as illustrated below. Enlarging the range of m is also pertinent if ties are present in the sample. Indeed, the presence of ties is particularly inappropriate for performing VS tests, because some spacings X (i+m) -X (i-m) can be null. The window size m has thus to be greater than the maximal number of ties in the sample. Hence, if the upper-bound n 1/3-δ is less than the maximal number of ties, the test statistic can not be computed. Setting extend to TRUE can avoid this behavior, as illustrated below. Finally, Vasicek's estimate V mn may exceed the parametric estimate of the entropy of the null distribution for all m between 1 and n 1/3-δ . Then, no window size exists satisfying (15), as illustrated below.

). Note that B(a, b) = ∫ 1 0 x a-1 (1 -x) b-1 dx is the Beta function.
set.seed(84) ech <-rpareto(20, mu = 1/2, c = 1) vs.test(x = ech, densfun = 'dpareto', param = c(1/2, 1))

Error in vs.estimate(x, densfun, ESTIM, extend, delta, relax): The sample entropy is greater than empirical maximal entropy for all possible window sizes; the sample may be too small or is unlikely to be drawn from the null distribution.

Enlarging the possible window sizes by setting extend to TRUE may enable Vasicek estimates to be smaller than empirical entropy.

Note that when computing the p-value by Monte-Carlo simulation, the constraint (14) may not be satisfied for some replicates, whatever be the window size. These replicates are then ignored and the p-value is computed from the remaining replicates. A warning message is added to the output, informing on the number of ignored replicates. data(contaminants) #load data from package vsgoftest; see ?contaminants set.seed(1) vs.test(x = aluminium2, densfun = 'dpareto') Warning in vs.test(x = aluminium2, densfun = "dpareto"): For 176 simulations (over 5000 ), entropy estimate is greater than empirical maximum entropy for all window sizes.

Vasicek-Song GOF test for the Pareto distribution data: aluminium2 Test statistic = 1.3676, Optimal window = 2, p-value < 2.2e-16 sample estimates: mu c 0.3288148 360.0000000

A large proportion of such ignored replicates may indicate that the original sample is too small or the null distribution does not fit it.

The function vs.test also allows to avoid the constraint ( 14) when computing the optimal window size, by setting the optional argument relax to TRUE. This however should be used with special care, even when the p-value is computed by Monte-carlo simulation, because it may lead to spurious conclusions. Some examples will be discussed in Section 5. This option is to recover the non-parametric likelihood ratio GOF test developed by [START_REF] Vexler | Empirical likelihood ratios applied to goodness-of-fit tests based on sample entropy[END_REF] and performed by dbEmpLikeGOF; see Section 4.2.

Technical information on the internal structure of the vsgoftest package

While entropy.estimate is a stand-alone function -depending only on the base and stats packages, vs.test is supported by a set of internal functions -not available for users; the structure of the package and connections between functions are described in the organisational chart presented in Figure 1. Functions available for users are depicted by rectangles while internal functions are depicted by ellipses. An arrow connecting a function to another means that the first function (say master function) calls the second (slave) during execution. When such a call is optional (depending on arguments given in the master function), the arrow is dashed and annotated with the corresponding argument settings. The function fitdist depicted by a dashed rectangle is a function implemented in the fitdistrplus package [START_REF] Delignette-Muller | fitdistrplus: An R Package for Fitting Distributions[END_REF]. The double-lined ellipse depicts a C++ encoded function that has been integrated via the package Rcpp [START_REF] Eddelbuettel | Rcpp: Seamless R and C++ Integration[END_REF].

The vsgoftest package is structured in such a way so as to:

• Allow easy access to the code source. Especially, the master function vs.test calls four slave functions corresponding to the following tasks (enumerated according to the organisational chart of Figure 1):

1. computing the MLE of the parameter θ of the null distribution through the function

MLE.param.

2. Computing Vasicek estimate V mn of Shannon entropy for the sample with the optimal window m given by ( 15). (2010), POT in [START_REF] Ribatet | POT: Generalized Pareto Distribution and Peaks Over Threshold. R package version 1.1-6[END_REF] and smoothmest in [START_REF] Hennig | smoothmest: Smoothed M-estimators for 1-dimensional location. R package version 0.1-2[END_REF]. The MLE θ n of the parameter of the null distribution is computed thanks to the function fitdist of the fitdistrplus package only if no closed form expression is known for it, i.e., for Gamma, Weibull, Beta and Fisher distributions. Otherwise, the closed form expression is used.

Computing

• Optimize time and resources, especially for Monte-Carlo simulation. To this end, the most time-consuming part of the procedure -namely, the computation of Vasicek estimate for all possible window sizes, has been converted to C++ and integrated to the package via Rcpp, in the internal function vestimates.

Performance of Vasicek-Song tests

First, a review of power studies of VS tests available in literature is presented in Section 4.1. Then, power comparisons of VS tests and classical GOF tests are proposed when applied to discriminate between close distributions, such as Pareto versus shifted log-normal and Exponential versus Weibull. Finally, the features of packages vsgoftest and dbEmpLikeGOF are compared in Section 4.2; the methodological differences are highlighted, the higher performance of vsgoftest both in terms of power and computational time is pointed out and illustrated.

Power computation

Comparisons of the power properties of VS tests are widely discussed in the literature. Various choices of null and alternative distribution families are considered. VS tests are shown to generally outperform classical GOF tests. A comprehensive list of these references is given in this section, with main conclusions summarized in Table 2. Especially, power properties of the VS test for normality have been discussed by [START_REF] Vasicek | A test for normality based on sample entropy[END_REF], [START_REF] Arizono | A test for normality based on Kullback-Leibler information[END_REF] and [START_REF] Gurevich | Standardized forms of Kullback-Leibler information based statistics for normality and exponentiality[END_REF] among many others. Compared with many tests, including Kolmogorov-Smirnov (KS), Cramér-von Mises (CvM), Anderson-Darling (AD) and Shapiro-Wilk (SW), the VS test exhibits higher power for most of alternative distributions. When the null distribution is an exponential distribution, the VS test is also shown in [START_REF] Ebrahimi | Testing exponentiality based on Kullback-Leibler information[END_REF] to be more powerful than the Van-Soest and Finkelstein and Schafer tests, which are modified versions of respectively CvM and KS tests, for various alternative distributions such as Weibull, gamma and log-normal. [START_REF] Choi | Testing goodness-of-fit for laplace distribution based on maximum entropy[END_REF] for Laplace and [START_REF] Lequesne | Tests statistiques basés sur la théorie de l'information[END_REF] for Pareto show that the VS test is more powerful than EDF tests, for various alternative distributions. The uniform VS GOF test is shown to outperform many other tests for alternative distributions having most of their mass near 0.5, but remains less powerful than CvM and Watson tests for other alternative distributions.

On the basis of power computation in the literature, we choose to compare the power of the VS test to the KS, CvM and AD tests, for close null and alternative distributions. In particular, difficulties in distinguishing a Pareto tail from that of a log-normal is an issue; see for example [START_REF] Malevergne | Testing the Pareto against the lognormal distributions with the uniformly most powerful unbiased test applied to the distribution of cities[END_REF]. For illustration, we estimate through Monte-Carlo simulation the power of VS, KS, CvM and AD of Pareto distributions applied to samples drawn from a (shifted) log-normal distribution. We simulate 10000 replicates x n 1 of a random sample X n 1 drawn from a shifted log-normal distribution LN (0, σ) with support [1, ∞[ and σ = 1, 1.25, for n ∈ {20, 30, 50, 100}. Then, we apply the tests for the simple null hypothesis H 0 : P = Par(1, µ) , for µ = 1 when σ = 1 and µ = 0.8 when σ = 1.25; the power is estimated by the proportion of rejections of the null hypothesis among the 10000 replicates. The following code chunk illustrates the procedure, for σ = µ = 1 and n = 20, using the VS test. This procedure immediately adapts to other values of σ, µ and n and to other tests 1 . Results are presented in Table 3

(top). N <-10000 n <-20 mu <-1 set.seed(54)
1 The seed of the pseudo-random number generator has been changed for each couple of µ and n; the whole procedure yielding Table 3

vsgoftest versus dbEmpLikeGOF for testing uniformity and normality

The package dbEmpLikeGOF in [START_REF] Miecznikowski | dbEmpLikeGOF: An R Package for Nonparametric Likelihood Ratio Tests for Goodness-of-Fit and Two-Sample Comparisons Based on Sample Entropy[END_REF] performs uniformity and normality tests based on empirical likelihood ratios (ELR) -say ELR tests. These tests are strongly linked to VS tests. Precisely, for testing the normality of a sample X 1 , . . . X n , the ELR test statistic is log V n , where

V n = min 1≤m<n 1/2   (2πe Š2 ) n/2 n ∏ i=1 2m n [ X (i+m)-X (i-m) ]   and Š2 = 1 n -1 n ∑ i=1   X i - 1 n n ∑ j=1 X j   2 .
Mere algebra yields log

V n = nI mn + 1 2 , ( 17

Application to real data

The vs.test package contains environmental data originating from a guidance report edited by the Technology Support Center of the United States Environmental Protection Agency; see [START_REF] Singh | The lognormal distribution in environmental applications[END_REF]. According to [START_REF] Singh | The lognormal distribution in environmental applications[END_REF], environmental scientists take remediation decisions at suspected sites based on organic and inorganic contaminant concentration measurements. These decisions usually derive from the computation of confidence upper bounds for contaminant concentrations. Testing the goodness-of-fit of specified models hence appears of prior interest. [START_REF] Singh | The lognormal distribution in environmental applications[END_REF] also points out that contaminant concentration data from sites often appear to follow a skewed probability distribution, making the log-normal family a frequently-used model. The authors illustrate their purpose by applying Shapiro-Wilk test to the log-transformed of the samples aluminium1, manganese, aluminium2 and toluene (stored in the present package)5 ; see the empirical skewness computed in the following chunk.

data(contaminants) #Load environmental data from package #Package DescTools required for this chunk unlist(lapply(X = list(aluminium1, manganese, aluminium2, toluene), FUN = DescTools::Skew))

[1] 2.323343 1.698686 1.996607 3.961129

The following code chunks intend to illustrate the use and behavior of the function vs.test for these environmental data. The significant level is fixed to 0.1 as in [START_REF] Singh | The lognormal distribution in environmental applications[END_REF].

Note that warning messages notifying that there are ties in the samples have been dropped out from outputs.

Conclusion

Vasicek-Song tests constitute powerful GOF tests for classical parametric families of distributions, relying on an information theoretical framework. They can be easily performed by using the vsgoftest package for R. Default and optional settings of the functions provided by the package make the procedure both intuitive and flexible. Its application to real datasets manages to illustrate its practical usage.

The package allows for testing GOF of a significant list of parametric models; this list could be extended in further releases. New entropy-based GOF tests could also be considered by using Rényi entropy and divergence -see Lequesne (2015a), thus extending even more the class of possible distributions, e.g., Student distributions.

  library('vsgoftest') Loading required package: fitdistrplus Loading required package: MASS Loading required package: survival set.seed(2) #set seed of PRNG samp <-rnorm(n = 100, mean = 0, sd = 1) #sampling from normal distribution entropy.estimate(x = samp, window = 8) #estimating entropy with window = 8 [1] 1.394728 log(2*pi*exp(1))/2 #the exact value of entropy [1] 1.418939

  Let us consider a sample drawn from a Pareto distribution with density p(x; c, µ) = µc µ x µ+1 , x ≥ c, where c > 0 and µ > 0, which can be obtained by making use of the function rpareto as illustrated below. Its Shannon entropy is S(p(.; c, µ)) = -ln µ + ln c + 1 µ + 1. vsgoftest: package for GOF testing based on KL-divergence set.seed(5) n <-100 #Sample size samp <-rpareto(n, c = 1, mu = 2) #sampling from Pareto distribution entropy.estimate(x = samp, window = 3) [1] 0.8480204 -log(2) + 3/2 #Exact value of entropy [1] 0.8068528

  vs.test(x = samp, densfun = 'dnorm') Vasicek-Song GOF test for the normal distribution data: samp Test statistic = 0.21655, Optimal window = 2, p-value = 0

  vsgoftest: package for GOF testing based on KL-divergence set.seed(26) vs.test(x = samp, densfun = 'dnorm', param = c(2,3)) Vasicek-Song GOF test for the normal distribution with Mean=2, St. dev.=3 data: samp Test statistic = 0.22196, Optimal window = 2, p-value = 0.331 If param is not consistent with the specified distribution -e.g., standard deviation for testing a normal distribution is missing or negative, the execution is stopped and an error message is returned. set.seed(2) samp <-rnorm(50, -2, 1) vs.test(samp, densfun = 'dnorm', param = -2) Error in vs.test(samp, densfun = "dnorm", param = -2): "param": invalid parameter (not consistent with the specified distribution) One can estimate the p-value of the sample by Monte-Carlo simulation, even when sample size is larger than 80, by setting the optional argument simulate.p.value to TRUE (NULL by default). The number of Monte-Carlo replicates can be fixed through the optional argument B (default is B = 5000). set.seed(1) samp <-rweibull(200, shape = 1.05, scale = 1) set.seed(2) vs.test(samp, densfun = 'dexp', simulate.p.value = TRUE, B = 10000) Vasicek-Song GOF test for the exponential distribution data: samp Test statistic = 0.10907, Optimal window = 3, p-value = 0

  set.seed(63) vs.test(samp, densfun = 'dexp', delta = 5/30) Vasicek-Song GOF test for the exponential distribution data: samp Test statistic = 0.16517, Optimal window = 2, p-value = 0

  set.seed(8) samp <-rexp(30, rate = 3) vs.test(x = samp, densfun = "dlnorm") Vasicek-Song GOF test for the log-normal distribution data: samp Test statistic = 0.30717, Optimal window = 2, p-value = 0.1206 sample estimates: Location Scale -2.162290 1.683868 vs.test(x = samp, densfun = "dlnorm", extend = TRUE) Vasicek-Song GOF test for the log-normal distribution data: samp Test statistic = 0.3029, Optimal window = 3, p-value = 0

  samp <-c(samp, rep(4,3)) #add ties in the previous sample vs.test(x = samp, densfun = "dexp") Warning in vs.estimate(x, densfun, ESTIM, extend, delta, relax): Ties should not be present for Vasicek-Song test Error in vs.estimate(x, densfun, ESTIM, extend, delta, relax): Too many ties to compute Vasicek estimate. vs.test(x = samp, densfun = "dexp", extend = TRUE) Warning in vs.estimate(x, densfun, ESTIM, extend, delta, relax): Ties should not be present for Vasicek-Song test Vasicek-Song GOF test for the exponential distribution data: samp Test statistic = 0.025702, Optimal window = 16, p-value = 0

Figure 1 :

 1 Figure 1: Organisational chart of the structure of the package and connections between functions available for users (depicted by rectangles) and internal functions (depicted by ellipses).

  simulate.p.value = TRUE, B = 1000)$p.value) The power of VS, KS, CvM and AD tests is similarly computed for null exponential and alternative Weibull distributions. The Weibull distribution W(a, b) reduces to an exponential distribution E(1/a) when b = 1. The main aim is thus to determine which test better discriminates between these distributions when the shape parameter of the Weibull distribution is close to 1, precisely b = 1.2 and b = 1.3. Results are given in

Figure 2 :

 2 Figure 2: Distributions (violin plots) of computation time of 100 iterations of the functions dbEmpLikeGOF (from package dbEmpLikeGOF) and vs.test for testing normality of a sample of size 50 drawn from the standard normal distribution (top) and for testing uniformity of a sample of size 50 drawn from the uniform distribution on [0, 1] (bottom). The p-values are estimated by Monte-Carlo simulation, based on 1000 replicates. The red dashed lines are the median computation times.

  set.seed(1) vs.test(x = aluminium1, densfun = 'dlnorm')Vasicek-Song GOF test for the log-normal distribution data: aluminium1 Test statistic = 0.31232, Optimal window = 2, p-value = 0The log-normal hypothesis is not rejected for aluminium1. Similar results are obtained for manganese. Log-normality is rejected for aluminium2. set.seed(1) vs.test(x = aluminium2, densfun = 'dlnorm') Vasicek-Song GOF test for the log-normal distribution set.seed(1) vs.test(x = aluminium2, densfun = 'dpareto') Vasicek-Song GOF test for the Pareto distribution data: aluminium2 Test statistic = 1.3676, Optimal window = 2, p-value < 2.2e-16 sample estimates: mu c 0.3288148 360.0000000 Applying vs.test to toluene with default settings yields no result because of numerous ties and the violation of (14). Uniformity of the sample transformed by the cumulative density function of the Pareto distribution can be tested as follows. Goodness-of-fit of the Pareto distribution is not rejected for toluene. #Compute the MLE of parameters of Pareto dist. res.test <-vs.test(x = toluene, densfun = 'dpareto', extend = TRUE, relax = TRUE) #Test uniformity of transformed data set.seed(5) vs.test(x = ppareto(toluene, mu = res.test$estimate[1], c = res.test$estimate[2]), densfun ='dunif', param = c(0,1), extend = TRUE) Vasicek-Song GOF test for the uniform distribution with Min=0, Max=1 data: ppareto(toluene, mu = res.test$estimate[1], c = res.test$estimate[2]) Test statistic = 0.25383, Optimal window = 10, p-value = 0.2496

  is available in the file vsgoftest performances.R, in the directory inst/doc of the package source file.

	Reference			Null distrib.	Alt. distrib.			VS compared with Most powerful
	Vasicek (1976)			exponential, gamma, uni-	KS, CvM, Kuiper,	AD (for Cauchy), VS
				Normal	form, beta, Cauchy			Watson, AD, SW	(for others)
	Arizono and Ohta		log-normal, uniform, χ 2 ,	KS, CvM, χ 2	VS
	(1989)				student
	Gurevich and David-		log-normal, χ 2 , Student,	KS	KS (for Student and
	son (2008)				uniform,	exponential,	Cauchy), VS (for
					gamma, beta, Cauchy		others)
	Ebrahimi	et	al.	Exponential	Weibull, gamma, log-normal Van-Soest,	Finkel-	VS
	(1992)									stein and Schafer
	Dudewicz and Van	Uniform	Distributions defined on	KS, CvM, Kuiper,	VS (for alternative
	Der Meulen (1981)		[0, 1]					Watson, AD, log-	having most of its
										statistic, χ 2	mass near 0.5), CvM
										or Watson (for oth-
										ers)
	Choi and Kim (2006) Laplace	Normal, Student, logis-	KS,	CvM,	AD,	VS
					tic,	Cauchy,	uniform,	Kuiper and Watson
					chi-squared, Weibull, log-
					normal, extreme value and
					inverse Gaussian
	Lequesne (2015b)		Pareto	Weibull,	gamma,	log-	KS, AD	VS
					normal,	two-parameter
					exponential
	Mudholkar and Tian	Inv. Gaussian	exponential,		uniform,	KS	VS (uniform and
	(2002)				Weibull and log-normal		Weibull)
	Alizadeh Noughabi,	Rayleigh	Weibull,	gamma,	log-	KS,	CvM,	AD,	VS (for uniform) and
	Alizadeh Noughabi,		normal, half-normal, uni-	Kuiper and Watson	AD (for other alter-
	and	Ebrahimi		form,	modified extreme	natives)
	Moghaddam	Be-		value,	linear increasing
	habadi (2014)			failure rate law, Dhillon's
					law and Chen's distribution
	Perez-Rodriguez,		Gumbel	Weibull, log-normal, nor-	KS,	CvM,	AD,	AD (for heavy tails),
	Vaquera-Huerta,			mal, logistic, Cauchy, Stu-	Kuiper and Kinnison	VS (for others)
	and Villaseñor-Alva		dent, gamma and Fréchet
	(2009)
	Tsujitani, Ohta, and	Extreme-Value Normal and 3-parameter	KS,	CvM,	AD,	VS
	Kase (1980)			log-Weibull				Kuiper, Mann et
										modified Mann test
	Lund and Rao Jam-	von Mises	Mixtures of Von Mises dis-	Watson and inte-	Watson (for long-
	malamadaka (2000)		tributions (bimodal, skewed,	grated squared error	tailed), VS (for half)
					long-tailed and half), the	test	and fairly equal for
					cardoid and triangular dis-	other alternatives
					tributions

Table 2 :

 2 Power studies of VS tests performed in the literature.

	res.pow <-replicate(n = N,
	expr = vs.test(x = 1 + rlnorm(n,
	meanlog = 0,
	sdlog = 1),
	densfun = 'dpareto',
	param = c(1,1),

Table 3 :

 3 Power (expressed as percentage of true rejection) of VS, KS, CvM and AD tests for testing Pareto and exponential distributions against shifted log-normal (top) and Weibull distributions (bottom).

Table 3 (

 3 bottom), clearly showing that the VS test outperforms EDF tests. Note that the above procedure for comparing the power of GOF tests adapts easily to other sets of null and alternative distributions.

Table 4 :

 4 Power comparisons between VS and ELR normality tests for samples drawn from a Laplace distribution (left) and a Student distribution (right). Power is estimated by means of Monte-Carlo simulation based on 1000 replicates of the samples. and the constraint (14) is not taken into account by the EL test. Enlarging the upper bound from n 1/4 to n 1/2 may lead to a more powerful decision rule, as mentioned and illustrated in Section 3.2. Still, practically, the normality VS test tends to outperform the ELR test when applied to heavy tailed samples, as illustrated by Table42 . Moreover, the upper bound n 1/4 legitimates the use of the asymptotic distribution of I mn in vs.test, which is not performed by dbEmpLikeGOF. As previously mentioned in Section 3.2, disabling the constraint (14) may lead to spurious conclusions.The ELR test can be performed using vs.test, by suitably setting its arguments, as illustrated by the following code chunk 3 . Additionnaly, from a computational view point, some significant differences exist between the functions dbEmpLikeGOF and vs.test. Precisely, the p-value returned by dbEmpLikeGOF is computed by default by linear interpolation from a table of pre-computed p-values for various sample sizes (from n = 10 to 10000) and test statistic values; the p-value can be approximated by Monte-Carlo simulation by setting pvl.Table = FALSE. By default, vs.test computes the p-value by means of Monte-Carlo simulation or uses the asymptotic distribution (12), depending on the sample size. In both cases, vs.test is approximately five times faster than dbEmpLikeGOF, as illustrated by Figure 2 4 .

	VS ELR	VS ELR
	H 1 : L(0, 1)	H 1 : S(4)
	n = 50 17.8 16.0 16.1 14.6
	n = 200 86.2 64.9 71.2 35.8

) with m ∈ argmax 1≤m<n 1/2 V mn . Hence, ELR and VS tests differ only in the window size choice: the upper bound is n 1/2 for the ELR test while it is (by default) n 1/4 for the VS test

The comparison procedure is available in the file vsgoftest performances.R, in the directory inst/doc of the package source file.

Some slight difference remains between the two computed values, due to numerical inaccuracy in computation procedures: the estimated entropy of the null distribution is computed from the closed form expression (4) in dbEmpLikeGOF while it is computed as the empirical mean of the log-likelihood of the sample in vs.test.

Simulations have been performed on a Dell Lattitude E5580 laptop, equipped with an Intel ® Core™ i7-7600U CPU at 2.80GHz x 4, with 16GB RAM. R code for generating these figures is available in the file vsgoftest performances.R, in the directory doc of the package source file.

A succinct description of these data is available by executing the following R command: ?contaminants

data: aluminium2

Test statistic = 0.48369, Optimal window = 2, p-value = 0.0256 sample estimates: Location Scale 8.9273293 0.8264409 Due to numerous ties in toluene, vs.test can not compute Vasicek entropy estimate unless extend is set to TRUE. Still, vs.test notifies that the constraint ( 14) is violated for all window sizes, which suggests that data are not likely to be drawn from the log-normal distribution; see Section 2. Turning relax to TRUE yields the following result. The log-normal hypothesis is not rejected for aluminium1 and manganese while it is rejected for aluminium2 and toluene. These results are consistent with those obtained by [START_REF] Singh | The lognormal distribution in environmental applications[END_REF]. Further, the goodness-of-fit to the Pareto distributions is performed for aluminium2 and toluene. Log-normal and Pareto distributions usually compete with closely related generating processes and hard to distinguish tail properties; see for example [START_REF] Malevergne | Testing the Pareto against the lognormal distributions with the uniformly most powerful unbiased test applied to the distribution of cities[END_REF]. Goodness-of-fit of Pareto distribution is rejected for aluminium2.