Installation

Install PBGoF 0.2.0 from a local source archive with:

install.packages("PBGoF_0.2.0.tar.gz", repos = NULL, type = "source")

PBGoF provides goodness-of-fit procedures, robust parameter fitting, and graphical diagnostics for skew-normal and skew-t models. Version 0.2.0 adds skew-t fitting, parameter-bootstrap KS and CvM tests, and a corresponding graphical diagnostic to the package’s skew-normal procedures. After installation, load the package with:

library(PBGoF)

Overview

PBGoF assesses whether a numeric sample is compatible with a fitted univariate skew-normal (SN) or skew-t (ST) distribution. In both cases the parameters are estimated from the same observations used to calculate the empirical distribution function (EDF) statistic. This makes the null hypothesis composite, so an ordinary Kolmogorov–Smirnov p-value for a completely specified distribution is not appropriate.

The package offers two complementary approaches:

1. Parametric bootstrap tests simulate a reference distribution for the observed sample and re-estimate all parameters in every replicate.
2. SN precomputed-quantile tests use tables obtained from 100,000 Monte Carlo replicates for each combination of sample size and centered skewness represented in the bundled tables.

Parametric bootstrap KS and CvM tests are available for both families. The precomputed tables apply only to SN. Parameter estimation is performed by sn.fit.robust() for SN and st.fit.robust() for ST; each function uses a sequence of increasingly stabilized fitting methods.

Distribution models

gamma1_used <- round(abs(gamma1_hat), 2)

Robust parameter estimation

library(PBGoF)

set.seed(2026)
x <- sn::rsn(100, xi = 0, omega = 1, alpha = 4)
x_st <- sn::rst(150, xi = 0, omega = 1, alpha = 3, nu = 8)

sn.fit.robust(x, para_form = "DP")
sn.fit.robust(x, para_form = "CP")
st.fit.robust(x_st, para_form = "DP")

sn.fit.robust() first attempts ordinary maximum likelihood estimation (MLE). If the fit does not yield finite estimates and standard errors, it tries maximum penalized likelihood estimation (MPLE) with the default penalty and then MPLE with the matching-prior penalty. If every attempt fails, the function returns a named vector of missing values and issues a warning.

st.fit.robust() applies the same idea to the four-parameter skew-t family. It tries MLE first and then Q-penalized MPLE with nlminb and BFGS fallbacks. The DP output contains xi, omega, alpha, and nu, together with their standard errors. The CP output contains mean, standard deviation, skewness, and kurtosis when those moments exist.

Here, robust refers to protection against numerical fitting failures. It does not mean that the estimator is resistant to outliers or contamination. Users should still inspect their data and assess whether the selected distribution is a scientifically reasonable model.

The input must be a numeric vector containing at least 10 finite observations and at least two distinct values. Matrices, factors, missing values, and infinite values are rejected rather than silently converted.

Graphical assessment of the fitted models

sn.plot.check() provides a quick visual comparison between the observed sample and its fitted skew-normal distribution. It draws a density-scale histogram, adds the fitted SN density, and optionally marks individual observations with a rug plot.

The following example generates a sample directly with the sn package and then checks the fitted model:

set.seed(123)
x_plot <- sn::rsn(
  n = 200,
  xi = 1,
  omega = 2,
  alpha = 5
)

plot_result <- sn.plot.check(x_plot)
plot_result$parameters

The function fits the model through sn.fit.robust(), so its parameter estimates use the same numerical fallback strategy as the formal PBGoF tests. It returns the fitted DP parameters and plotted coordinates invisibly, allowing the underlying values to be inspected or tested programmatically.

Visual inspection can reveal features that a single p-value does not describe, including isolated outliers, multimodality, tail discrepancies, and systematic differences between the histogram and fitted density. The appearance of a histogram depends on its breaks, so it is often useful to compare several choices:

sn.plot.check(x_plot, breaks = "FD")
sn.plot.check(x_plot, breaks = "Scott")
sn.plot.check(x_plot, breaks = 20)

For a skew-t fit, st.plot.check() uses the same visual design and calls st.fit.robust() internally:

st.plot.check(x_st, breaks = "FD")

The plots are diagnostic aids rather than formal decision rules. Use sn.plot.check() with SN tests and st.plot.check() with ST tests.

Parametric bootstrap tests

The bootstrap functions are:

sn.para.bootstrap.ks.test(x, B = 1000, seed = 103)
sn.para.bootstrap.cvm.test(x, B = 1000, seed = 103)
st.para.bootstrap.ks.test(x_st, B = 1000, seed = 103)
st.para.bootstrap.cvm.test(x_st, B = 1000, seed = 104)

For each test, PBGoF performs these steps:

1. Fit the selected SN or ST distribution to the observed data.
2. Calculate the observed EDF statistic using the fitted DP parameters.
3. Generate B samples of the same size from the fitted model.
4. Refit the selected distribution separately in every bootstrap sample; the ST tests re-estimate all four parameters, including nu.
5. Recalculate the statistic using each bootstrap sample’s fitted parameters.
6. Estimate the p-value from the proportion of valid bootstrap statistics at least as large as the observed statistic.

The finite-simulation correction adds one to both the exceedance count and the number of valid replicates. It prevents a p-value of exactly zero. Failed fits are excluded. The returned numeric p-value has attributes recording the observed statistic, requested number of replicates, and numbers of valid and failed replicates:

p <- sn.para.bootstrap.ks.test(x, B = 999, seed = 103)
p
attributes(p)

The seed argument makes the bootstrap reproducible. PBGoF restores the caller’s previous random-number state when the test finishes. Set seed = NULL to continue from the current random-number stream.

Larger values of B give finer and more stable p-values but require more computation. Values such as 999 or 1999 are useful during analysis; substantially larger values may be preferable for final inference near a decision threshold.

Skew-normal tests based on precomputed quantiles

The fast lookup functions are:

PBGoF_ks_test(x)
PBGoF_cvm_test(x)

They avoid running a new bootstrap for each dataset. For observed data, the functions:

1. estimate DP parameters for evaluating the fitted distribution;
2. estimate CP parameters and calculate abs(gamma1), so negative and positive fitted skewness of the same magnitude use the same reference distribution;
3. round the absolute skewness to two decimal places and restrict it to the table range 0.01–0.99;
4. select the table row for the sample size and matched skewness;
5. compare the observed statistic with the stored 0.01–0.99 quantiles.

The bundled tables cover sample sizes 30 through 500. For a sample larger than 500, all observations are retained when fitting the model and constructing the EDF, but the lookup sample size and external statistic scaling are set to 500. The result reports both values:

set.seed(1)
x_large <- sn::rsn(750, alpha = 3)
result <- PBGoF_ks_test(x_large)

result$n       # 750: actual number of observations
result$n_used  # 500: scaling and lookup value

Mateu-Figueras, Puig, and Pewsey (2007) reported that the quantiles of the EDF statistics for sample sizes above 500 were almost identical to those for n = 500, and recommended using the n = 500 critical values. Replacing the external scaling value by 500 is the PBGoF convention used to remain consistent with the construction of its bundled tables; the full empirical distribution always uses all observations.

The lookup result contains:

• statistic: observed scaled test statistic;
• n: actual sample size;
• n_used: sample size used for scaling and table lookup;
• gamma1_hat: fitted signed CP skewness;
• gamma1_used: rounded absolute skewness used for lookup;
• p.value: approximate upper-tail p-value.

Because the stored probability grid advances in increments of 0.01, lookup p-values are conservative step-function approximations restricted to 0.01–0.99. They should not be interpreted as having greater precision than the table permits.

Advanced users can supply a compatible custom table:

PBGoF_ks_test(x, ks_table = my_ks_table)
PBGoF_cvm_test(x, cvm_table = my_cvm_table)

A custom table must be a data frame containing n, gamma1, and probability columns named q_0.01 through q_0.99. Quantiles within a selected row must be finite and nondecreasing.

Choosing between the approaches

Use the parametric bootstrap when maximum fidelity to the observed sample size is important, when the sample size is below the table range, or when a more finely resolved p-value is needed. Its main cost is computation because every bootstrap sample must be refitted.

Use an SN precomputed-quantile test for rapid repeated screening when the table’s sample-size and skewness grid is appropriate. It is especially useful when many datasets must be tested, but its p-values are discretized and depend on the simulation design used to create the tables. Precomputed tables are available only for SN; use parameter bootstrap for ST.

For important conclusions, it can be useful to run the lookup test first and confirm borderline results with a sufficiently large parametric bootstrap.

Interpretation and limitations

A small p-value indicates that the observed EDF discrepancy is unusually large under the fitted SN or ST model. It is evidence against the selected model, not a measure of the practical importance of the discrepancy. Conversely, a large p-value does not prove that the fitted family is correct; it means that the test did not detect a departure at the available sample size and resolution.

The tests assume independent observations from a continuous distribution. Serial dependence, clustering, censoring, rounding, or many ties can alter the null distribution. The procedures also inherit limitations of parameter estimation, especially near symmetry, extreme skewness, or an ST fit whose estimated nu lies near a numerical boundary. Graphical checks such as histograms, fitted densities, and probability plots should be used alongside the formal tests.

References

Azzalini, A. (1985). A class of distributions which includes the normal ones. Scandinavian Journal of Statistics, 12(2), 171–178.

Azzalini, A., and Capitanio, A. (2014). The Skew-Normal and Related Families. Cambridge University Press.

Babu, G. J., and Rao, C. R. (2004). Goodness-of-fit tests when parameters are estimated. Sankhya, 66, 63–74.

Mateu-Figueras, G., Puig, P., and Pewsey, A. (2007). Goodness-of-fit tests for the skew-normal distribution when the parameters are estimated from the data. Communications in Statistics—Theory and Methods, 36(9), 1735–1755. doi:10.1080/03610920601126217.