Kernel Distribution Estimation with Mixed Data Types

npudist computes kernel unconditional cumulative distribution estimates on evaluation data, given a set of training data and a bandwidth specification (a dbandwidth object or a bandwidth vector, bandwidth type, and kernel type) using the method of Li, Li and Racine (2017).

Usage

npudist(bws, ...)

# S3 method for class 'formula'
npudist(bws, data = NULL, newdata = NULL, ...)

# S3 method for class 'dbandwidth'
npudist(bws,
        tdat = stop("invoked without training data 'tdat'"),
        edat,
        ...,
        se = FALSE)


# Default S3 method
npudist(bws, tdat, ..., se = FALSE)

Arguments

se

Logical; compute and retain asymptotic standard errors when TRUE. The default FALSE skips uncertainty-only computation without changing the point estimates.

Data, Bandwidth Inputs And Formula Interface

These arguments identify the bandwidth specification, formula/data interface, and training data.

bws

a dbandwidth specification. This can be set as a dbandwidth object returned from an invocation of npudistbw, or as a \(p\)-vector of bandwidths, with an element for each variable in the training data. If specified as a vector, then additional arguments will need to be supplied as necessary to change them from the defaults to specify the bandwidth type, kernel types, training data, and so on.

data

This argument belongs to the formula interface (including supported formula-bearing bandwidth-object refits). Native calls must supply their vectors, matrices or data frames directly, or use with(data, ...). Supplying data, including explicit data = NULL, to a native route that does not consume it is an error; see np-data-interfaces.

an optional data frame, list or environment (or object coercible to a data frame by as.data.frame) containing the variables in the model. If not found in data, the variables are taken from environment(bws), typically the environment from which npudistbw was called.

tdat

a \(p\)-variate data frame of sample realizations (training data) used to estimate the cumulative distribution. Defaults to the corresponding training data recoverable from the bandwidth object in the current session.

Evaluation Data

These arguments control where the fitted cumulative distribution is evaluated.

edat

a \(p\)-variate data frame of cumulative distribution evaluation points. By default, evaluation takes place on the data provided by tdat.

newdata

An optional data frame in which to look for evaluation data. If omitted, the training data are used.

Additional Arguments

Further arguments are passed to npudistbw when bandwidths are computed internally, or used to interpret a numeric bws vector.

...

additional arguments supplied to npudistbw when npudist computes bandwidths internally, or arguments needed to interpret a numeric bws vector. This is where bandwidth selection controls such as bwmethod, bwtype, bwscaling, kernel/support controls such as ckertype, ckerorder, and ckerbound, categorical kernel controls such as okertype, and search controls such as nmulti, ngrid, and scale.factor.search.lower are supplied. See npudistbw for the complete bandwidth-selection argument surface.

Details

With generalized nearest-neighbor bandwidths, query identity matters when selecting a radius: a training observation excludes its own occurrence, whereas explicit evaluation data remain external queries even if their values equal training observations. Consequently, predict(fit, newdata = training.data) can differ from fitted(fit). This does not omit the training contribution from the fitted estimator. See np.kernels for the radius convention.

For native (non-formula) objects, named newdata columns are matched to the retained variable names and reordered as needed. Extra columns are ignored; missing, duplicate or ambiguous names are errors. An unnamed matrix or vector remains positional for a single evaluation role. Explicit native evaluation arguments retain their positional contract and take precedence over newdata.

For a native bandwidth object, newdata supplies edat; explicitly supplied native evaluation arguments take precedence. This evaluates at the requested rows without repeating bandwidth selection.

With bwtype="adaptive_nn" and se=TRUE, standard errors include the same-sample influence of each donor-centred nearest-neighbour radius. A two-sided rank-spacing pilot estimates the marginal density sum at the radius endpoints. Cross-coordinate and categorical-kernel covariance are retained in the joint observation influence. The selected NN ranks are treated as fixed; bandwidth-selection uncertainty is not included. This calculation supports Gaussian and Epanechnikov kernels of orders 2, 4, 6 and 8, and uniform kernels, without finite kernel bounds. Uniform density errors additionally include moving-support boundary influences with joint other-coordinate and categorical kernel marks. These use type-7 sample-quantile endpoints and deterministic two-sided rank-spacing pilots whose interior window grows as the two-thirds power of sample size. A strictly absent support face contributes zero; an exact probability endpoint or undefined finite-face pilot is unavailable. Finite kernel bounds (including beta kernels), extended or saturated NN ranks, and degenerate pilot spacings are not currently qualified: point estimates are retained and standard errors are NA, with one availability notice. No additional uncertainty computation is performed for se=FALSE. Requested adaptive-NN standard errors can take several times as long as the point fit because they include the donor-radius influence calculation; this cost is incurred only when standard errors are requested.

For ordinary kernels with mixed continuous and ordered categorical data, fixed and generalized-nearest-neighbour standard errors retain the leading continuous-indicator approximation and the squared integrated categorical kernel product. The continuous integrated kernel remains to power one in that second moment; this is not a finite-smoothing sample-contribution variance. Purely continuous and purely categorical conventions are unchanged. Generalized-NN errors condition on the current query radii. The separate adaptive-NN influence calculation is described above.

Documentation guide: see npudistbw for bandwidth selection and search controls, np.kernels for kernels, np.options for global options, and plot, plot.np for plotting options.

When bws is omitted, the formula and default methods call npudistbw first and pass bandwidth-selection arguments from ... to that call. When bws is already a dbandwidth object, npudist estimates with the stored bandwidth metadata in that object.

With finite continuous bounds, npudist uses the observation-centered truncated-kernel CDF documented in np.kernels. It is exactly zero and one at the declared lower and upper endpoints and has the linear uniform-CDF limit as the bandwidth diverges.

The beta associated CDF is available at orders 2, 4, 6, and 8 with manual or cross-validated metric bandwidths. Fixed, generalized-NN, and adaptive-NN modes are supported. Set ckertype="beta" and use either ckerbound="fixed" with finite ckerlb/ckerub, or ckerbound="range" for beta outer half-spacing sample bounds. The data must be continuous and all training and evaluation values must lie within the resulting supports. Categorical products and normal-reference selection are not supported.

For observation \(x_i\), target \(x\), support \([a,b]\), \(L=b-a\), \(z_i=(x_i-a)/L\), \(t=(x-a)/L\), and \(\tau=(L/h)^2\), the contribution is $$I_t(1+z_i\tau,1+(1-z_i)\tau).$$ The returned estimate is the sample mean of these contributions; its reported standard error is the sample standard deviation of the multivariate contribution divided by \(\sqrt{n}\). The implementation accumulates this variance in centered form to avoid cancellation near zero or one. This beta CDF is observation-centered and is not the target-variable integral of Chen's evaluation-centered beta density kernel. See np.kernels for the design rationale and references.

Fixed and adaptive order-2 beta CDF estimates are proper CDFs on the declared support. A generalized-NN order-2 curve retains exact endpoints and contribution-wise range but may be locally nonmonotone because its bandwidth changes with the target. Orders 4, 6, and 8 are signed and may overshoot \([0,1]\) or be locally nonmonotone in any bandwidth mode. Returned values are not silently clipped or monotonized.

For purely ordered data, let \(Z_i(x)\) be the complete multivariate CDF contribution and \(\bar Z(x)=n^{-1}\sum_i Z_i(x)\). Standard errors use the empirical-population variance of these contributions divided by sample size, $$\widehat{\mathrm{se}}\{\hat F(x)\}= \frac{1}{n}\sqrt{\sum_{i=1}^n\{Z_i(x)-\bar Z(x)\}^2}.$$ At zero smoothing this retains the usual \(\sqrt{\hat F(x)(1-\hat F(x))/n}\). At finite smoothing a contribution need not be an indicator, so the Bernoulli expression does not generally apply. Profile compression retains original observation multiplicities.

Argument groups for bandwidth selection are documented on npudistbw. The most common workflow is to choose data and bandwidth inputs first, then bandwidth criterion and representation, then kernel/support controls, distribution-specific integral/grid controls, and numerical search controls if defaults need to be changed.

For S3 plotting help, see plot.np. You can list available plot methods with methods("plot").

Typical usages are (see below for a complete list of options and also the examples at the end of this help file)

    Usage 1: first compute the bandwidth object via npudistbw and then
    compute the cumulative distribution:

    bw <- npudistbw(~y)
    Fhat <- npudist(bw)

    Usage 2: alternatively, compute the bandwidth object indirectly:

    Fhat <- npudist(~y)

    Usage 3: modify the default kernel and order:

    Fhat <- npudist(~y, ckertype="epanechnikov", ckerorder=4)

    Usage 4: use the data frame interface rather than the formula
    interface:

    Fhat <- npudist(tdat = y, ckertype="epanechnikov", ckerorder=4)
  

npudist implements a variety of methods for estimating multivariate cumulative distributions (\(p\)-variate) defined over a set of possibly continuous and/or discrete (ordered) data. The approach is based on Li and Racine (2003) who employ ‘generalized product kernels’ that admit a mix of continuous and discrete data types.

Three classes of kernel estimators for the continuous data types are available: fixed, adaptive nearest-neighbor, and generalized nearest-neighbor. Adaptive nearest-neighbor bandwidths change with each sample realization in the set, \(x_i\), when estimating the cumulative distribution at the point \(x\). Generalized nearest-neighbor bandwidths change with the point at which the cumulative distribution is estimated, \(x\). Fixed bandwidths are constant over the support of \(x\).

Data contained in the data frame tdat (and also edat) may be a mix of continuous (default) and ordered discrete (to be specified in the data frame tdat using the ordered command). Data can be entered in an arbitrary order and data types will be detected automatically by the routine (see np for details).

A variety of kernels may be specified by the user. Kernels implemented for continuous data types include the second, fourth, sixth, and eighth-order Gaussian and Epanechnikov kernels, and the uniform kernel. Ordered data types use a variation of the Wang and van Ryzin (1981) kernel.

Value

Standard errors are opt-in. To add them to an existing fit without repeating bandwidth selection, use fit <- npudist(bws=fit$bws, se=TRUE), supplying the original training data and desired evaluation data again when these cannot be recovered from the bandwidth object's call. The se() extractor reads retained errors only and gives a helpful refit message when they were not requested. Use predict(fit, se.fit=TRUE, ...) for prediction with standard errors.

npudist returns a npdistribution object. The generic accessor functions fitted and se extract estimated values and asymptotic standard errors on estimates, respectively, from the returned object. Furthermore, the functions predict, summary and plot support objects of both classes. The returned objects have the following components:

eval

the evaluation points.

dist

estimate of the cumulative distribution at the evaluation points

derr

Standard errors of the cumulative distribution estimates when se=TRUE, otherwise an empty numeric vector.

se

Whether standard-error computation was requested.

Book And Method Pointers

The unconditional distribution target is \(F(x)=\Pr(X\le x)\). The estimator replaces the density kernel factors by integrated distribution-kernel factors where appropriate, while retaining the same mixed-data product-kernel logic for continuous and ordered components. The returned standard errors are asymptotic standard errors for the estimated cumulative distribution values.

For book-length derivations, see Li and Racine (2007), Chapter 1 Density Estimation, especially Sections 1.4 and 1.5 for CDF estimation and bandwidth selection, together with Chapter 4 Kernel Estimation with Mixed Data. The later workflow treatment is Racine (2019), Chapter 2 Continuous Density and Cumulative Distribution Functions and Chapter 3 Mixed-Data Probability Density and Cumulative Distribution Functions.

References

Aitchison, J. and C.G.G. Aitken (1976), “ Multivariate binary discrimination by the kernel method,” Biometrika, 63, 413-420.

Li, Q. and J.S. Racine (2007), Nonparametric Econometrics: Theory and Practice, Princeton University Press.

Li, Q. and J.S. Racine (2003), “Nonparametric estimation of distributions with categorical and continuous data,” Journal of Multivariate Analysis, 86, 266-292.

Li, C. and H. Li and J.S. Racine (2017), “Cross-Validated Mixed Datatype Bandwidth Selection for Nonparametric Cumulative Distribution/Survivor Functions,” Econometric Reviews, 36, 970-987.

Chen, S.X. (1999), “Beta kernel estimators for density functions,” Computational Statistics & Data Analysis, 31(2), 131–145. doi:10.1016/S0167-9473(99)00010-9.

Igarashi, G. and Kakizawa, Y. (2020), “Higher-order bias corrections for asymmetric kernel density estimators,” Journal of Nonparametric Statistics, 32(3), 716–752. doi:10.1080/10485252.2020.1770754.

Mansouri, B., Rastin, A., and Allah Mombeni, H. (2024), “Beta Kernel Estimator for a Cumulative Distribution Function with Bounded Support,” Journal of Sciences, Islamic Republic of Iran, 34(4), 349–361 (volume dated 2023; accepted and DOI registered in 2024). doi:10.22059/jsciences.2024.367933.1007836.

Ouyang, D. and Q. Li and J.S. Racine (2006), “Cross-validation and the estimation of probability distributions with categorical data,” Journal of Nonparametric Statistics, 18, 69-100.

Pagan, A. and A. Ullah (1999), Nonparametric Econometrics, Cambridge University Press.

Scott, D.W. (1992), Multivariate Density Estimation. Theory, Practice and Visualization, New York: Wiley.

Silverman, B.W. (1986), Density Estimation, London: Chapman and Hall.

Wang, M.C. and J. van Ryzin (1981), “A class of smooth estimators for discrete distributions,” Biometrika, 68, 301-309.

Author

Tristen Hayfield tristen.hayfield@gmail.com, Jeffrey S. Racine racinej@mcmaster.ca

Usage Issues

If you are using data of mixed types, then it is advisable to use the data.frame function to construct your input data and not cbind, since cbind will typically not work as intended on mixed data types and will coerce the data to the same type.

For formula fits, character-valued missing-data policies are resolved through the stats model-frame owner; pass a caller-local policy as a function object rather than its name. With replacement training data, an explicit na.action, including NULL, takes precedence over the saved policy. Native fitting still requires complete rows, so a permissive policy does not guarantee the same retained sample as stats::model.frame.

When reusing a bandwidth object without explicit training or evaluation data, newly constructed bandwidth objects use retained training values. A fit with explicit replacement training data retains that sample in its returned bandwidth object; evaluation grids do not replace training data. Older objects without retained samples use their stored call environment when available and may require explicit training inputs after saveRDS() and readRDS(). Custom formula transformations used on new evaluation data must still be available.

See also

np.kernels, np.options, plot, plot.np npudistbw , density

Examples

if (FALSE) { # \dontrun{
# EXAMPLE 1 (INTERFACE=FORMULA): For this example, we load Giovanni
# Baiocchi's Italian GDP panel (see Italy for details), then create a
# data frame in which year is an ordered factor, GDP is continuous,
# compute bandwidths using cross-validation, then create a grid of data
# on which the cumulative distribution will be evaluated for plotting
# purposes.

data("Italy")
with(Italy, {

# Compute bandwidths using cross-validation (default).

bw <- npudistbw(formula=~ordered(year)+gdp)

# At this stage you could use npudist() to do a variety of things. Here
# we compute the npudist() object and place it in Fhat.

Fhat <- npudist(bws=bw)

# Note that simply typing the name of the object returns some useful
# information. For more info, one can call summary:

summary(Fhat)

# Next, we illustrate how to create a grid of `evaluation data' and feed
# it to the perspective plotting routines in R, among others.

# Create an evaluation data matrix

year.seq <- sort(unique(year))
gdp.seq <- seq(1,36,length=50)
data.eval <- expand.grid(year=year.seq,gdp=gdp.seq)

# Generate the estimated cumulative distribution computed for the
# evaluation data

Fhat <- fitted(npudist(bws=bw, newdata=data.eval))

# Coerce the data into a matrix for plotting with persp()

Fhat <- matrix(Fhat, length(unique(year)), 50)

# Next, create a 3D perspective plot of the CDF F, and a 2D
# contour plot.

persp(as.integer(levels(year.seq)), gdp.seq, Fhat, col="lightblue",
      ticktype="detailed", ylab="GDP", xlab="Year", zlab="Density",
      theta=300, phi=50)

# Sleep for 5 seconds so that we can examine the output...

if (interactive()) Sys.sleep(5)

contour(as.integer(levels(year.seq)),
        gdp.seq,
        Fhat,
        xlab="Year",
        ylab="GDP",
        main = "Cumulative Distribution Contour Plot",
        col=topo.colors(100))

# Sleep for 5 seconds so that we can examine the output...

if (interactive()) Sys.sleep(5)

# Alternatively, you could use the plot() command (<ctrl>-C will
# interrupt on *NIX systems, <esc> will interrupt on MS Windows
# systems).

if (interactive()) plot(bw)

})

# EXAMPLE 1 (INTERFACE=DATA FRAME): For this example, we load Giovanni
# Baiocchi's Italian GDP panel (see Italy for details), then create a
# data frame in which year is an ordered factor, GDP is continuous,
# compute bandwidths using cross-validation, then create a grid of data
# on which the cumulative distribution will be evaluated for plotting
# purposes.

data("Italy")
with(Italy, {

data <- data.frame(year=ordered(year), gdp)

# Compute bandwidths using cross-validation (default).

bw <- npudistbw(dat=data)

# At this stage you could use npudist() to do a variety of
# things. Here we compute the npudist() object and place it in Fhat.

Fhat <- npudist(bws=bw)

# Note that simply typing the name of the object returns some useful
# information. For more info, one can call summary:

summary(Fhat)

# Next, we illustrate how to create a grid of `evaluation data' and feed
# it to the perspective plotting routines in R, among others.

# Create an evaluation data matrix

year.seq <- sort(unique(year))
gdp.seq <- seq(1,36,length=50)
data.eval <- expand.grid(year=year.seq,gdp=gdp.seq)

# Generate the estimated cumulative distribution computed for the
# evaluation data

Fhat <- fitted(npudist(edat = data.eval, bws=bw))

# Coerce the data into a matrix for plotting with persp()

Fhat <- matrix(Fhat, length(unique(year)), 50)

# Next, create a 3D perspective plot of the CDF F, and a 2D
# contour plot.

persp(as.integer(levels(year.seq)), gdp.seq, Fhat, col="lightblue",
      ticktype="detailed", ylab="GDP", xlab="Year",
      zlab="Cumulative Distribution",
      theta=300, phi=50)

# Sleep for 5 seconds so that we can examine the output...

if (interactive()) Sys.sleep(5)

contour(as.integer(levels(year.seq)),
        gdp.seq,
        Fhat,
        xlab="Year",
        ylab="GDP",
        main = "Cumulative Distribution Contour Plot",
        col=topo.colors(100))

# Sleep for 5 seconds so that we can examine the output...

if (interactive()) Sys.sleep(5)

# Alternatively, you could use the plot() command (<ctrl>-C will
# interrupt on *NIX systems, <esc> will interrupt on MS Windows
# systems).

if (interactive()) plot(bw)

})

# EXAMPLE 2 (INTERFACE=FORMULA): For this example, we load the old
# faithful geyser data and compute the cumulative distribution function.

library("datasets")
data("faithful")
with(faithful, {

# Note - this may take a few minutes depending on the speed of your
# computer...

bw <- npudistbw(formula=~eruptions+waiting)

summary(bw)

# Plot the cumulative distribution function (<ctrl>-C will interrupt on
# *NIX systems, <esc> will interrupt on MS Windows systems). Note that
# we use xtrim = -0.2 to extend the plot outside the support of the data
# (i.e., extend the tails of the estimate to meet the horizontal axis).

if (interactive()) plot(bw, xtrim=-0.2)

})

# EXAMPLE 2 (INTERFACE=DATA FRAME): For this example, we load the old
# faithful geyser data and compute the cumulative distribution function.

library("datasets")
data("faithful")
with(faithful, {

# Note - this may take a few minutes depending on the speed of your
# computer...

bw <- npudistbw(dat=faithful)

summary(bw)

# Plot the cumulative distribution function (<ctrl>-C will interrupt on
# *NIX systems, <esc> will interrupt on MS Windows systems). Note that
# we use xtrim = -0.2 to extend the plot outside the support of the data
# (i.e., extend the tails of the estimate to meet the horizontal axis).

if (interactive()) plot(bw, xtrim=-0.2)

})
} # }