`statspai.nonparametric`¶

nonparametric ¶

Nonparametric estimation methods.

Provides local polynomial regression (lpoly), kernel density estimation (kdensity), nonparametric regression (npregress), and conditional inference trees (ctree).

CTreeResult `dataclass` ¶

Bases: ResultProtocolMixin

Result of :func:ctree.

Attributes:

Name	Type	Description
`nodes`	`DataFrame`	One row per node, numbered depth first with the left child before the right one (the numbering of `partykit`): `parent`, `depth`, `n`, `variable` and `split` (inner nodes), `statistic` and `p_value` of the test that selected the variable, `prediction`, and `error` (share misclassified, or the within-node variance for a numeric outcome).
`tests`	`DataFrame`	Every test that was run: `node`, `variable`, `statistic`, `df`, `p_value_raw` and `p_value` (adjusted for the number of covariates tested in the node).
`fitted`	`ndarray`	Prediction for each row of the data.
`node_id`	`ndarray`	Terminal node of each row.
`task`	`{'classification', 'regression'}`
`classes`	`list or None`	Outcome classes, in the column order of `predict(type='prob')`.
`n_leaves, n_obs`	`int`

Examples:

>>> import numpy as np, pandas as pd, statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"x": rng.uniform(size=400), "z": rng.normal(size=400)})
>>> df["y"] = np.where(df["x"] > 0.5, 2.0, 0.0) + rng.normal(size=400)
>>> fit = sp.ctree(df, y="y", covariates=["x", "z"])
>>> isinstance(fit, sp.CTreeResult)
True
>>> fit.nodes.loc[1, "variable"]
'x'
>>> fit.predict(df.head(3)).shape
(3,)

predict ¶

predict(data: DataFrame, type: str = 'response') -> ndarray

Predict for new data.

type='response': the majority class of the terminal node, or its mean for a numeric outcome. 'prob': class shares in the terminal node, one column per class. 'node': the terminal node itself.

KDensityResult ¶

Bases: ResultProtocolMixin

Results from kernel density estimation.

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(3)
>>> df = pd.DataFrame({"income": rng.normal(10.0, 2.0, 500)})
>>> res = sp.kdensity(df, x="income", kernel="gaussian")
>>> type(res).__name__
'KDensityResult'
>>> bool(res.bandwidth > 0)
True
>>> bool((res.density >= 0).all())
True

plot ¶

plot(ax: Any = None, hist: bool = False, rug: bool = False, **kwargs: Any) -> Any

Plot the kernel density estimate.

LPolyResult ¶

Bases: ResultProtocolMixin

Results from local polynomial regression.

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(5)
>>> x = rng.uniform(0, 10, 300)
>>> y = np.sin(x) + rng.normal(0, 0.3, 300)
>>> df = pd.DataFrame({"y": y, "x": x})
>>> res = sp.lpoly(df, y="y", x="x", degree=1)
>>> type(res).__name__
'LPolyResult'
>>> res.degree
1
>>> bool(res.bandwidth > 0)
True

plot ¶

plot(ax: Any = None, scatter: bool = True, ci: bool = True, **kwargs: Any) -> Any

Plot the local polynomial fit with optional scatter and CI.

LProbustPoint `dataclass` ¶

One evaluation point of :func:lprobust_at_point.

Attributes:

Name	Type	Description
`eval_point`	`float`	The evaluation point `x0`.
`h, b`	`float`	Main and bias bandwidths.
`n_eff`	`int`	Observations with positive kernel weight in the `h` window or the `b` window (the rows that enter either fit).
`tau_us`	`float`	Conventional estimate of `m^(deriv)(x0)`.
`tau_bc`	`float`	Bias-corrected estimate.
`se_us`	`float`	Standard error of `tau_us`.
`se_rb`	`float`	Robust standard error of `tau_bc`.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> x = rng.uniform(-1, 1, 500)
>>> y = np.sin(2 * x) + rng.normal(0, 0.3, 500)
>>> point = sp.lprobust_at_point(x=x, y=y, eval_point=0.0, h=0.5)
>>> type(point).__name__
'LProbustPoint'
>>> point.n_eff                      # observations inside the bandwidth
248

The robust bias-corrected standard error is wider than the conventional one, because it also carries the variance of the estimated bias:

>>> bool(point.se_rb > point.se_us)
True

bias `property` ¶

bias: float

Estimated bias of the conventional estimate, tau_us - tau_bc.

ctree ¶

ctree(data: DataFrame, y: str, covariates: Sequence[str], *, task: str = 'auto', alpha: float = 0.05, test: str = 'bonferroni', min_split: int = 20, min_bucket: int = 7, min_prob: float = 0.01, max_depth: Optional[int] = None, split_stat: str = 'quadratic') -> CTreeResult

Conditional inference tree (Hothorn, Hornik and Zeileis 2006).

In each node every covariate is tested for association with the outcome by a permutation test in its asymptotic chi-squared form. The covariate with the smallest p-value is split on, at the cut point with the largest two-sample statistic. A node whose smallest adjusted p-value is not below alpha becomes a leaf, so the size of the tree is set by the significance level and there is no pruning step.

Parameters:

Name	Type	Description	Default
`data`	`DataFrame`		required
`y`	`str`	Outcome column.	required
`covariates`	`sequence of str`	Numeric, boolean and ordered-categorical columns are split at a cut point. String and unordered-categorical columns are split into two sets of levels (at most 12 levels).	required
`task`	`('auto', 'classification', 'regression')`	`'auto'`: classification for a string, categorical or boolean outcome, regression for a numeric one. A numeric 0/1 outcome gives the same tree either way; classification predicts the majority class, regression the share.	`'auto'`
`alpha`	`float`	Significance level a split must reach (`1 - mincriterion` in `party`).	`0.05`
`test`	`('bonferroni', 'univariate')`	`'bonferroni'`: the p-value is adjusted to `1 - (1 - p)^m` for the `m` covariates that vary in the node. `'univariate'`: no adjustment.	`'bonferroni'`
`min_split`	`int`	A node with fewer observations is not split.	`20`
`min_bucket`	`int`		`7`
`min_prob`	`float`	A child must hold at least `max(min_bucket, min_prob * n)` observations, `n` being the size of the node.	`0.01`
`max_depth`	`int`	Depth at which growth stops; the root has depth 0.	`None`
`split_stat`	`('quadratic', 'maximum')`	Statistic that ranks cut points: the quadratic form (the default of `partykit`) or the largest standardised component (the paper). They coincide for a numeric or a two-class outcome.	`'quadratic'`

Returns:

Type	Description
`CTreeResult`

Notes

Variable selection uses the quadratic-form statistic c = (t - mu)' Sigma^+ (t - mu) of the linear statistic t = sum_i g(x_i) h(y_i)', where h is the outcome itself or its class indicators and g is the covariate itself or its level indicators. Its conditional mean and covariance under permutation are those of eq. 2 of the paper, and c is referred to a chi-squared distribution with rank(Sigma) degrees of freedom. The maximum-type test statistic and Monte Carlo p-values of the R packages are not implemented.

The p-values are those of the tests that grew the tree. They are not valid for inference on the fitted tree, whose splits were chosen by looking at the same data.

Rows with a missing outcome or covariate are refused; surrogate splits are not implemented. There are no case weights.

With the same data and settings the tree, the test statistics and the p-values equal those of partykit::ctree (and of party::ctree, the implementation used by Varian 2014).

Examples:

>>> import numpy as np, pandas as pd, statspai as sp
>>> rng = np.random.default_rng(1)
>>> df = pd.DataFrame({"age": rng.uniform(1, 70, 500),
...                    "noise": rng.normal(size=500)})
>>> df["survived"] = (rng.uniform(size=500)
...                   < np.where(df["age"] < 20, 0.7, 0.35)).astype(int)
>>> fit = sp.ctree(df, y="survived", covariates=["age", "noise"],
...                task="classification")
>>> fit.nodes.loc[1, "variable"]
'age'
>>> fit.predict(df.head(2), type="prob").shape
(2, 2)

References

[@hothorn2006unbiased], [@varian2014big]

lpbwselect_ce_rot ¶

lpbwselect_ce_rot(x: ndarray, y: ndarray, eval_point: float, *, kernel: str = 'epanechnikov', p: int = 1, deriv: int = 0, n_neighbors: int = 3, bwcheck: int = 21, bwregul: float = 1.0) -> Dict[str, float]

Coverage-error-optimal rule of thumb, as bwselect('ce-rot').

The :func:lpbwselect_mse_dpi bandwidths shrunk by the sample-size factors of Calonico, Cattaneo and Farrell (2019, Section 2.6): h * n^(-p/((2p+3)(p+3))) and b * n^(-(q+2)/((2q+5)(q+3))) with q = p + 1 (b estimates m^(p+1) with an order-q fit, an even case). n counts the rows actually used. The floor bwcheck acts before the rescaling, not after.

Parameters:

Name	Type	Description	Default
`x`	`array - like`	Running variable and outcome; rows with a non-finite value are dropped.	required
`y`	`array - like`	Running variable and outcome; rows with a non-finite value are dropped.	required
`eval_point`	`float`	Evaluation point.	required
`kernel`	`('epanechnikov', 'triangular', 'uniform', 'gaussian')`	Kernel function. Default epanechnikov.	`'epanechnikov'`
`p`	`int`	Polynomial order. Default 1.	`1`
`deriv`	`int`	Derivative order; `p - deriv` must be odd. Default 0.	`0`
`n_neighbors`	`int`	Nearest neighbours in the variance estimates. Default 3.	`3`
`bwcheck`	`int`	Minimum number of observations in every window. Default 21.	`21`
`bwregul`	`float`	Scale of the regularization term. Default 1.	`1.0`

Returns:

Type	Description
`dict`	`h` and `b`.

Raises:

Type	Description
`MethodIncompatibility`	For `p - deriv` even or invalid options.
`DataInsufficient`	When fewer than `bwcheck` rows are available.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> x = rng.uniform(0, 2, 400)
>>> y = np.sin(2 * x) + rng.normal(0, 0.3, 400)
>>> ce = sp.lpbwselect_ce_rot(x, y, eval_point=1.0)
>>> mse = sp.lpbwselect_mse_dpi(x, y, eval_point=1.0)
>>> bool(ce["h"] < mse["h"])
True

References

[@calonico2019nprobust], [@calonico2018effect]

lpbwselect_imse_dpi ¶

lpbwselect_imse_dpi(x: ndarray, y: ndarray, *, kernel: str = 'epanechnikov', p: int = 1, deriv: int = 0, imsegrid: int = 30, n_neighbors: int = 3, bwcheck: int = 21, bwregul: float = 1.0) -> Dict[str, float]

IMSE-optimal direct-plug-in bandwidths, as bwselect('imse-dpi').

On imsegrid equally spaced points from min(x) to max(x):

b puts the grid averages of the stage-2 variance constant and unregularized squared bias of :func:lpbwselect_mse_dpi into the stage-2 closed form;
h does the same with stage 3, where the bias derivative at each grid point is estimated at that point's own mse-dpi b (regularized by bwregul) and the squared bias is unregularized.

No cap or floor is applied to the two averages; bwcheck acts on the pointwise pilots.

Parameters:

Name	Type	Description	Default
`x`	`array - like`	Running variable and outcome; rows with a non-finite value are dropped.	required
`y`	`array - like`	Running variable and outcome; rows with a non-finite value are dropped.	required
`kernel`	`('epanechnikov', 'triangular', 'uniform', 'gaussian')`	Kernel function. Default epanechnikov.	`'epanechnikov'`
`p`	`int`	Polynomial order. Default 1.	`1`
`deriv`	`int`	Derivative order; `p - deriv` must be odd. Default 0.	`0`
`imsegrid`	`int`	Number of grid points, at least 2. Default 30.	`30`
`n_neighbors`	`int`	Nearest neighbours in the variance estimates. Default 3.	`3`
`bwcheck`	`int`	Minimum number of observations in every pilot window. Default 21.	`21`
`bwregul`	`float`	Scale of the regularization term of the pointwise `b`. Default 1.	`1.0`

Returns:

Type	Description
`dict`	`h` and `b`.

Raises:

Type	Description
`MethodIncompatibility`	For `p - deriv` even or invalid options.
`DataInsufficient`	When fewer than `bwcheck` rows are available.
`NumericalInstability`	When a pilot fit at a grid point is singular.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> x = rng.uniform(0, 2, 400)
>>> y = np.sin(2 * x) + rng.normal(0, 0.3, 400)
>>> bw = sp.lpbwselect_imse_dpi(x, y)
>>> bool(bw["h"] > 0 and bw["b"] > 0)
True

References

[@calonico2019nprobust]

lpbwselect_imse_rot ¶

lpbwselect_imse_rot(x: ndarray, y: ndarray, *, kernel: str = 'epanechnikov', p: int = 1, deriv: int = 0) -> Dict[str, float]

IMSE-optimal rule-of-thumb bandwidth, as bwselect('imse-rot').

The pointwise mse-rot constants V and B^2 are averaged over the x quantiles at probabilities 0.05, 0.075, ..., 0.95 and put into the same closed form. The result is global, so there is no evaluation point.

Parameters:

Name	Type	Description	Default
`x`	`array - like`	Running variable and outcome; rows with a non-finite value are dropped.	required
`y`	`array - like`	Running variable and outcome; rows with a non-finite value are dropped.	required
`kernel`	`('epanechnikov', 'triangular', 'uniform')`	Kernel function (compact kernels only). Default epanechnikov.	`'epanechnikov'`
`p`	`int`	Polynomial order. Default 1.	`1`
`deriv`	`int`	Derivative order; `p - deriv` must be odd. Default 0.	`0`

Returns:

Type	Description
`dict`	`h`. The bias bandwidth is not reported, as for `mse-rot`.

Raises:

Type	Description
`MethodIncompatibility`	For the gaussian kernel, `p - deriv` even or invalid options.
`DataInsufficient`	When the global polynomial has too few observations.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> x = rng.uniform(0, 2, 400)
>>> y = np.sin(2 * x) + rng.normal(0, 0.3, 400)
>>> bool(sp.lpbwselect_imse_rot(x, y)["h"] > 0)
True

References

[@calonico2019nprobust]

lpbwselect_mse_dpi ¶

lpbwselect_mse_dpi(x: ndarray, y: ndarray, eval_point: float, *, kernel: str = 'epanechnikov', p: int = 1, deriv: int = 0, n_neighbors: int = 3, bwcheck: int = 21, bwregul: float = 1.0) -> Dict[str, float]

MSE-optimal direct-plug-in bandwidths, as bwselect('mse-dpi').

Three closed-form stages (Calonico, Cattaneo and Farrell 2019, Section 2.3, case p - deriv odd), all with variance and bias constants evaluated at a normal-reference pilot bandwidth:

a bandwidth for m^(p+2) (order-p+2 fit), whose bias uses an order-p+3 fit over the whole sample;
b, the MSE-optimal bandwidth for m^(p+1) with an order-p+1 fit, its bias derivative estimated at stage 1;
h, the MSE-optimal bandwidth for m^(deriv) with an order-p fit, its bias derivative estimated at b.

Stages 2 and 3 add bwregul times three times the variance of the estimated bias to the squared bias. Every bandwidth is capped at the largest distance from eval_point to the data and floored at the distance to the bwcheck-th nearest observation.

Parameters:

Name	Type	Description	Default
`x`	`array - like`	Running variable and outcome; rows with a non-finite value are dropped.	required
`y`	`array - like`	Running variable and outcome; rows with a non-finite value are dropped.	required
`eval_point`	`float`	Evaluation point.	required
`kernel`	`('epanechnikov', 'triangular', 'uniform', 'gaussian')`	Kernel function. Default epanechnikov.	`'epanechnikov'`
`p`	`int`	Order of the local polynomial for the point estimate. Default 1.	`1`
`deriv`	`int`	Derivative order; `p - deriv` must be odd. Default 0.	`0`
`n_neighbors`	`int`	Nearest neighbours in the variance estimates. Default 3.	`3`
`bwcheck`	`int`	Minimum number of observations every window must contain. Default 21. `1` imposes no minimum.	`21`
`bwregul`	`float`	Scale of the regularization term; `0` removes it. Default 1.	`1.0`

Returns:

Type	Description
`dict`	`h` and `b`, plus the pieces of the closed forms in coefficient units: `V_h` / `B_h` (stage 3) and `V_b` / `B_b` (stage 2), where `bw^(2 order + 3) = V / (n (B + regularization))`.

Raises:

Type	Description
`MethodIncompatibility`	For `p - deriv` even or invalid options.
`DataInsufficient`	When fewer than `bwcheck` rows are available.
`NumericalInstability`	When a pilot fit is singular.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> x = rng.uniform(0, 2, 400)
>>> y = np.sin(2 * x) + rng.normal(0, 0.3, 400)
>>> bw = sp.lpbwselect_mse_dpi(x, y, eval_point=1.0)
>>> bool(0 < bw["h"] <= 1.0 + 1e-12)
True

References

[@calonico2019nprobust], [@calonico2014robust]

lpbwselect_mse_rot ¶

lpbwselect_mse_rot(x: ndarray, y: ndarray, eval_point: float, *, kernel: str = 'epanechnikov', p: int = 1, deriv: int = 0) -> Dict[str, float]

MSE-optimal rule-of-thumb bandwidth, as bwselect('mse-rot').

h = ((1 + 2 deriv) V / (2 (p + 1 - deriv) n B^2))^(1 / (2p + 3)) with B = C1 * m^(p+1)(x0) / (p+1)! and V = C2 * sigma^2 / f(x0): C1, C2 are the equivalent-kernel bias and variance constants on [0, 1]; the derivative comes from an order-p+3 kernel-weighted polynomial fit over the whole sample; sigma^2 is the residual variance of a global polynomial of degree p + 3; f is a uniform-kernel density estimate at the normal-reference pilot bandwidth. No cap or bwcheck floor is applied.

Parameters:

Name	Type	Description	Default
`x`	`array - like`	Running variable and outcome; rows with a non-finite value are dropped.	required
`y`	`array - like`	Running variable and outcome; rows with a non-finite value are dropped.	required
`eval_point`	`float`	Evaluation point.	required
`kernel`	`('epanechnikov', 'triangular', 'uniform')`	Kernel function (compact kernels only). Default epanechnikov.	`'epanechnikov'`
`p`	`int`	Polynomial order. Default 1.	`1`
`deriv`	`int`	Derivative order; `p - deriv` must be odd. Default 0.	`0`

Returns:

Type	Description
`dict`	`h` and the pieces `B`, `V`, `C1`, `C2`. The selector's bias bandwidth is not reported (it would need the `p - deriv` even branch); consumers that run `rho = 1` use `h` for both.

Raises:

Type	Description
`MethodIncompatibility`	For the gaussian kernel, `p - deriv` even or invalid options.
`DataInsufficient`	When the global polynomial has too few observations.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> x = rng.uniform(0, 2, 400)
>>> y = np.sin(2 * x) + rng.normal(0, 0.3, 400)
>>> bool(sp.lpbwselect_mse_rot(x, y, eval_point=1.0)["h"] > 0)
True

References

[@calonico2019nprobust]

lprobust_at_point ¶

lprobust_at_point(x: ndarray, y: ndarray, eval_point: float, h: float, b: Optional[float] = None, *, kernel: str = 'epanechnikov', p: int = 1, deriv: int = 0, n_neighbors: int = 3) -> LProbustPoint

Local polynomial fit at eval_point with robust bias correction.

Estimates m^(deriv)(eval_point) with an order-p local polynomial at bandwidth h, removes the estimated leading bias with an order-p+1 fit at bandwidth b, and reports both standard errors of Calonico, Cattaneo and Farrell (2018). Rows with positive kernel weight in either window form the estimation sample; the nearest-neighbour variance estimate is computed within that sample.

No minimum window size is imposed: the bandwidths are used as given.

Parameters:

Name	Type	Description	Default
`x`	`array - like`	Running variable (the regressor).	required
`y`	`array - like`	Outcome, same length as `x`. Rows where either is non-finite are dropped.	required
`eval_point`	`float`	Point at which the regression function (or its derivative) is estimated.	required
`h`	`float`	Main bandwidth, a positive finite number.	required
`b`	`float`	Bias bandwidth, for the order-`p+1` fit that estimates the leading bias. Defaults to `h` (`rho = 1`), the reference's behaviour when only `h` is supplied.	`None`
`kernel`	`('epanechnikov', 'triangular', 'uniform', 'gaussian')`	Kernel function. Default epanechnikov, as in the reference.	`'epanechnikov'`
`p`	`int`	Order of the local polynomial for the point estimate. Default 1.	`1`
`deriv`	`int`	Order of the derivative estimated, `0 <= deriv <= p`. Default 0.	`0`
`n_neighbors`	`int`	Minimum number of nearest neighbours in the variance estimate; ties with the last neighbour are all included. Default 3.	`3`

Returns:

Type	Description
`LProbustPoint`	Bandwidths, effective sample size, the two estimates and the two standard errors.

Raises:

Type	Description
`MethodIncompatibility`	For an invalid bandwidth, kernel or order, or when the window holds too few distinct values of `x` for the polynomial fits.
`DataInsufficient`	When `p + 2` or fewer observations fall inside the bandwidths.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(1)
>>> x = rng.uniform(0, 2, 800)
>>> y = 1 + x - 0.5 * x**2 + rng.normal(0, 0.2, 800)
>>> fit = sp.lprobust_at_point(x, y, eval_point=0.0, h=0.6)
>>> bool(abs(fit.tau_bc - 1.0) < 3 * fit.se_rb)
True

References

[@calonico2018effect], [@calonico2019nprobust]

statspai.nonparametric¶

nonparametric ¶

CTreeResult dataclass ¶

predict ¶

KDensityResult ¶

plot ¶

LPolyResult ¶

plot ¶

LProbustPoint dataclass ¶

bias property ¶

ctree ¶

lpbwselect_ce_rot ¶

lpbwselect_imse_dpi ¶

lpbwselect_imse_rot ¶

lpbwselect_mse_dpi ¶

lpbwselect_mse_rot ¶

lprobust_at_point ¶

`statspai.nonparametric`¶

CTreeResult `dataclass` ¶

LProbustPoint `dataclass` ¶

bias `property` ¶