Kernel Conditional Distribution Estimation with Mixed Data Types

npcdist computes kernel cumulative conditional distribution estimates on $p+q$-variate evaluation data, given a set of training data (both explanatory and dependent) and a bandwidth specification (a condbandwidth object or a bandwidth vector, bandwidth type, and kernel type) using the method of Li and Racine (2008) and Li, Lin, and Racine (2013). The data may be continuous, discrete (unordered and ordered factors), or some combination thereof.

Usage

npcdist(bws, ...)

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


# S3 method for class 'condbandwidth'
npcdist(bws,
        txdat = stop("invoked without training data 'txdat'"),
        tydat = stop("invoked without training data 'tydat'"),
        exdat,
        eydat,
        gradients = FALSE,
        gradient.order = 1L,
        proper = FALSE,
        proper.method = c("isotonic"),
        proper.control = list(),
        ...,
        se = FALSE)

# Default S3 method
npcdist(bws, txdat, tydat, nomad = FALSE, ...)

Arguments

Data, Bandwidth Inputs And Formula Interface

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

bws

a bandwidth specification. This can be set as a condbandwidth object returned from a previous invocation of npcdistbw, or as a $p+q$-vector of bandwidths, with each element $i$ up to $i=q$ corresponding to the bandwidth for column $i$ in tydat, and each element $i$ from $i=q+1$ to $i=p+q$ corresponding to the bandwidth for column $i-q$ in txdat. If specified as a vector, then additional arguments will need to be supplied as necessary 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 npcdistbw was called.

txdat

a $p$-variate data frame of sample realizations of explanatory data (training data). Defaults to the corresponding training data recoverable from the bandwidth object in the current session.

tydat

a $q$-variate data frame of sample realizations of dependent data (training data). Defaults to the corresponding training data recoverable from the bandwidth object in the current session.

Local-Polynomial Degree And Bandwidth Search

This argument controls the recommended automatic local-polynomial NOMAD route, which jointly selects continuous polynomial degree and bandwidths when these are computed inside npcdist.

nomad: logical shortcut passed through to npcdistbw when bandwidths are computed inside npcdist. When TRUE, the bandwidth route fills any missing values among regtype, search.engine, degree.select, bernstein.basis, degree.min, degree.max, degree.verify, and bwtype with the recommended automatic local-polynomial degree-and-bandwidth NOMAD preset documented in npcdistbw. Additional NOMAD tuning arguments such as nomad.nmulti may also be supplied through ...; nmulti remains the outer restart count while nomad.nmulti controls inner native crs NOMAD API multistarts within each outer restart. After fitting, inspect fit$bws$nomad.shortcut on the returned object fit to see the normalized shortcut metadata.

Evaluation Data And Returned Quantities

These arguments control where the fitted conditional distribution is evaluated and which estimates are returned.

exdat: a $p$-variate data frame of explanatory data on which cumulative conditional distributions will be evaluated. By default, evaluation takes place on the data provided by txdat.
eydat: a $q$-variate data frame of dependent data on which cumulative conditional distributions will be evaluated. By default, evaluation takes place on the data provided by tydat.
se: logical scalar, default FALSE. If TRUE, compute and retain available asymptotic standard errors of fitted values and, when gradients = TRUE, of the requested gradients. If FALSE, skip uncertainty-only computation; point estimates and gradients are unchanged. Unavailable derivative errors remain NA when requested. To obtain errors without repeating bandwidth selection, refit using npcdist(bws = fit$bws, se = TRUE); also supply gradients = TRUE for gradient errors. Explicit training data may be supplied if they are no longer recoverable from the bandwidth object. The argument must be named.
gradients: a logical value specifying whether to return estimates of the gradients at the evaluation points. Defaults to FALSE.
gradient.order: derivative order for continuous explanatory-variable gradients when gradients = TRUE and regtype = "lp". A scalar is recycled across continuous explanatory variables, or one value may be supplied per continuous explanatory variable. The default 1L returns first derivatives. Higher orders are available for continuous explanatory variables only and must not exceed the corresponding local-polynomial degree, apart from the established all-degree-zero first-derivative compatibility route. An unavailable continuous component is returned as NA; plotting retains it as an empty panel when another continuous derivative or categorical effect is available. A request for which no component is available is rejected. Unordered and ordered predictors retain their usual first-order discrete effects.
newdata: An optional data frame in which to look for evaluation data. If omitted, the training data are used.

Fit Properization Controls

These arguments control optional post-estimation properization of the fitted conditional distribution.

proper: a logical value specifying whether to apply post-estimation properization to the conditional distribution estimate. Defaults to FALSE. Local-constant or degree-zero fits bypass correction only when their kernels and bandwidth geometry guarantee a convex mixture of proper response CDFs. Signed kernels and generalized-nearest-neighbour response bandwidths do not have this guarantee. Requested correction then uses the existing supported grid or slice procedure, with its additional evaluation cost and support requirements.
proper.control: a list of controls for properization. Supported entries are tol, grid.check, store.raw, and fail.on.unsupported. Resolved controls are retained in fitted objects and inherited by predict() and plot() unless explicitly replaced. An explicit control list replaces the retained list; unspecified entries then use their defaults. Prediction still errors on unsupported correction unless fail.on.unsupported=FALSE is explicitly supplied.
proper.method: a character string specifying the properization method. Currently "isotonic" is supported.

Additional Arguments

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

...: additional arguments supplied to npcdistbw when npcdist computes bandwidths internally, or arguments needed to interpret a numeric bws vector. This is where bandwidth selection controls such as bwmethod, bwtype, kernel/support controls such as cxkertype, cykertype, cxkerbound, and cykerbound, search controls such as nmulti, scale.factor.search.lower, and nomad.nmulti, and local-polynomial controls such as regtype, degree, basis, and bernstein.basis are supplied. See npcdistbw for the complete bandwidth-selection argument surface.

Details

Conditional fitting and replay use the literal categorical kernels for both explanatory and response variables. Near-upper categorical bandwidths are not replaced by approximate constant factors. The same target is used for levels, gradients, standard errors, conditional hats, and reconstructed plot uncertainty. This conditional-family policy does not change direct regression or unconditional kernel-sum evaluation and does not modify the caller's np.largelambda option.

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.

With a formula-derived bandwidth object, supplying only tydat replaces the evaluated response in retained-training-sample order: one row per observation retained after the original subset and missing-value selection. Original excluded rows are not restored; row names or equal lengths do not imply original-data order. The effective na.action is then applied jointly to the retained predictors and replacement response; any new omission indices refer to this retained sample. The selected bandwidths remain fixed, and the returned object retains the replacement sample for subsequent refits. Use data for full formula-data replacement, not together with a partial response replacement. Explicit native training arguments retain their separate positional contract.

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 is split into exdat and eydat using the retained predictor and response column names. Names must identify both roles unambiguously; explicitly supplied native evaluation arguments take precedence. This evaluates at the requested rows without repeating bandwidth selection.

For explicit external evaluation rows, a completed explanatory-kernel row whose weights are all zero has no conditional estimate. Such rows return NA for the fitted value and requested uncertainty or derivatives, with one notice. A missing categorical counterfactual affects only its contrast; it does not invalidate an otherwise supported base estimate. A zero response-kernel contribution is not a missing explanatory row. Supported rows retain their usual computation. Required internal evaluations remain strict. With proper = TRUE, complete unsupported fixed-X slices retain NA; projection still applies to supported slices.

Documentation guide: see npcdistbw 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 npcdistbw first and pass bandwidth-selection arguments from ... to that call. When bws is already a condbandwidth object, npcdist estimates with the stored bandwidth metadata in that object.

Beta associated kernels may be selected independently for explanatory X coordinates (cxkertype="beta") and response Y coordinates (cykertype="beta"). Orders 2, 4, 6, and 8, fixed/generalized-NN/ adaptive-NN bandwidths, explicit finite bounds, sample-range bounds, and automatic cv.ls bandwidth selection are supported. Beta kernels act on continuous coordinates alongside the kernels for categorical coordinates. Local-constant, local-linear and local-polynomial fits are supported. For local-constant fits, analytic first derivatives with respect to continuous X coordinates are available with gradients=TRUE; the Y-side CDF factor itself is not differentiated. Fixed/adaptive order-2 Y CDF factors are proper on their declared support; generalized-NN order-2 curves may be locally nonmonotone, and signed orders 4, 6, and 8 may also leave $[0,1]$. No clipping or monotone repair is applied. At exact beta-X endpoints, numerator and denominator jumps are cancelled before classifying the conditional-distribution gradient. A constant Y CDF factor (in particular at its exact support endpoints) gives a structural zero gradient. Genuine ratio jumps are reported as signed infinities and indeterminate numerical cancellations as NA, with warnings. Returned infinities retain that extended-real meaning. Plot methods convert them to missing geometry only in a temporary drawing copy, leaving fitted and plot-data objects unchanged. For ordinary conditional-distribution level plots, bootstrap intervals support different continuous X- and Y-side kernel families or orders, including beta on either side. Beta-containing bootstrap levels use the same signed log-domain conditional ratio as the fitted estimator. With non-fixed bandwidths, np_boot_control(nonfixed="exact") reconstructs each duplicate-row resample and recomputes its nearest-neighbour bandwidths, while nonfixed="frozen" keeps the original fitted bandwidths. In both modes, resolved beta "range" support bounds are replayed from the bandwidth object rather than recomputed from bootstrap counts. This mixed-family level-bootstrap support does not alter the separate restrictions on gradient, quantile, or bias-corrected centers.

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

Argument groups for bandwidth selection are documented on npcdistbw. The most common workflow is to choose data and bandwidth inputs first, then bandwidth criterion and representation, then kernel/support controls, numerical search controls, and finally local-polynomial/NOMAD controls for polynomial-adaptive fits.

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

npcdist implements a variety of methods for estimating multivariate conditional cumulative distributions ($p+q$-variate) defined over a set of possibly continuous and/or discrete (unordered, ordered) data. The approach is based on Li and Racine (2004) 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 conditional distribution at the point $x$. Generalized nearest-neighbor bandwidths change with the point at which the cumulative conditional distribution is estimated, $x$. Fixed bandwidths are constant over the support of $x$.

Training and evaluation input data may be a mix of continuous (default), unordered discrete (to be specified in the data frames using factor), and ordered discrete (to be specified in the data frames using ordered). 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. Unordered discrete data types use a variation on Aitchison and Aitken's (1976) kernel, while ordered data types use a variation of the Wang and van Ryzin (1981) kernel.

For practitioners who want the recommended automatic local-polynomial degree-and-bandwidth NOMAD route without spelling out all LP tuning arguments, npcdist(..., nomad=TRUE) and npcdistbw(..., nomad=TRUE) expand missing settings to the same documented preset. Explicit incompatible settings fail fast rather than being silently rewritten.

With a raw polynomial basis, high powers of poorly scaled continuous predictors (such as uncentred calendar years) can produce severely ill-conditioned systems. Conditional local-polynomial fits do not apply the design-conditioning rejection used by npreg; they may return inaccurate finite values without a warning, with or without a ridge adjustment. Centre or scale such predictors before fitting and apply the same transformation to all evaluation or prediction data. The explicit bernstein.basis = TRUE representation is another option for improving conditioning. Check numerical stability for the chosen data, degree and bandwidths. proper = TRUE enforces probability-shape constraints; it does not repair an ill-conditioned polynomial solve.

With regtype = "lp", gradients = TRUE and gradient.order greater than one expose higher-order derivative estimates with respect to continuous explanatory variables. These derivatives use the same local-polynomial basis, degree, Bernstein option, bandwidths, and kernels as the fitted conditional distribution. See npreg for the complete heterogeneous GLP multi-index and Bernstein-representation definitions shared by these routes. With se = TRUE, available higher-order derivative columns for fixed and adaptive nearest-neighbour bandwidths also have standard errors. These combine the existing local pseudo-response variance with the squared norm of the actual derivative weighting row, including its accepted ridge adjustment. They condition on the selected bandwidths, neighbour radii and local-polynomial map; they do not account for their selection uncertainty. Generalized nearest-neighbour higher-order point derivatives remain subject to the existing hat-operator restriction. Bootstrap inference remains recommended.

With positive-degree local-polynomial fits (including local-linear fits), numerical ridging follows the policies described in npreg (Seifert and Gasser, 2000). The conditional pseudo-response uncertainty formulas here are distinct from ordinary regression's residual normalization.

With positive-degree local-polynomial fits (including local-linear fits), categorical first-difference standard errors use paired endpoint influence calculations. The covariance between the two category endpoints is retained, not replaced by the sum of their marginal variances. These calculations condition on the selected basis, smoothing parameters and realized neighbour radii for all three bandwidth types. They run only when both gradients = TRUE and se = TRUE request them. A nonsmooth regularization boundary or nonfinite influence gives NA uncertainty with a notice, without changing the point contrast.

For ordinary scalar fits, leading interior standard errors use the squared X-kernel moment; conditional density errors additionally use the independent Y-kernel moment. Analytic first-derivative errors for fixed and generalized nearest-neighbour bandwidths with Gaussian or Epanechnikov kernels use the squared derivative-kernel moment and the X bandwidth. These formulas do not establish uniform-kernel derivative inference or full adaptive-NN inference with estimated random radii. Positive-degree local-polynomial fits use their separate covariance calculation.

For fixed and generalized nearest-neighbour scalar fits with categorical variables, these leading moments retain squared categorical kernel factors and the joint training-category frequencies. Conditional distributions, and conditional densities with no continuous response, retain the numerator–denominator covariance of the ratio estimator. In particular, a constant response-kernel contribution has zero conditional variance. Purely categorical fits use this ratio covariance without a continuous bandwidth factor, including nearest-neighbour-labelled bandwidth objects. Beta-kernel and positive-degree local-polynomial covariance calculations retain their separate targets.

For scalar adaptive nearest-neighbour fits with any continuous variable, requested level and available continuous first-derivative standard errors use the complete empirical influence of the conditional ratio. They retain the actual normalized X weights and continuous/categorical Y contributions, including numerator and denominator derivatives. This working covariance holds the selected model, observed X geometry and realized response nearest-neighbour radii fixed; it does not account for their sampling or selection uncertainty. Local-constant-equivalent degree-zero fits share this calculation, including continuous responses with categorical-only X. Positive-degree local-polynomial, beta-kernel and paired categorical- contrast errors retain their separate covariance calculations. With se = FALSE, the additional uncertainty pass is not performed.

Value

The se component records whether uncertainty was requested. conderr and congerr are NULL when se = FALSE. se() and gradients(..., se = TRUE) only extract stored results and give a no-search refit message when the errors were omitted. predict(..., se.fit = TRUE) explicitly requests the needed errors.

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

xbw: bandwidth(s), scale factor(s) or nearest neighbours for the explanatory data, txdat
ybw: bandwidth(s), scale factor(s) or nearest neighbours for the dependent data, tydat
xeval: the evaluation points of the explanatory data
yeval: the evaluation points of the dependent data
condist: estimates of the conditional cumulative distribution at the evaluation points. An evaluation row whose required explanatory normalization is exactly zero or non-finite is returned as NA; other rows remain available.
conderr: when se = TRUE, standard errors of the cumulative conditional distribution estimates. Undefined normalization rows are returned as NA.
congrad: if invoked with gradients = TRUE, estimates of the gradients at the evaluation points. For local-polynomial fits, continuous-coordinate derivative orders are controlled by gradient.order. Undefined normalization rows are returned as NA.
congerr: if invoked with gradients = TRUE and se = TRUE, standard errors of the gradients at the evaluation points, including available higher-order continuous derivatives under the working-covariance assumptions described above. Undefined normalization or uncertainty rows are returned as NA.
log_likelihood: log likelihood of the cumulative conditional distribution estimate

Book And Method Pointers

The conditional distribution target is $F(y\mid x)=\Pr(Y\le y\mid X=x)$. The estimator uses the selected mixed-data conditional distribution bandwidths and kernels for the response and conditioning coordinates; local-polynomial routes use the selected continuous-coordinate polynomial degree. These fitted conditional CDFs are also the object inverted by npqreg when computing conditional quantiles.

Scalar conditional fits involving beta kernels use the empirical observation-influence covariance. For normalized explanatory weights $a_i$, response-kernel contributions $z_i$, and fitted value $m$, the level variance is $\{n/(n-1)\}\sum_i [a_i(z_i-m)]^2$, using the full training count $n$. Supported continuous-derivative errors use the corresponding differentiated influence. Categorical contrasts use the difference of the two endpoints' observation influences inside the same squared norm, thereby retaining their covariance; self-contrasts have zero variance. These calculations hold the selected bandwidth state fixed; they do not add bandwidth-selection or random-NN-radius uncertainty, and do not replace positive-degree LP covariance.

For book-length derivations, see Li and Racine (2007), Chapter 6 Conditional CDF and Quantile Estimation, especially Sections 6.1, 6.2, and 6.5, together with Chapter 5 Conditional Density Estimation. The later workflow treatment is Racine (2019), Chapter 4 Conditional 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.

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. doi:10.22059/jsciences.2024.367933.1007836.

Hall, P. and J.S. Racine and Q. Li (2004), “Cross-validation and the estimation of conditional probability densities,” Journal of the American Statistical Association, 99, 1015-1026.

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

Li, Q. and J.S. Racine (2008), “Nonparametric estimation of conditional CDF and quantile functions with mixed categorical and continuous data,” Journal of Business and Economic Statistics, 26, 423-434.

Li, Q. and J. Lin and J.S. Racine (2013), “Optimal bandwidth selection for nonparametric conditional distribution and quantile functions”, Journal of Business and Economic Statistics, 31, 57-65.

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.

Examples

if (FALSE) { # \dontrun{
# EXAMPLE 1 (INTERFACE=FORMULA): For this example, we load Giovanni
# Baiocchi's Italian GDP panel (see Italy for details), and compute the
# cross-validated bandwidths (default) using a second-order Gaussian
# kernel (default). Note - this may take a minute or two depending on
# the speed of your computer.

data("Italy")
Italy <- Italy[seq_len(min(300, nrow(Italy))), ]
with(Italy, {

# First, compute the bandwidths.

bw <- npcdistbw(formula=gdp~ordered(year), nmulti=1)

# Next, compute the condistribution object...

Fhat <- npcdist(bws=bw)

# The object Fhat now contains results such as the estimated cumulative
# conditional distribution function (Fhat$condist) and so on...

summary(Fhat)

# Call the plot() function to visualize the results (<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), and compute the
# cross-validated bandwidths (default) using a second-order Gaussian
# kernel (default). Note - this may take a minute or two depending on
# the speed of your computer.

data("Italy")
Italy <- Italy[seq_len(min(300, nrow(Italy))), ]
with(Italy, {

# First, compute the bandwidths.

# Note - we cast `X' and `y' as data frames so that plot() can
# automatically grab names (this looks like overkill, but in
# multivariate settings you would do this anyway, so may as well get in
# the habit).

X <- data.frame(year=ordered(year))
y <- data.frame(gdp)

bw <- npcdistbw(xdat=X, ydat=y, nmulti=1)

# Next, compute the condistribution object...

Fhat <- npcdist(bws=bw)

# The object Fhat now contains results such as the estimated cumulative
# conditional distribution function (Fhat$condist) and so on...

summary(Fhat)

# Call the plot() function to visualize the results (<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 from the R `datasets' library and compute the
# conditional distribution function.

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

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

bw <- npcdistbw(formula=eruptions~waiting, nmulti=1)

summary(bw)

# Plot the conditional cumulative distribution function (<ctrl>-C will
# interrupt on *NIX systems, <esc> will interrupt on MS Windows
# systems).

if (interactive()) plot(bw)

})

# EXAMPLE 2 (INTERFACE=DATA FRAME): For this example, we load the old
# faithful geyser data from the R `datasets' library and compute the
# cumulative conditional distribution function.

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

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

# Note - we cast `X' and `y' as data frames so that plot() can
# automatically grab names (this looks like overkill, but in
# multivariate settings you would do this anyway, so may as well get in
# the habit).

X <- data.frame(waiting)
y <- data.frame(eruptions)

bw <- npcdistbw(xdat=X, ydat=y, nmulti=1)

summary(bw)

# Plot the conditional cumulative distribution function (<ctrl>-C will
# interrupt on *NIX systems, <esc> will interrupt on MS Windows systems)

if (interactive()) plot(bw)

})

# EXAMPLE 3: Variations on local polynomial conditional distribution
# estimation with proper = TRUE.

data("Italy")

# Work in years relative to 1975 to condition the raw polynomial basis.
# Calendar year = model year + year.offset, also for prediction data.
year.offset <- 1975
Italy2 <- within(Italy, {
  year <- as.numeric(as.character(year)) - year.offset
})

# Plot only: make the plotted surface proper on the plot evaluation grid.

Fhat <- npcdist(gdp ~ year, data = Italy2,
                regtype = "lp", degree = 3, nmulti = 1)

plot(Fhat, proper = TRUE)

# Fit an object whose fitted values are themselves proper.

ctrl_fit <- list(
  mode = "slice",
  apply = "fitted",
  slice.grid.size = 101L,
  slice.extend.factor = 0.1
)

Fhat_fit <- npcdist(
  gdp ~ year,
  data = Italy2,
  regtype = "lp",
  degree = 3,
  nmulti = 1,
  proper = TRUE,
  proper.control = ctrl_fit
)

fit_proper <- fitted(Fhat_fit)
fit_raw <- Fhat_fit$condist.raw

# Predict on a common explicit y-grid for several years, and render
# those predictions proper.

g.grid <- seq(min(Italy2$gdp), max(Italy2$gdp), length.out = 200)

nd_grid <- expand.grid(
  gdp = g.grid,
  year = c(1955, 1975, 1995) - year.offset
)

pred_grid <- predict(Fhat, newdata = nd_grid, proper = TRUE)

# Predict on paired rows with different gdp grids by year, and still
# make the predictions proper via slice mode.

g1 <- seq(quantile(Italy2$gdp, 0.10),
          quantile(Italy2$gdp, 0.60), length.out = 60)
g2 <- seq(quantile(Italy2$gdp, 0.30),
          quantile(Italy2$gdp, 0.90), length.out = 35)

nd_slice <- rbind(
  data.frame(gdp = g1, year = rep(1960 - year.offset, length(g1))),
  data.frame(gdp = g2, year = rep(1985 - year.offset, length(g2)))
)

pred_slice <- predict(
  Fhat,
  newdata = nd_slice,
  proper = TRUE,
  proper.control = list(mode = "slice")
)

# One object that carries properization for fitted values and for later
# predict() calls.

ctrl_both <- list(
  mode = "slice",
  apply = "both",
  slice.grid.size = 101L,
  slice.extend.factor = 0.1
)

Fhat_both <- npcdist(
  gdp ~ year,
  data = Italy2,
  regtype = "lp",
  degree = 3,
  nmulti = 1,
  proper = TRUE,
  proper.control = ctrl_both
)

fit_both <- fitted(Fhat_both)
pred_both <- predict(
  Fhat_both,
  newdata = nd_slice,
  proper.control = ctrl_both
)

} # }