---
title: "Introduction to PBGoF"
author: "Hongxiang Li, Chenglin Xu, and Tsung Fei Khang"
date: "2026-09-16"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Introduction to PBGoF}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

<style>
#MathJax_Message { display: none !important; }
h1.title { hyphens: none; overflow-wrap: normal; word-break: normal; }
</style>

# 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

## Skew-normal model

Under the direct parameterization (DP), the skew-normal model has location
xi, positive scale omega, and shape alpha. Setting alpha to zero gives a normal
distribution; positive and negative values produce right- and left-skewed
densities, respectively.

The centered parameterization (CP) expresses the same model through:

- mean: the population mean;
- sd: the population standard deviation;
- gamma1: the population coefficient of skewness.

DP is convenient for evaluating the fitted SN distribution function. CP is
convenient for matching fitted skewness to the precomputed simulation tables.
The signs of alpha and gamma1 agree: positive alpha corresponds to positive
gamma1, and negative alpha corresponds to negative gamma1.

## Skew-t model

The skew-t model can represent asymmetry and heavy tails simultaneously. Its
direct parameterization is $(\xi,\omega,\alpha,\nu)$, where $\xi$ is location,
$\omega>0$ is scale, $\alpha$ controls asymmetry, and $\nu>0$ is the degrees of
freedom controlling tail weight. As $\nu$ becomes large, the skew-t model
approaches the skew-normal model.

PBGoF fits the model with `st.fit.robust()`, evaluates it graphically with
`st.plot.check()`, and provides parameter-bootstrap KS and CvM tests through
`st.para.bootstrap.ks.test()` and `st.para.bootstrap.cvm.test()`. The bootstrap
procedure re-estimates all four parameters in every valid replicate.

## Why lookup uses the absolute skewness

Suppose X follows an SN distribution with DP parameters (xi, omega, alpha).
The reflected variable -X is also skew-normal, with shape -alpha. Reflection
reverses the sign of the CP skewness gamma1 but preserves the magnitude and
does not change the sampling distribution of a reflection-invariant EDF
goodness-of-fit statistic.

Mateu-Figueras, Puig, and Pewsey (2007) explicitly report that the
distributions of the five EDF statistics they study are invariant to changes
in the sign of the skew-normal shape parameter. Their test instructions state
that when the fitted shape parameter is negative, the table corresponding to
its positive counterpart should be used.

PBGoF indexes its simulation tables by the CP coefficient gamma1 rather than
the DP shape alpha. Because gamma1 changes sign together with alpha, PBGoF
implements the same symmetry through

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

followed by restriction to the available table range 0.01--0.99. Thus, for
example, fitted values gamma1_hat = -0.63 and gamma1_hat = 0.63 both use the
gamma1 = 0.63 row. The original sign is retained in gamma1_hat in the returned
result, while gamma1_used records the non-negative lookup value.

# 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.