andrey.pnl

andrey.pnl(data, *, alpha=None)

Test the direction of a variable pair with the post-nonlinear (PNL) model.

PNL [1] [2] fits y = f2(f1(x) + e) in each direction. An HSIC test checks whether the estimated noise e is independent of the proposed cause. The two p-values determine the edge, and the two test statistics decide a tie.

Parameters:

• data (ArrayLike) – The two variables, one per column: column 0 is x and column 1 is y. Coerced to float64. A DataFrame’s column names become the structure’s labels.

• alpha (float | None) – None chooses the direction with the larger p-value, or with the smaller test statistic when the p-values are equal. A significance level such as 0.05 draws an edge only when exactly one p-value exceeds it. See Notes.

Returns:

Contains a two-node GraphStructure of kind dag with the chosen edge, or no edge. metadata["algorithm"] is "PNL"; pval_forward is the independence p-value for x -> y and pval_backward is the p-value for y -> x; stat_forward and stat_backward are the HSIC statistics of the same two tests, each divided by its mean under independence.

Return type:

StructureOutput

Raises:

ValueError – If data does not have 2 columns, or if alpha is not in (0, 1).

Notes

With alpha=None, the edge points in the direction with the larger p-value. On a large sample, or a pair that fits neither direction, both p-values can be zero; equal p-values leave the choice to the HSIC statistics, each divided by its mean under independence, and the edge points in the direction with the smaller one, the weaker dependence between the proposed cause and its estimated noise. Equal statistics give no edge. The rule therefore produces an edge even for unrelated variables or a pair whose direction cannot be identified, such as a linear Gaussian pair. Mooij et al. (2016) [4] compare p-values to avoid choosing a significance level and found the comparison among the best methods in their benchmark of additive-noise models.

With a significance level, only the p-values decide; a direction passes when its p-value exceeds alpha:

• only x -> y passes: the edge x -> y;

• only y -> x passes: the edge y -> x;

• both pass: no edge, since neither direction is rejected;

• neither passes: no edge, since both directions are rejected.

Zhang and Hyvarinen (2010) [2] use a level of 0.01 for PNL; Hoyer et al. (2008) [3] use 0.02 for additive-noise models. The result depends on alpha. Estimated noise can appear dependent on the true cause, so a standard level can reject the true direction; no theory establishes the right level (Mooij et al., 2016). Both p-values and both statistics remain in metadata.

PNL fits each transform as a sum of smooth basis functions, so the fit is deterministic but can miss a sharply varying nonlinearity.

The kernel matrices of the two HSIC tests run on a CUDA GPU, in float64, under the cuda backend, or under auto from ANDREY_PNL_GPU_THRESHOLD rows squared once the GPU is running in the process; the fits stay on the CPU, and the p-values then match the CPU’s to rounding.

On the CPU, num_workers above 1 (andrey.config or ANDREY_NUM_WORKERS) fits the two directions on two threads, each limited to one BLAS thread; with one worker the BLAS threads are not limited. On a machine with several cores the fits then round differently, so the p-values and statistics with one worker and with several can differ in their last digits, and the direction can differ when the two p-values nearly tie. Each thread holds three n x n float64 matrices; when the threads asked for would hold more than 4 GiB, fewer run, with a PerformanceWarning.

Limitations. Measured on simulated cause-effect pairs with a post-nonlinear mechanism, a linear inner transform and an invertible nonlinear outer one, and Laplace noise, from 200 to 5,000 rows. It named the true direction on three pairs in four at 200 rows and on every pair from 1,000 rows. The median fit took 0.8 s at 5,000 rows, on 4 cores. Each HSIC test builds kernel matrices with one row and one column per sample, so time and memory grow with the square of the row count: each worker thread holds three such float64 matrices, 600 MB at 5,000 rows.

References

[1]

Zhang, K., and Hyvarinen, A. (2009). On the identifiability of the post-nonlinear causal model. UAI, 647-655.

[2] (1,2)

Zhang, K., and Hyvarinen, A. (2010). Distinguishing causes from effects using nonlinear acyclic causal models. JMLR Workshop and Conference Proceedings 6, 157-164.

[3]

Hoyer, P. O., Janzing, D., Mooij, J. M., Peters, J., and Scholkopf, B. (2008). Nonlinear causal discovery with additive noise models. NeurIPS 21, 689-696.

[4]

Mooij, J. M., Peters, J., Janzing, D., Zscheischler, J., and Scholkopf, B. (2016). Distinguishing cause from effect using observational data: methods and benchmarks. JMLR 17(32), 1-102.

Examples

>>> import numpy as np
>>> from andrey import pnl
>>> rng = np.random.default_rng(2)
>>> x = rng.uniform(-1, 1, 300)
>>> y = np.tanh(x**3 + 0.3 * rng.uniform(-1, 1, 300))
>>> out = pnl(np.column_stack([x, y]))
>>> out.structure.oriented_edges()
[(0, 1, 'directed')]
>>> out.metadata["pval_forward"] > 0.05 >= out.metadata["pval_backward"]
True