API Reference

beta_kde.BetaKDE

Bases: DensityMixin, BaseEstimator

Beta Kernel Density Estimation with Scikit-learn API compatibility.

This estimator is designed for data strictly bounded within a fixed support (default [0, 1]). It addresses the Boundary Bias problem common in Gaussian KDEs by using Beta distributions as kernels.

Parameters:

Name	Type	Description	Default
bandwidth	float, str, or None	The bandwidth selection method for the MARGINALS. Options: float, 'beta-reference', 'LCV', 'LSCV'.	None
bounds	tuple of float	The strict support of the data (min, max).	(0.0, 1.0)
bandwidth_bounds	tuple of float	The search range (min_h, max_h) used when bandwidth is set to 'LCV' or 'LSCV'.	(0.01, 0.2)
selection_grid_points	int	Points for LSCV grid search.	30
heuristic_factor	float	Expansion factor for LSCV heuristic search.	4.0
integration_points	int	Points used for numerical integration in LSCV.	200
copula_grid_size	int	Resolution of the grid used for Copula transformation.	1000
verbose	int	Verbosity level.	0

Source code in src/beta_kde/estimator.py

class BetaKDE(DensityMixin, BaseEstimator):
    r"""
    Beta Kernel Density Estimation with Scikit-learn API compatibility.

    This estimator is designed for data strictly bounded within a fixed support
    (default [0, 1]). It addresses the **Boundary Bias** problem common in
    Gaussian KDEs by using Beta distributions as kernels.

    Parameters
    ----------
    bandwidth : float, str, or None, default=None
        The bandwidth selection method for the MARGINALS.
        Options: float, 'beta-reference', 'LCV', 'LSCV'.

    bounds : tuple of float, default=(0.0, 1.0)
        The strict support of the data (min, max).

    bandwidth_bounds : tuple of float, default=(0.01, 0.2)
        The search range (min_h, max_h) used when `bandwidth` is set to 'LCV' or 'LSCV'.

    selection_grid_points : int, default=30
        Points for LSCV grid search.

    heuristic_factor : float, default=4.0
        Expansion factor for LSCV heuristic search.

    integration_points : int, default=200
        Points used for numerical integration in LSCV.

    copula_grid_size : int, default=1000
        Resolution of the grid used for Copula transformation.

    verbose : int, default=0
        Verbosity level.
    """

    VALID_SELECTION_METHODS = ["LCV", "LSCV", "beta-reference"]

    def __init__(
        self,
        bandwidth: float | str | None = None,
        bounds: tuple[float, float] = (0.0, 1.0),
        bandwidth_bounds: tuple[float, float] = (0.01, 0.2),
        selection_grid_points: int = 30,
        heuristic_factor: float = 4.0,
        integration_points: int = 200,
        copula_grid_size: int = 1000,
        verbose: int = 0,
    ):
        self.bandwidth = bandwidth
        self.bounds = bounds
        self.bandwidth_bounds = bandwidth_bounds
        self.selection_grid_points = selection_grid_points
        self.heuristic_factor = heuristic_factor
        self.integration_points = integration_points
        self.copula_grid_size = copula_grid_size
        self.verbose = verbose

    def __sklearn_tags__(self):
        tags = super().__sklearn_tags__()
        tags.input_tags.positive_only = self.bounds[0] >= 0
        tags.input_tags.one_d_array = False
        tags.input_tags.two_d_array = True
        tags.target_tags.required = False
        return tags

    # --- Class attributes (type hints for mypy) ---
    # Attributes set after fitting
    n_samples_: int
    n_features_: int
    n_features_in_: int
    training_data_: np.ndarray
    scale_factor_: float
    shift_: float
    _epsilon: float
    data_clipped_: np.ndarray
    marginal_bandwidths_: list[float]
    x_grids_: list[np.ndarray]
    cdf_grids_: list[np.ndarray]
    bandwidth_: float | None
    is_fallback_: bool | None
    copula_bandwidth_: float | None
    normalization_constant_: float | None
    is_fitted_: bool

    def fit(
        self,
        X: ArrayLike,
        y: Any = None,
        compute_normalization: bool = False,
    ) -> BetaKDE:
        """
        Fit the Beta Kernel Density model to the training data.

        Parameters
        ----------
        X : array-like of shape (n_samples, n_features)
            Training data.
        y : Ignored
        compute_normalization : bool, default=False
            If True, triggers the lazy calculation of the normalization constant
            immediately. Useful for performance benchmarking or if you know
            you will need normalized scores later.

        Returns
        -------
        self : BetaKDE
            Fitted estimator.
        """
        # Reset attributes
        self.bandwidth_ = None
        self.is_fallback_ = None
        # Reset normalization constant to None (invalidating previous fit)
        self.normalization_constant_ = None

        X = check_array(X, ensure_2d=True, order="C", dtype=np.float64)

        self.n_samples_, self.n_features_ = X.shape
        self.training_data_ = X
        self.n_features_in_ = self.n_features_

        lower, upper = self.bounds
        if lower >= upper:
            raise ValueError(f"Bounds must be strictly increasing. Got {self.bounds}")

        if not np.all((X >= lower) & (X <= upper)):
            raise ValueError(
                f"All data points must be within the interval {self.bounds}. "
                f"Found range [{X.min():.3f}, {X.max():.3f}]."
            )

        if isinstance(self.bandwidth, (float, int)) and not isinstance(
            self.bandwidth, bool
        ):
            if self.bandwidth <= 0:
                raise ValueError("Bandwidth must be positive.")
        elif isinstance(self.bandwidth, str):
            if self.bandwidth not in self.VALID_SELECTION_METHODS:
                raise ValueError(
                    f"Unknown bandwidth selection method: '{self.bandwidth}'"
                )

        # Scale Data to [0, 1]
        self.scale_factor_ = upper - lower
        self.shift_ = lower

        X_scaled = (X - self.shift_) / self.scale_factor_
        self._epsilon = 1e-10
        self.data_clipped_ = np.clip(X_scaled, self._epsilon, 1.0 - self._epsilon)

        # Fit Marginals
        self.marginal_bandwidths_ = []
        fallback_statuses = []
        self.cdf_grids_ = []
        self.x_grids_ = []

        for d in range(self.n_features_):
            data_d = self.data_clipped_[:, d]

            if self.verbose > 0 and self.n_features_ > 1:
                print(f"Fitting Dimension {d+1}/{self.n_features_}...")

            h, is_fb = self._select_bandwidth_for_dim(data_d)
            self.marginal_bandwidths_.append(h)
            fallback_statuses.append(is_fb)

            if self.verbose > 0:
                if is_fb:
                    print(
                        f"  Dim {d+1}: MISE rule failed constraints. Using fallback: h = {h:.4f}"
                    )
                elif self.n_features_ > 1:
                    print(f"  Dim {d+1}: Bandwidth selected: h = {h:.4f}")

            # Pre-compute CDF for Copula transform
            if self.n_features_ > 1:
                grid = np.linspace(0, 1, self.copula_grid_size)
                log_pdf = self._score_samples_1d(grid, data_d, h)
                pdf = np.exp(log_pdf)
                cdf: np.ndarray = np.cumsum(pdf)  # Type annotation for mypy
                cdf = cdf / cdf[-1]  # Normalize
                self.x_grids_.append(grid)
                self.cdf_grids_.append(cdf)

        if self.n_features_ == 1:
            self.bandwidth_ = self.marginal_bandwidths_[0]
            self.is_fallback_ = fallback_statuses[0]
            if self.verbose > 0:
                if self.is_fallback_:
                    print(
                        f"MISE rule failed constraints. Using fallback: h = {self.bandwidth_:.4f}"
                    )
                else:
                    print(f"Bandwidth selected by MISE rule: h = {self.bandwidth_:.4f}")

        # Copula Bandwidth
        if self.n_features_ > 1:
            self.U_train_ = self._transform_to_uniform(self.data_clipped_)
            self.copula_bandwidth_ = self.n_samples_ ** (-1.0 / (self.n_features_ + 4))

            if self.verbose > 0:
                print(f"Copula Bandwidth (Scott's Rule): {self.copula_bandwidth_:.4f}")

        self.is_fitted_ = True

        # Trigger lazy computation if explicitly requested
        if compute_normalization:
            _ = self.normalization_constant

        return self

    @property
    def normalization_constant(self) -> float:
        """
        The normalization constant of the density.
        Computed lazily via numerical integration upon first access.
        """
        check_is_fitted(self)
        if self.normalization_constant_ is None:
            self.normalization_constant_ = self._compute_normalization_constant()
        return self.normalization_constant_

    def _normalization_integrand(self, x_val, h, data_d):
        """Integrand helper method."""
        if np.ndim(x_val) == 0:
            x_val = np.array([x_val])

        mask = (x_val > 0) & (x_val < 1)
        if not np.any(mask):
            return 0.0 if np.ndim(x_val) == 0 else np.zeros_like(x_val)

        x_valid = x_val[mask]
        k_mat = self._kernel_matrix(x_valid, data_d, h)

        res = np.zeros_like(x_val)
        res[mask] = np.mean(k_mat, axis=1)

        return res if res.size > 1 else res.item()

    def _compute_normalization_constant(self) -> float:
        """Internal worker to compute and cache the normalization constant."""
        if self.n_samples_ < 2:
            raise ValueError(
                f"Cannot compute normalization constant with only {self.n_samples_} sample(s). "
                "At least 2 samples are required for numerical integration. "
                "Use score_samples(normalized=False) for unnormalized densities."
            )

        marginal_constants = []
        for d in range(self.n_features_):
            h = self.marginal_bandwidths_[d]
            data_d = self.data_clipped_[:, d]

            # Check for constant data in this dimension
            if np.std(data_d) < 1e-10:
                raise ValueError(
                    f"Dimension {d+1} has constant data (std < 1e-10). "
                    f"Cannot compute normalization constant for constant data. "
                    f"Use score_samples(normalized=False) instead."
                )

            integral, _ = scipy.integrate.quad(
                self._normalization_integrand,
                0,
                1,
                args=(h, data_d),
                epsabs=1e-4,
                limit=50,
            )

            if not np.isfinite(integral) or integral <= 0:
                raise RuntimeError(
                    f"Dimension {d+1}: Integration failed with result={integral}. "
                    f"Bandwidth: h={h:.6f}. "
                    "This may be due to numerical issues or extreme parameter values. "
                    "Try using a different bandwidth method or fewer dimensions."
                )

            marginal_constants.append(integral)

        result = np.prod(marginal_constants)
        if not np.isfinite(result) or result <= 0:
            raise RuntimeError(
                f"Normalization constant computation failed: result={result}. "
                "This may be due to numerical issues with the integration. "
                "Try using a different bandwidth method or fewer dimensions."
            )

        return result

    def score_samples(self, X: ArrayLike, normalized: bool = False) -> np.ndarray:
        """
        Compute the log-likelihood of each sample.

        Parameters
        ----------
        X : array-like of shape (n_samples, n_features)
            Data to score.
        normalized : bool, default=False
            If True, ensures the density integrates to 1.0. 
            This triggers numerical integration if not yet computed.

        Returns
        -------
        log_density : ndarray of shape (n_samples,)
            Log-likelihood of each sample.
        """
        check_is_fitted(self)

        if np.ndim(X) == 0:
            X = np.array([[X]])
        X = check_array(X, ensure_2d=True, order="C", dtype=np.float64)

        if hasattr(self, "n_features_in_"):
            if X.shape[1] != self.n_features_in_:
                raise ValueError(
                    f"X has {X.shape[1]} features, but BetaKDE is expecting "
                    f"{self.n_features_in_} features as input."
                )
        elif X.shape[1] != self.n_features_:
            raise ValueError(
                f"Mismatch in dimensions. Model: {self.n_features_}, Data: {X.shape[1]}"
            )

        X_scaled = (X - self.shift_) / self.scale_factor_
        X_safe = np.clip(X_scaled, self._epsilon, 1.0 - self._epsilon)

        n_test = X.shape[0]
        log_density = np.zeros(n_test)

        # 1. Marginal Log-Likelihoods
        for d in range(self.n_features_):
            h = self.marginal_bandwidths_[d]
            train_d = self.data_clipped_[:, d]
            log_pdf_scaled = self._score_samples_1d(X_safe[:, d], train_d, h)
            log_pdf = log_pdf_scaled - np.log(self.scale_factor_)
            log_density += log_pdf

        # 2. Copula Log-Likelihood
        if self.n_features_ > 1:
            U_test = self._transform_to_uniform(X_safe)
            log_copula = self._score_copula(
                U_test, self.U_train_, self.copula_bandwidth_
            )
            log_density += log_copula

        # 3. Normalization (Lazy)
        if normalized:
            log_norm = np.log(self.normalization_constant)
            log_density -= log_norm

        return log_density

    def score(self, X: ArrayLike, y: Any = None) -> float:
        """
        Compute the total log-likelihood under the model.

        **Note:** This method explicitly forces `normalized=True` to ensure
        statistical validity when used in cross-validation (e.g., GridSearchCV).
        To get raw scores, use `score_samples(X, normalized=False).sum()`.

        Parameters
        ----------
        X : array-like of shape (n_samples, n_features)
            Data to score.
        y : Ignored

        Returns
        -------
        score : float
            Total log-likelihood.
        """
        return np.sum(self.score_samples(X, normalized=True))

    def pdf(self, X: ArrayLike, normalized: bool = False) -> float | np.ndarray:
        """
        Convenience method returning the probability density (exp(score_samples)).

        Parameters
        ----------
        X : array-like of shape (n_samples, n_features)
            Data to evaluate.
        normalized : bool, default=False
            If True, ensures the density integrates to 1.0.

        Returns
        -------
        pdf : float or ndarray of shape (n_samples,)
            Probability density values.
        """
        is_scalar = np.ndim(X) == 0
        if is_scalar:
            X_arg = np.array([[X]])
        elif np.ndim(X) == 1:
            X_arg = X.reshape(-1, 1)
        else:
            X_arg = X

        log_pdf = self.score_samples(X_arg, normalized=normalized)
        pdf_vals = np.exp(log_pdf)

        if is_scalar:
            return float(pdf_vals[0])
        return pdf_vals

    def plot(
        self,
        eval_points: np.ndarray | None = None,
        show_histogram: bool = True,
        bins: int = 20,
        normalized: bool = False,
        ax: Any | None = None,
        label: str | None = None,
        **kwargs: Any,
    ) -> Any | tuple[Any, Any]:
        """
        Plots the estimated Marginal Probability Density Functions (PDFs).
        """
        import matplotlib.pyplot as plt

        check_is_fitted(self)
        lower, upper = self.bounds
        n_dims = self.n_features_

        if ax is None:
            if n_dims == 1:
                fig, ax = plt.subplots(figsize=(10, 6))
                axes_list = [ax]
            else:
                cols = int(np.ceil(np.sqrt(n_dims)))
                rows = int(np.ceil(n_dims / cols))
                fig, axs = plt.subplots(rows, cols, figsize=(5 * cols, 4 * rows))
                axes_list = axs.flatten()
        else:
            if n_dims == 1:
                try:
                    fig = ax.figure
                except AttributeError:
                    fig = None
                axes_list = [ax]
            else:
                if not isinstance(ax, (list, np.ndarray)):
                    warnings.warn(
                        "Multivariate plot requested but single axis provided. Plotting 1st dimension only."
                    )
                    axes_list = [ax]
                else:
                    axes_list = np.array(ax).flatten()
                    try:
                        fig = axes_list[0].figure
                    except AttributeError:
                        fig = None

        for d in range(len(axes_list)):
            if d >= n_dims:
                axes_list[d].axis("off")
                continue

            curr_ax = axes_list[d]
            data_d = self.training_data_[:, d]
            h = self.marginal_bandwidths_[d]

            if eval_points is None:
                x_plot = np.linspace(lower, upper, 1000)
            else:
                x_plot = eval_points

            x_scaled = (x_plot - self.shift_) / self.scale_factor_
            x_safe = np.clip(x_scaled, self._epsilon, 1.0 - self._epsilon)
            train_d = self.data_clipped_[:, d]

            log_pdf_scaled = self._score_samples_1d(x_safe, train_d, h)
            log_pdf = log_pdf_scaled - np.log(self.scale_factor_)
            pdf_vals = np.exp(log_pdf)

            # NOTE: For plotting, we use trapezoidal rule normalization.
            # This is robust for multivariate marginals where the total
            # normalization constant != marginal constant.
            if normalized:
                integral = np.trapezoid(pdf_vals, x_plot)
                if integral > 0:
                    pdf_vals /= integral

            if n_dims == 1:
                plot_label = f"Beta KDE (h={h:.3f})" if label is None else label
                curr_ax.set_title("Beta Kernel Density Estimation")
            else:
                plot_label = f"Dim {d+1} (h={h:.3f})" if label is None else label
                curr_ax.set_title(f"Dimension {d+1}")

            if normalized:
                plot_label += " [Norm]"

            curr_ax.plot(x_plot, pdf_vals, label=plot_label, **kwargs)

            if show_histogram:
                curr_ax.hist(
                    data_d,
                    bins=bins,
                    density=True,
                    alpha=0.5,
                    color="gray",
                    edgecolor="none",
                    range=(lower, upper),
                )

            curr_ax.set_xlim(lower, upper)
            curr_ax.set_ylim(bottom=0)
            curr_ax.legend()

        if fig is not None:
            plt.tight_layout()
            return fig, axes_list if n_dims > 1 else axes_list[0]
        return axes_list if n_dims > 1 else axes_list[0]

    def _select_bandwidth_for_dim(self, data_1d):
        method = self.bandwidth if self.bandwidth else "beta-reference"
        if isinstance(method, (float, int)):
            return float(method), False
        if method == "LCV":
            return self._select_bandwidth_lcv(data_1d, self.bandwidth_bounds), False
        elif method == "LSCV":
            return (
                self._select_bandwidth_lscv(
                    data_1d,
                    self.bandwidth_bounds,
                    self.selection_grid_points,
                    self.heuristic_factor,
                    self.integration_points,
                ),
                False,
            )
        elif method == "beta-reference":
            return self._select_bandwidth_beta_reference(data_1d)
        raise ValueError(f"Unknown method: {method}")

    def _score_samples_1d(self, x_eval, data_train, h):
        k_mat = self._kernel_matrix(x_eval, data_train, h)
        pdf_vals = np.mean(k_mat, axis=1)
        return np.log(pdf_vals + 1e-300)

    def _transform_to_uniform(self, X_scaled):
        U = np.zeros_like(X_scaled)
        for d in range(self.n_features_):
            U[:, d] = np.interp(X_scaled[:, d], self.x_grids_[d], self.cdf_grids_[d])
        return np.clip(U, 1e-5, 1 - 1e-5)

    def _score_copula(self, U_test, U_train, h):
        n_test = U_test.shape[0]
        n_train = U_train.shape[0]
        d_dims = U_test.shape[1]
        log_weights = np.zeros((n_test, n_train))
        for j in range(d_dims):
            k_mat_j = self._kernel_matrix(U_test[:, j], U_train[:, j], h)
            log_weights += np.log(k_mat_j + 1e-300)
        max_log = np.max(log_weights, axis=1)
        sum_exp = np.sum(np.exp(log_weights - max_log[:, None]), axis=1)
        log_copula = max_log + np.log(sum_exp + 1e-300) - np.log(n_train)
        return log_copula

    def _lcv_objective(self, bandwidth, data):
        if not (0 < bandwidth < 1):
            return np.inf
        n = len(data)
        K_mat = self._kernel_matrix(data, data, bandwidth)
        row_sums = K_mat.sum(axis=1)
        diag_elems = np.diag(K_mat)
        f_hat_loo = (row_sums - diag_elems) / (n - 1)
        f_hat_loo = np.maximum(f_hat_loo, 1e-10)
        return -np.sum(np.log(f_hat_loo))

    def _select_bandwidth_lcv(self, data, bounds):
        res = scipy.optimize.minimize_scalar(
            lambda h: self._lcv_objective(h, data), bounds=bounds, method="bounded"
        )
        if res.success:
            return float(res.x)
        raise RuntimeError("LCV failed")

    def _lscv_objective(self, bandwidth, data, integration_points):
        if not (0 < bandwidth < 1):
            return np.inf
        n = len(data)
        x_grid = np.linspace(1e-5, 1.0 - 1e-5, integration_points)
        K_grid = self._kernel_matrix(x_grid, data, bandwidth)
        pdf_grid = K_grid.mean(axis=1)
        term1 = scipy.integrate.trapezoid(pdf_grid**2, x_grid)
        K_data = self._kernel_matrix(data, data, bandwidth)
        term2 = (np.sum(K_data) - np.sum(np.diag(K_data))) * (-2 / (n * (n - 1)))
        return term1 + term2

    def _select_bandwidth_lscv(
        self, data, bounds, grid_points, heuristic_factor, integration_points
    ):
        std_dev = np.std(data, ddof=0)
        n = len(data)
        search_bounds = bounds
        if std_dev > 1e-8:
            h_rule = 0.9 * std_dev * (n ** (-0.2))
            search_bounds = (
                max(bounds[0], h_rule / heuristic_factor),
                min(bounds[1], h_rule * heuristic_factor),
            )
        h_grid = np.linspace(search_bounds[0], search_bounds[1], grid_points)
        scores = [self._lscv_objective(h, data, integration_points) for h in h_grid]
        best_h = h_grid[np.nanargmin(scores)]

        step = h_grid[1] - h_grid[0] if grid_points > 1 else 0.01
        refine_bounds = (max(bounds[0], best_h - step), min(bounds[1], best_h + step))

        res = scipy.optimize.minimize_scalar(
            lambda h: self._lscv_objective(h, data, integration_points),
            bounds=refine_bounds,
            method="bounded",
        )
        return float(res.x) if res.success else best_h

    def _select_bandwidth_beta_reference(self, data):
        X_filtered = data[(data > 0) & (data < 1)]
        h_final = 0.1
        is_fallback = False

        try:
            ahat, bhat = self._estimate_beta_params(X_filtered)
            if not (ahat > 1.5 and bhat > 1.5 and (ahat + bhat) > 3):
                raise ValueError("Parameters too small for MISE rule.")

            a, b, n = ahat, bhat, len(data)
            log_num = (
                np.log(2 * a + 2 * b - 5)
                + np.log(2 * a + 2 * b - 3)
                + sp.gammaln(2 * a + 2 * b - 6)
                + sp.gammaln(a)
                + sp.gammaln(b)
                + sp.gammaln(a - 0.5)
                + sp.gammaln(b - 0.5)
            )
            denom_term_1 = (a - 1) * (b - 1)
            denom_term_2 = 6 - 4 * b + a * (3 * b - 4)
            if denom_term_1 <= 0 or denom_term_2 <= 0:
                raise ValueError("Denominator factor non-positive.")
            log_denom = (
                np.log(denom_term_1)
                + np.log(denom_term_2)
                + sp.gammaln(2 * a - 3)
                + sp.gammaln(2 * b - 3)
                + sp.gammaln(a + b)
                + sp.gammaln(a + b - 1)
            )
            log_factor = np.log(2) + np.log(n) + 0.5 * np.log(np.pi)
            log_h = (2 / 5) * (log_num - log_denom - log_factor)
            h_final = np.exp(log_h)
            if not (0 < h_final < 1):
                raise ValueError("Calculated bandwidth outside (0, 1).")

        except (ValueError, RuntimeError) as e:
            if ("Sample variance is zero" in str(e) or "too large" in str(e)) and len(
                data
            ) > 1:
                raise

            if not (hasattr(self, "ahat_") and hasattr(self, "bhat_")):
                try:
                    self._estimate_beta_params(X_filtered)
                except ValueError:
                    return 1.0 * (len(data) ** (-0.4)), True
            h_final = self._calculate_hybrid_fallback(self.ahat_, self.bhat_, len(data))
            is_fallback = True
            if self.verbose > 0:
                warnings.warn(f"MISE Rule failed: {e}. Using fallback.", RuntimeWarning)
        return h_final, is_fallback

    def _estimate_beta_params(self, X_filtered: np.ndarray) -> tuple[float, float]:
        if X_filtered.size == 0:
            raise ValueError(
                "No data strictly within (0, 1). "
                "Beta kernel bandwidth selection requires at least some data "
                "points strictly inside the open interval (0, 1). "
                "Try clipping your data or using a different bandwidth method."
            )
        mean_x = np.mean(X_filtered)
        var_x = np.var(X_filtered, ddof=0)

        if var_x == 0:
            raise ValueError(
                f"Sample variance is zero. "
                f"All {len(X_filtered)} data points have the same value: {X_filtered[0]:.6f}. "
                "Cannot estimate Beta parameters from constant data. "
                "Try using a fixed bandwidth value instead."
            )

        max_var = mean_x * (1 - mean_x)
        if var_x >= max_var:
            raise ValueError(
                f"Sample variance ({var_x:.6f}) exceeds the maximum possible "
                f"for Beta distribution with mean {mean_x:.6f}. "
                f"Maximum allowed variance: {max_var:.6f}. "
                "This typically occurs when data is concentrated at boundaries. "
                "Consider using the 'fallback' bandwidth method."
            )

        common = ((mean_x * (1 - mean_x)) / var_x) - 1
        a, b = mean_x * common, (1 - mean_x) * common
        if a <= 0 or b <= 0:
            raise ValueError(
                f"Estimated Beta parameters are not positive: a={a:.6f}, b={b:.6f}. "
                f"Mean: {mean_x:.6f}, Variance: {var_x:.6f}. "
                "This indicates numerical instability in parameter estimation. "
                "Try using a different bandwidth method or increasing your dataset."
            )
        self.ahat_, self.bhat_ = a, b
        return a, b

    def _calculate_hybrid_fallback(self, a, b, n):
        s = np.sqrt(self._variance(a, b))
        correction = 1 + abs(self._skewness(a, b)) + abs(self._kurtosis(a, b))
        return (s / correction) * (n ** (-0.4)) if s > 0 else 1e-5

    @staticmethod
    def _skewness(a, b):
        return (2 * (b - a) * np.sqrt(a + b + 1)) / ((a + b + 2) * np.sqrt(a * b))

    @staticmethod
    def _kurtosis(a, b):
        num = 6 * ((a - b) ** 2 * (a + b + 1) - a * b * (a + b + 2))
        den = a * b * (a + b + 2) * (a + b + 3)
        return num / den

    @staticmethod
    def _variance(a, b):
        return (a * b) / ((a + b) ** 2 * (a + b + 1))

    def _rho_vec(self, x_arr, bandwidth):
        h = bandwidth
        term2 = np.maximum(4 * h**4 + 6 * h**2 + 2.25 - x_arr**2 - x_arr / h, 0.0)
        return (2 * h**2 + 2.5) - np.sqrt(term2)

    def _kernel_matrix(self, x_eval, data_pts, bandwidth):
        n_eval = x_eval.shape[0]
        x_col = x_eval.reshape(n_eval, 1)
        h = bandwidth
        lower_thresh, upper_thresh = 2 * h, 1 - 2 * h
        alpha = x_col / h
        beta_p = (1 - x_col) / h
        alpha = np.where(x_col < lower_thresh, self._rho_vec(x_col, h), alpha)
        beta_p = np.where(x_col > upper_thresh, self._rho_vec(1 - x_col, h), beta_p)
        return beta_dist.pdf(data_pts[np.newaxis, :], alpha, beta_p)

normalization_constant property

The normalization constant of the density. Computed lazily via numerical integration upon first access.

fit(X, y=None, compute_normalization=False)

Fit the Beta Kernel Density model to the training data.

Parameters:

Name	Type	Description	Default
X	array-like of shape (n_samples, n_features)	Training data.	required
y	Ignored		None
compute_normalization	bool	If True, triggers the lazy calculation of the normalization constant immediately. Useful for performance benchmarking or if you know you will need normalized scores later.	False

Returns:

Name	Type	Description
self	BetaKDE	Fitted estimator.

Source code in src/beta_kde/estimator.py

def fit(
    self,
    X: ArrayLike,
    y: Any = None,
    compute_normalization: bool = False,
) -> BetaKDE:
    """
    Fit the Beta Kernel Density model to the training data.

    Parameters
    ----------
    X : array-like of shape (n_samples, n_features)
        Training data.
    y : Ignored
    compute_normalization : bool, default=False
        If True, triggers the lazy calculation of the normalization constant
        immediately. Useful for performance benchmarking or if you know
        you will need normalized scores later.

    Returns
    -------
    self : BetaKDE
        Fitted estimator.
    """
    # Reset attributes
    self.bandwidth_ = None
    self.is_fallback_ = None
    # Reset normalization constant to None (invalidating previous fit)
    self.normalization_constant_ = None

    X = check_array(X, ensure_2d=True, order="C", dtype=np.float64)

    self.n_samples_, self.n_features_ = X.shape
    self.training_data_ = X
    self.n_features_in_ = self.n_features_

    lower, upper = self.bounds
    if lower >= upper:
        raise ValueError(f"Bounds must be strictly increasing. Got {self.bounds}")

    if not np.all((X >= lower) & (X <= upper)):
        raise ValueError(
            f"All data points must be within the interval {self.bounds}. "
            f"Found range [{X.min():.3f}, {X.max():.3f}]."
        )

    if isinstance(self.bandwidth, (float, int)) and not isinstance(
        self.bandwidth, bool
    ):
        if self.bandwidth <= 0:
            raise ValueError("Bandwidth must be positive.")
    elif isinstance(self.bandwidth, str):
        if self.bandwidth not in self.VALID_SELECTION_METHODS:
            raise ValueError(
                f"Unknown bandwidth selection method: '{self.bandwidth}'"
            )

    # Scale Data to [0, 1]
    self.scale_factor_ = upper - lower
    self.shift_ = lower

    X_scaled = (X - self.shift_) / self.scale_factor_
    self._epsilon = 1e-10
    self.data_clipped_ = np.clip(X_scaled, self._epsilon, 1.0 - self._epsilon)

    # Fit Marginals
    self.marginal_bandwidths_ = []
    fallback_statuses = []
    self.cdf_grids_ = []
    self.x_grids_ = []

    for d in range(self.n_features_):
        data_d = self.data_clipped_[:, d]

        if self.verbose > 0 and self.n_features_ > 1:
            print(f"Fitting Dimension {d+1}/{self.n_features_}...")

        h, is_fb = self._select_bandwidth_for_dim(data_d)
        self.marginal_bandwidths_.append(h)
        fallback_statuses.append(is_fb)

        if self.verbose > 0:
            if is_fb:
                print(
                    f"  Dim {d+1}: MISE rule failed constraints. Using fallback: h = {h:.4f}"
                )
            elif self.n_features_ > 1:
                print(f"  Dim {d+1}: Bandwidth selected: h = {h:.4f}")

        # Pre-compute CDF for Copula transform
        if self.n_features_ > 1:
            grid = np.linspace(0, 1, self.copula_grid_size)
            log_pdf = self._score_samples_1d(grid, data_d, h)
            pdf = np.exp(log_pdf)
            cdf: np.ndarray = np.cumsum(pdf)  # Type annotation for mypy
            cdf = cdf / cdf[-1]  # Normalize
            self.x_grids_.append(grid)
            self.cdf_grids_.append(cdf)

    if self.n_features_ == 1:
        self.bandwidth_ = self.marginal_bandwidths_[0]
        self.is_fallback_ = fallback_statuses[0]
        if self.verbose > 0:
            if self.is_fallback_:
                print(
                    f"MISE rule failed constraints. Using fallback: h = {self.bandwidth_:.4f}"
                )
            else:
                print(f"Bandwidth selected by MISE rule: h = {self.bandwidth_:.4f}")

    # Copula Bandwidth
    if self.n_features_ > 1:
        self.U_train_ = self._transform_to_uniform(self.data_clipped_)
        self.copula_bandwidth_ = self.n_samples_ ** (-1.0 / (self.n_features_ + 4))

        if self.verbose > 0:
            print(f"Copula Bandwidth (Scott's Rule): {self.copula_bandwidth_:.4f}")

    self.is_fitted_ = True

    # Trigger lazy computation if explicitly requested
    if compute_normalization:
        _ = self.normalization_constant

    return self

pdf(X, normalized=False)

Convenience method returning the probability density (exp(score_samples)).

Parameters:

Name	Type	Description	Default
X	array-like of shape (n_samples, n_features)	Data to evaluate.	required
normalized	bool	If True, ensures the density integrates to 1.0.	False

Returns:

Name	Type	Description
pdf	float or ndarray of shape (n_samples,)	Probability density values.

Source code in src/beta_kde/estimator.py

def pdf(self, X: ArrayLike, normalized: bool = False) -> float | np.ndarray:
    """
    Convenience method returning the probability density (exp(score_samples)).

    Parameters
    ----------
    X : array-like of shape (n_samples, n_features)
        Data to evaluate.
    normalized : bool, default=False
        If True, ensures the density integrates to 1.0.

    Returns
    -------
    pdf : float or ndarray of shape (n_samples,)
        Probability density values.
    """
    is_scalar = np.ndim(X) == 0
    if is_scalar:
        X_arg = np.array([[X]])
    elif np.ndim(X) == 1:
        X_arg = X.reshape(-1, 1)
    else:
        X_arg = X

    log_pdf = self.score_samples(X_arg, normalized=normalized)
    pdf_vals = np.exp(log_pdf)

    if is_scalar:
        return float(pdf_vals[0])
    return pdf_vals

plot(eval_points=None, show_histogram=True, bins=20, normalized=False, ax=None, label=None, **kwargs)

Plots the estimated Marginal Probability Density Functions (PDFs).

Source code in src/beta_kde/estimator.py

def plot(
    self,
    eval_points: np.ndarray | None = None,
    show_histogram: bool = True,
    bins: int = 20,
    normalized: bool = False,
    ax: Any | None = None,
    label: str | None = None,
    **kwargs: Any,
) -> Any | tuple[Any, Any]:
    """
    Plots the estimated Marginal Probability Density Functions (PDFs).
    """
    import matplotlib.pyplot as plt

    check_is_fitted(self)
    lower, upper = self.bounds
    n_dims = self.n_features_

    if ax is None:
        if n_dims == 1:
            fig, ax = plt.subplots(figsize=(10, 6))
            axes_list = [ax]
        else:
            cols = int(np.ceil(np.sqrt(n_dims)))
            rows = int(np.ceil(n_dims / cols))
            fig, axs = plt.subplots(rows, cols, figsize=(5 * cols, 4 * rows))
            axes_list = axs.flatten()
    else:
        if n_dims == 1:
            try:
                fig = ax.figure
            except AttributeError:
                fig = None
            axes_list = [ax]
        else:
            if not isinstance(ax, (list, np.ndarray)):
                warnings.warn(
                    "Multivariate plot requested but single axis provided. Plotting 1st dimension only."
                )
                axes_list = [ax]
            else:
                axes_list = np.array(ax).flatten()
                try:
                    fig = axes_list[0].figure
                except AttributeError:
                    fig = None

    for d in range(len(axes_list)):
        if d >= n_dims:
            axes_list[d].axis("off")
            continue

        curr_ax = axes_list[d]
        data_d = self.training_data_[:, d]
        h = self.marginal_bandwidths_[d]

        if eval_points is None:
            x_plot = np.linspace(lower, upper, 1000)
        else:
            x_plot = eval_points

        x_scaled = (x_plot - self.shift_) / self.scale_factor_
        x_safe = np.clip(x_scaled, self._epsilon, 1.0 - self._epsilon)
        train_d = self.data_clipped_[:, d]

        log_pdf_scaled = self._score_samples_1d(x_safe, train_d, h)
        log_pdf = log_pdf_scaled - np.log(self.scale_factor_)
        pdf_vals = np.exp(log_pdf)

        # NOTE: For plotting, we use trapezoidal rule normalization.
        # This is robust for multivariate marginals where the total
        # normalization constant != marginal constant.
        if normalized:
            integral = np.trapezoid(pdf_vals, x_plot)
            if integral > 0:
                pdf_vals /= integral

        if n_dims == 1:
            plot_label = f"Beta KDE (h={h:.3f})" if label is None else label
            curr_ax.set_title("Beta Kernel Density Estimation")
        else:
            plot_label = f"Dim {d+1} (h={h:.3f})" if label is None else label
            curr_ax.set_title(f"Dimension {d+1}")

        if normalized:
            plot_label += " [Norm]"

        curr_ax.plot(x_plot, pdf_vals, label=plot_label, **kwargs)

        if show_histogram:
            curr_ax.hist(
                data_d,
                bins=bins,
                density=True,
                alpha=0.5,
                color="gray",
                edgecolor="none",
                range=(lower, upper),
            )

        curr_ax.set_xlim(lower, upper)
        curr_ax.set_ylim(bottom=0)
        curr_ax.legend()

    if fig is not None:
        plt.tight_layout()
        return fig, axes_list if n_dims > 1 else axes_list[0]
    return axes_list if n_dims > 1 else axes_list[0]

score(X, y=None)

Compute the total log-likelihood under the model.

Note: This method explicitly forces normalized=True to ensure statistical validity when used in cross-validation (e.g., GridSearchCV). To get raw scores, use score_samples(X, normalized=False).sum().

Parameters:

Name	Type	Description	Default
X	array-like of shape (n_samples, n_features)	Data to score.	required
y	Ignored		None

Returns:

Name	Type	Description
score	float	Total log-likelihood.

Source code in src/beta_kde/estimator.py

def score(self, X: ArrayLike, y: Any = None) -> float:
    """
    Compute the total log-likelihood under the model.

    **Note:** This method explicitly forces `normalized=True` to ensure
    statistical validity when used in cross-validation (e.g., GridSearchCV).
    To get raw scores, use `score_samples(X, normalized=False).sum()`.

    Parameters
    ----------
    X : array-like of shape (n_samples, n_features)
        Data to score.
    y : Ignored

    Returns
    -------
    score : float
        Total log-likelihood.
    """
    return np.sum(self.score_samples(X, normalized=True))

score_samples(X, normalized=False)

Compute the log-likelihood of each sample.

Parameters:

Name	Type	Description	Default
X	array-like of shape (n_samples, n_features)	Data to score.	required
normalized	bool	If True, ensures the density integrates to 1.0. This triggers numerical integration if not yet computed.	False

Returns:

Name	Type	Description
log_density	ndarray of shape (n_samples,)	Log-likelihood of each sample.

Source code in src/beta_kde/estimator.py

def score_samples(self, X: ArrayLike, normalized: bool = False) -> np.ndarray:
    """
    Compute the log-likelihood of each sample.

    Parameters
    ----------
    X : array-like of shape (n_samples, n_features)
        Data to score.
    normalized : bool, default=False
        If True, ensures the density integrates to 1.0. 
        This triggers numerical integration if not yet computed.

    Returns
    -------
    log_density : ndarray of shape (n_samples,)
        Log-likelihood of each sample.
    """
    check_is_fitted(self)

    if np.ndim(X) == 0:
        X = np.array([[X]])
    X = check_array(X, ensure_2d=True, order="C", dtype=np.float64)

    if hasattr(self, "n_features_in_"):
        if X.shape[1] != self.n_features_in_:
            raise ValueError(
                f"X has {X.shape[1]} features, but BetaKDE is expecting "
                f"{self.n_features_in_} features as input."
            )
    elif X.shape[1] != self.n_features_:
        raise ValueError(
            f"Mismatch in dimensions. Model: {self.n_features_}, Data: {X.shape[1]}"
        )

    X_scaled = (X - self.shift_) / self.scale_factor_
    X_safe = np.clip(X_scaled, self._epsilon, 1.0 - self._epsilon)

    n_test = X.shape[0]
    log_density = np.zeros(n_test)

    # 1. Marginal Log-Likelihoods
    for d in range(self.n_features_):
        h = self.marginal_bandwidths_[d]
        train_d = self.data_clipped_[:, d]
        log_pdf_scaled = self._score_samples_1d(X_safe[:, d], train_d, h)
        log_pdf = log_pdf_scaled - np.log(self.scale_factor_)
        log_density += log_pdf

    # 2. Copula Log-Likelihood
    if self.n_features_ > 1:
        U_test = self._transform_to_uniform(X_safe)
        log_copula = self._score_copula(
            U_test, self.U_train_, self.copula_bandwidth_
        )
        log_density += log_copula

    # 3. Normalization (Lazy)
    if normalized:
        log_norm = np.log(self.normalization_constant)
        log_density -= log_norm

    return log_density