nw_predict etc., another v. bump

Author: soodoku

hbw/__init__.py

@@ -10,20 +10,33 @@
 """
 
 from ._kernels import _KERNELS
-from .kde import kde_bandwidth, kde_bandwidth_mv, lscv, lscv_grad, lscv_mv, lscv_score
+from .kde import (
+    kde_bandwidth,
+    kde_bandwidth_mv,
+    kde_evaluate,
+    kde_evaluate_mv,
+    lscv,
+    lscv_grad,
+    lscv_mv,
+    lscv_score,
+)
 from .nw import (
     loocv_mse,
     loocv_mse_grad,
     loocv_mse_mv,
     loocv_mse_score,
     nw_bandwidth,
     nw_bandwidth_mv,
+    nw_predict,
+    nw_predict_mv,
 )
 
 __all__ = [
     "_KERNELS",
     "kde_bandwidth",
     "kde_bandwidth_mv",
+    "kde_evaluate",
+    "kde_evaluate_mv",
     "loocv_mse",
     "loocv_mse_grad",
     "loocv_mse_mv",
@@ -34,5 +47,7 @@
     "lscv_score",
     "nw_bandwidth",
     "nw_bandwidth_mv",
+    "nw_predict",
+    "nw_predict_mv",
 ]
-__version__ = "0.5.0"
+__version__ = "0.6.0"

hbw/kde.py

@@ -377,6 +377,132 @@ def _newton_armijo_mv_numba(
     return h
 
 
+def kde_evaluate(
+    x_train: ArrayLike,
+    x_eval: ArrayLike,
+    h: float,
+    kernel: str = "gauss",
+) -> NDArray[Any]:
+    """Evaluate kernel density estimate at given points.
+
+    Parameters
+    ----------
+    x_train
+        Training sample data (1D array-like).
+    x_eval
+        Points at which to evaluate the density (1D array-like).
+    h
+        Bandwidth (obtained from kde_bandwidth).
+    kernel
+        Kernel function: "gauss", "epan", "unif", "biweight", "triweight", or "cosine".
+
+    Returns
+    -------
+    NDArray
+        Estimated density values at x_eval locations.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> x = np.random.randn(1000)
+    >>> h = kde_bandwidth(x)
+    >>> x_grid = np.linspace(-3, 3, 100)
+    >>> density = kde_evaluate(x, x_grid, h)
+
+    Notes
+    -----
+    This function returns point estimates only. No confidence intervals or
+    standard errors are provided. Key assumptions: smooth underlying density,
+    IID observations, continuous density (no point masses).
+
+    For inference, consider bootstrap resampling or see statsmodels.nonparametric.
+    """
+    x_tr = np.asarray(x_train, dtype=float).ravel()
+    x_ev = np.asarray(x_eval, dtype=float).ravel()
+
+    if kernel not in _KERNELS:
+        raise ValueError(f"kernel must be one of {list(_KERNELS.keys())}, got {kernel!r}")
+
+    K, _, _, _, _, _ = _KERNELS[kernel]
+    n = len(x_tr)
+
+    u = (x_ev[:, None] - x_tr[None, :]) / h
+    density = K(u).sum(axis=1) / (n * h)
+
+    return density
+
+
+def kde_evaluate_mv(
+    data_train: ArrayLike,
+    data_eval: ArrayLike,
+    h: float,
+    kernel: str = "gauss",
+) -> NDArray[Any]:
+    """Evaluate multivariate kernel density estimate at given points using product kernel.
+
+    Parameters
+    ----------
+    data_train
+        Training sample data, shape (n_train, d).
+    data_eval
+        Points at which to evaluate the density, shape (n_eval, d).
+    h
+        Bandwidth (scalar, applied to all dimensions; obtained from kde_bandwidth_mv).
+    kernel
+        Kernel function: "gauss", "epan", "unif", "biweight", "triweight", or "cosine".
+
+    Returns
+    -------
+    NDArray
+        Estimated density values at data_eval locations.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> data = np.random.randn(500, 2)
+    >>> h = kde_bandwidth_mv(data)
+    >>> data_grid = np.column_stack([np.linspace(-3, 3, 50), np.linspace(-3, 3, 50)])
+    >>> density = kde_evaluate_mv(data, data_grid, h)
+
+    Notes
+    -----
+    This function returns point estimates only. No confidence intervals or
+    standard errors are provided. Key assumptions: smooth underlying density,
+    IID observations, continuous density (no point masses). Uses product kernel
+    with isotropic bandwidth; data should be standardized for best results.
+
+    For inference, consider bootstrap resampling or see statsmodels.nonparametric.
+    """
+    data_tr = np.asarray(data_train, dtype=float)
+    data_ev = np.asarray(data_eval, dtype=float)
+
+    if data_tr.ndim == 1:
+        data_tr = data_tr.reshape(-1, 1)
+    if data_ev.ndim == 1:
+        data_ev = data_ev.reshape(-1, 1)
+
+    if data_tr.ndim != 2:
+        raise ValueError(f"data_train must be 2D array, got shape {data_tr.shape}")
+    if data_ev.ndim != 2:
+        raise ValueError(f"data_eval must be 2D array, got shape {data_ev.shape}")
+    if data_tr.shape[1] != data_ev.shape[1]:
+        raise ValueError(f"data_train and data_eval must have same number of dimensions, got {data_tr.shape[1]} and {data_ev.shape[1]}")
+    if kernel not in _KERNELS:
+        raise ValueError(f"kernel must be one of {list(_KERNELS.keys())}, got {kernel!r}")
+
+    K, _, _, _, _, _ = _KERNELS[kernel]
+    n = len(data_tr)
+    d = data_tr.shape[1]
+
+    U = (data_ev[:, None, :] - data_tr[None, :, :]) / h
+    K_vals = K(U)
+    K_prod = np.prod(K_vals, axis=2)
+
+    density = K_prod.sum(axis=1) / (n * h**d)
+
+    return density
+
+
 def kde_bandwidth_mv(
     data: ArrayLike,
     kernel: str = "gauss",

hbw/nw.py

@@ -598,6 +598,156 @@ def _newton_armijo_mv_nw_numba(
     return h
 
 
+def nw_predict(
+    x_train: ArrayLike,
+    y_train: ArrayLike,
+    x_test: ArrayLike,
+    h: float,
+    kernel: str = "gauss",
+) -> NDArray[Any]:
+    """Nadaraya-Watson kernel regression predictions.
+
+    Parameters
+    ----------
+    x_train
+        Training predictor values (1D array-like).
+    y_train
+        Training response values (1D array-like).
+    x_test
+        Test predictor values where predictions are desired (1D array-like).
+    h
+        Bandwidth (obtained from nw_bandwidth).
+    kernel
+        Kernel function: "gauss", "epan", "unif", "biweight", "triweight", or "cosine".
+
+    Returns
+    -------
+    NDArray
+        Predicted values at x_test locations.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> x = np.linspace(-2, 2, 200)
+    >>> y = np.sin(x) + 0.1 * np.random.randn(len(x))
+    >>> h = nw_bandwidth(x, y)
+    >>> y_pred = nw_predict(x, y, x, h)  # in-sample predictions
+
+    Notes
+    -----
+    This function returns point estimates only. No confidence intervals or
+    standard errors are provided. Key assumptions: smooth regression function,
+    IID observations, design density bounded away from zero in region of interest.
+
+    For inference, consider bootstrap resampling or see statsmodels.nonparametric.
+    """
+    x_tr = np.asarray(x_train, dtype=float).ravel()
+    y_tr = np.asarray(y_train, dtype=float).ravel()
+    x_te = np.asarray(x_test, dtype=float).ravel()
+
+    if len(x_tr) != len(y_tr):
+        raise ValueError(f"x_train and y_train must have same length, got {len(x_tr)} and {len(y_tr)}")
+    if kernel not in _KERNELS:
+        raise ValueError(f"kernel must be one of {list(_KERNELS.keys())}, got {kernel!r}")
+
+    K, _, _, _, _, _ = _KERNELS[kernel]
+
+    u = (x_te[:, None] - x_tr[None, :]) / h
+    w = K(u)
+
+    w_sum = w.sum(axis=1)
+    w_sum_safe = np.where(w_sum == 0, np.finfo(float).eps, w_sum)
+    y_pred = (w @ y_tr) / w_sum_safe
+
+    zero_weight_mask = w_sum == 0
+    if np.any(zero_weight_mask):
+        y_pred[zero_weight_mask] = np.mean(y_tr)
+
+    return y_pred
+
+
+def nw_predict_mv(
+    data_train: ArrayLike,
+    y_train: ArrayLike,
+    data_test: ArrayLike,
+    h: float,
+    kernel: str = "gauss",
+) -> NDArray[Any]:
+    """Multivariate Nadaraya-Watson kernel regression predictions using product kernel.
+
+    Parameters
+    ----------
+    data_train
+        Training predictor values, shape (n_train, d).
+    y_train
+        Training response values (1D array of length n_train).
+    data_test
+        Test predictor values where predictions are desired, shape (n_test, d).
+    h
+        Bandwidth (scalar, applied to all dimensions; obtained from nw_bandwidth_mv).
+    kernel
+        Kernel function: "gauss", "epan", "unif", "biweight", "triweight", or "cosine".
+
+    Returns
+    -------
+    NDArray
+        Predicted values at data_test locations.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> data = np.random.randn(500, 2)
+    >>> y = np.sin(data[:, 0]) + 0.5 * data[:, 1] + 0.3 * np.random.randn(500)
+    >>> h = nw_bandwidth_mv(data, y)
+    >>> y_pred = nw_predict_mv(data, y, data, h)  # in-sample predictions
+
+    Notes
+    -----
+    This function returns point estimates only. No confidence intervals or
+    standard errors are provided. Key assumptions: smooth regression function,
+    IID observations, design density bounded away from zero in region of interest.
+    Uses product kernel with isotropic bandwidth; data should be standardized
+    for best results.
+
+    For inference, consider bootstrap resampling or see statsmodels.nonparametric.
+    """
+    data_tr = np.asarray(data_train, dtype=float)
+    y_tr = np.asarray(y_train, dtype=float).ravel()
+    data_te = np.asarray(data_test, dtype=float)
+
+    if data_tr.ndim == 1:
+        data_tr = data_tr.reshape(-1, 1)
+    if data_te.ndim == 1:
+        data_te = data_te.reshape(-1, 1)
+
+    if data_tr.ndim != 2:
+        raise ValueError(f"data_train must be 2D array, got shape {data_tr.shape}")
+    if data_te.ndim != 2:
+        raise ValueError(f"data_test must be 2D array, got shape {data_te.shape}")
+    if len(data_tr) != len(y_tr):
+        raise ValueError(f"data_train and y_train must have same length, got {len(data_tr)} and {len(y_tr)}")
+    if data_tr.shape[1] != data_te.shape[1]:
+        raise ValueError(f"data_train and data_test must have same number of dimensions, got {data_tr.shape[1]} and {data_te.shape[1]}")
+    if kernel not in _KERNELS:
+        raise ValueError(f"kernel must be one of {list(_KERNELS.keys())}, got {kernel!r}")
+
+    K, _, _, _, _, _ = _KERNELS[kernel]
+
+    U = (data_te[:, None, :] - data_tr[None, :, :]) / h
+    K_vals = K(U)
+    w = np.prod(K_vals, axis=2)
+
+    w_sum = w.sum(axis=1)
+    w_sum_safe = np.where(w_sum == 0, np.finfo(float).eps, w_sum)
+    y_pred = (w @ y_tr) / w_sum_safe
+
+    zero_weight_mask = w_sum == 0
+    if np.any(zero_weight_mask):
+        y_pred[zero_weight_mask] = np.mean(y_tr)
+
+    return y_pred
+
+
 def nw_bandwidth_mv(
     data: ArrayLike,
     y: ArrayLike,

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "hatchling.build"
 
 [project]
 name = "hbw"
-version = "0.5.0"
+version = "0.6.0"
 description = "Fast kernel bandwidth selection via analytic Hessian Newton optimization"
 readme = "README.md"
 license = "MIT"