diff --git a/benchmark/benchmark.cpp b/benchmark/benchmark.cpp
index 281983d..a1251b4 100644
--- a/benchmark/benchmark.cpp
+++ b/benchmark/benchmark.cpp
@@ -10,6 +10,9 @@ TEST_CASE("continuous KDE performance", "[!benchmark]")
unbounded.array() += 0.35 * (4.0 * unbounded.array()).sin();
Eigen::VectorXd bounded =
(1.0 / (1.0 + (-unbounded.array()).exp())).matrix();
+ Eigen::VectorXd probabilities = Eigen::VectorXd::LinSpaced(
+ unbounded.size(), 0.5 / unbounded.size(), 1.0 - 0.5 / unbounded.size());
+ Eigen::VectorXd one_sided = (-(1.0 - probabilities.array()).log());
Eigen::VectorXd weights = Eigen::VectorXd::LinSpaced(10000, 0.5, 1.5);
BENCHMARK("fit/end-to-end/bounded/tll/n=10000")
@@ -19,6 +22,20 @@ TEST_CASE("continuous KDE performance", "[!benchmark]")
return fit.get_loglik();
};
+ BENCHMARK("fit/end-to-end/one-sided-finite-expert/tll/n=10000")
+ {
+ kde1d::Kde1d fit(0.0, NAN, "continuous", 1.0, NAN, 2, 400);
+ fit.fit(one_sided);
+ return fit.get_loglik();
+ };
+
+ BENCHMARK("fit/end-to-end/two-sided-finite-experts/tll/n=10000")
+ {
+ kde1d::Kde1d fit(0.0, 1.0, "continuous", 1.0, NAN, 2, 400);
+ fit.fit(probabilities);
+ return fit.get_loglik();
+ };
+
BENCHMARK("fit/fixed-bandwidth/unbounded/tll/n=10000")
{
kde1d::Kde1d fit(NAN, NAN, "continuous", 1.0, 0.3, 2, 400);
diff --git a/cmake/buildTargets.cmake b/cmake/buildTargets.cmake
index eff6110..00a69c6 100755
--- a/cmake/buildTargets.cmake
+++ b/cmake/buildTargets.cmake
@@ -17,10 +17,6 @@ if(BUILD_BENCHMARKS)
add_subdirectory(benchmark)
endif(BUILD_BENCHMARKS)
-if(BUILD_DIAGNOSTICS)
- add_subdirectory(diagnostic)
-endif(BUILD_DIAGNOSTICS)
-
# Related to exports for linux/mac and code coverage
####
# Installation
diff --git a/cmake/options.cmake b/cmake/options.cmake
index 5434743..da71e6c 100755
--- a/cmake/options.cmake
+++ b/cmake/options.cmake
@@ -9,5 +9,4 @@ option(WARNINGS_AS_ERRORS "Compiler warnings as errors" "OFF")
option(OPT_ASAN "Use adress sanitizer (debug)" "ON")
option(BUILD_TESTING "Build tests." "ON")
option(BUILD_BENCHMARKS "Build benchmarks." "OFF")
-option(BUILD_DIAGNOSTICS "Build numerical diagnostics." "OFF")
option(CODE_COVERAGE "Code coverage." "OFF")
diff --git a/cmake/printInfo.cmake b/cmake/printInfo.cmake
index 0b2c2ac..254fa14 100755
--- a/cmake/printInfo.cmake
+++ b/cmake/printInfo.cmake
@@ -14,6 +14,5 @@ message( STATUS "CMAKE_CXX_FLAGS_RELEASE= ${CMAKE_CXX_FLAGS_RELEASE}")
message( STATUS )
message( STATUS "BUILD_TESTING= ${BUILD_TESTING}")
message( STATUS "BUILD_BENCHMARKS= ${BUILD_BENCHMARKS}")
-message( STATUS "BUILD_DIAGNOSTICS= ${BUILD_DIAGNOSTICS}")
message( STATUS "CODE_COVERAGE= ${CODE_COVERAGE}")
message( STATUS )
diff --git a/diagnostic/CMakeLists.txt b/diagnostic/CMakeLists.txt
deleted file mode 100644
index b27293b..0000000
--- a/diagnostic/CMakeLists.txt
+++ /dev/null
@@ -1,6 +0,0 @@
-add_executable(kde1d-boundary-diagnostic boundary_diagnostic.cpp)
-target_link_libraries(kde1d-boundary-diagnostic kde1d)
-set_target_properties(
- kde1d-boundary-diagnostic
- PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}/bin
-)
diff --git a/diagnostic/PR_COMMENT.md b/diagnostic/PR_COMMENT.md
new file mode 100644
index 0000000..a2d0fdd
--- /dev/null
+++ b/diagnostic/PR_COMMENT.md
@@ -0,0 +1,240 @@
+## Boundary-repair estimator
+
+This PR changes continuous fits with known finite support in two steps.
+
+For a one-sided support, write $d(x)$ for distance from the finite endpoint,
+let $s$ be the weighted median observed distance, and set
+$\epsilon=10^{-5}s$. The transformed bulk estimator uses an endpoint-anchored,
+scaled member of the Box-Cox family with power parameter $\lambda=1/4$:
+
+$$
+T\{d(x)\}
+=4[\{\frac{d(x)+\epsilon}{s}\}^{1/4}
+-\{\frac{\epsilon}{s}\}^{1/4}],
+\qquad
+|T'(x)|
+=\frac{1}{s}\{\frac{d(x)+\epsilon}{s}\}^{-3/4}.
+$$
+
+This Box-Cox transform replaces the previous log transform and makes the
+construction scale equivariant. Two-sided fits retain the regularized probit
+bulk transform.
+
+Each finite endpoint is then classified independently from its first
+$k\simeq2\sqrt{n_e}$ order statistics. In the unweighted case,
+
+$$
+\widehat\beta
+=\frac{k}{\sum_{i=1}^k
+\log\{d_{(k+1)}/d_{(i)}\}}.
+$$
+
+Endpoints consistent with a finite nonzero limiting density receive a boundary
+expert. Exploding, vanishing, ambiguous, and degenerate endpoints retain the
+bulk. The expert is a Gaussian local-linear equivalent-kernel estimate with an
+automatically selected bandwidth on the original scale. The selector uses all
+positive-weight observations on bounded support and the closest 75% by row
+count on one-sided support. The two required Gaussian convolutions are
+evaluated on a shared 256-bin FFT grid.
+
+Let $F_B$ be the normalized bulk CDF,
+
+$$
+q=\min(1/4,n_e^{-1/2}),
+\qquad
+g(p)=\{1-\min(1,p/q)^2\}^2.
+$$
+
+Selected lower and upper experts receive the biweight fusion weights
+$w_L(x)=g\{F_B(x)\}$ and $w_U(x)=g\{1-F_B(x)\}$, giving
+
+$$
+\widehat f(x)
+=w_L(x)\widehat f_L(x)
++\{1-w_L(x)-w_U(x)\}\widehat f_B(x)
++w_U(x)\widehat f_U(x).
+$$
+
+The result is clamped at zero and normalized. The same boundary expert is used
+for bulk degrees 0, 1, and 2. EDF is a conditional influence-trace
+approximation that holds the classifier, selected bandwidths, fusion weights,
+clamps, and normalization fixed.
+
+## Validation design
+
+> **Before merging:** remove `diagnostic/`; it is retained in this branch only
+> to make the PR validation reproducible during review.
+
+The primary benchmark is the actual pre-PR package at parent commit `7faec53`
+with C++ commit `b16d514`: log transform for one-sided support and probit for
+two-sided support. The experiment used 200 paired replications for
+$n\in\{25,100,1000,2000\}$ and degrees 0, 1, and 2.
+
+- One-sided: Exponential, half-normal, Gamma$(0.75)$, Gamma$(2)$, and
+ Lomax$(2)$, with lower/upper reflections and scales $10^{-4},1,10^4$.
+- Two-sided: Uniform and Beta$(1,2)$, Beta$(2,1)$, Beta$(0.75,2)$,
+ Beta$(2,0.75)$, Beta$(0.75,0.75)$, and Beta$(2,2)$.
+- Global ISE and boundary ISE were integrated separately. The boundary region
+ is the outer 10% probability mass at each finite endpoint.
+
+The tables below use scale 1; reflected one-sided results agree to numerical
+precision. Entries are candidate/baseline ratios, so **values below 1 favor the
+PR**. Parentheses are paired Monte Carlo standard errors.
+
+### Global ISE ratio, averaged over families
+
+| support | degree | $n=25$ | $n=100$ | $n=1000$ | $n=2000$ |
+|---|---:|---:|---:|---:|---:|
+| one-sided | 0 | 1.810 (0.116) | 1.134 (0.049) | 0.586 (0.015) | 0.519 (0.011) |
+| one-sided | 1 | 0.789 (0.038) | 0.508 (0.011) | 0.566 (0.042) | 0.612 (0.037) |
+| one-sided | 2 | 0.550 (0.034) | 0.394 (0.015) | 0.231 (0.008) | 0.231 (0.008) |
+| two-sided | 0 | 0.946 (0.005) | 0.940 (0.005) | 0.943 (0.006) | 0.926 (0.006) |
+| two-sided | 1 | 0.856 (0.007) | 0.817 (0.007) | 0.846 (0.005) | 0.850 (0.004) |
+| two-sided | 2 | 0.881 (0.013) | 0.906 (0.006) | 0.943 (0.007) | 0.932 (0.007) |
+
+### Boundary ratio, averaged over families
+
+| support | degree | $n=25$ | $n=100$ | $n=1000$ | $n=2000$ |
+|---|---:|---:|---:|---:|---:|
+| one-sided | 0 | 3.149 (0.278) | 1.553 (0.097) | 0.563 (0.020) | 0.461 (0.014) |
+| one-sided | 1 | 0.931 (0.073) | 0.437 (0.017) | 0.519 (0.047) | 0.590 (0.042) |
+| one-sided | 2 | 0.612 (0.048) | 0.396 (0.024) | 0.190 (0.009) | 0.190 (0.008) |
+| two-sided | 0 | 0.918 (0.005) | 0.939 (0.005) | 0.955 (0.006) | 0.940 (0.006) |
+| two-sided | 1 | 0.817 (0.009) | 0.821 (0.007) | 0.878 (0.004) | 0.876 (0.004) |
+| two-sided | 2 | 0.877 (0.012) | 0.910 (0.005) | 0.958 (0.006) | 0.946 (0.006) |
+
+Against the unrepaired new bulk, averaged global ISE ratios range from
+0.83--0.94 for degree 0, 0.73--0.88 for degree 1, and 0.87--0.92 for degree 2
+on one-sided support. The corresponding bounded ratios are 0.93--0.95,
+0.82--0.86, and 0.88--0.94.
+
+
+Family-level global / boundary ISE ratios
+
+#### One-sided, degree 0
+
+| family | $n=25$ | $n=100$ | $n=1000$ | $n=2000$ |
+|---|---:|---:|---:|---:|
+| Exponential | 1.48 / 2.67 | 1.36 / 2.41 | 0.58 / 0.56 | 0.51 / 0.45 |
+| Gamma$(0.75)$ | 2.05 / 2.76 | 0.91 / 0.97 | 0.33 / 0.32 | 0.28 / 0.27 |
+| Gamma$(2)$ | 0.79 / 0.66 | 0.84 / 0.75 | 0.87 / 0.80 | 0.86 / 0.77 |
+| half-normal | 0.88 / 1.28 | 0.79 / 1.08 | 0.54 / 0.53 | 0.46 / 0.38 |
+| Lomax$(2)$ | 3.84 / 8.37 | 1.77 / 2.56 | 0.61 / 0.60 | 0.50 / 0.44 |
+
+#### One-sided, degree 1
+
+| family | $n=25$ | $n=100$ | $n=1000$ | $n=2000$ |
+|---|---:|---:|---:|---:|
+| Exponential | 0.70 / 0.83 | 0.47 / 0.43 | 0.37 / 0.25 | 0.35 / 0.24 |
+| Gamma$(0.75)$ | 0.72 / 0.88 | 0.33 / 0.30 | 0.32 / 0.30 | 0.36 / 0.34 |
+| Gamma$(2)$ | 0.74 / 0.59 | 0.84 / 0.73 | 0.93 / 1.03 | 1.00 / 1.21 |
+| half-normal | 0.61 / 0.65 | 0.50 / 0.40 | 0.47 / 0.29 | 0.44 / 0.26 |
+| Lomax$(2)$ | 1.17 / 1.71 | 0.40 / 0.32 | 0.74 / 0.72 | 0.91 / 0.90 |
+
+#### One-sided, degree 2
+
+| family | $n=25$ | $n=100$ | $n=1000$ | $n=2000$ |
+|---|---:|---:|---:|---:|
+| Exponential | 0.43 / 0.53 | 0.32 / 0.38 | 0.10 / 0.07 | 0.10 / 0.07 |
+| Gamma$(0.75)$ | 0.47 / 0.56 | 0.18 / 0.17 | 0.12 / 0.11 | 0.10 / 0.10 |
+| Gamma$(2)$ | 0.61 / 0.48 | 0.84 / 0.83 | 0.67 / 0.59 | 0.67 / 0.61 |
+| half-normal | 0.37 / 0.42 | 0.32 / 0.32 | 0.16 / 0.10 | 0.13 / 0.06 |
+| Lomax$(2)$ | 0.87 / 1.09 | 0.30 / 0.28 | 0.11 / 0.08 | 0.15 / 0.11 |
+
+#### Two-sided, degree 0
+
+| family | $n=25$ | $n=100$ | $n=1000$ | $n=2000$ |
+|---|---:|---:|---:|---:|
+| Beta$(0.75,0.75)$ | 1.01 / 1.02 | 1.03 / 1.04 | 1.05 / 1.05 | 1.06 / 1.06 |
+| Beta$(0.75,2)$ | 0.98 / 0.97 | 1.01 / 1.01 | 1.09 / 1.08 | 1.07 / 1.07 |
+| Beta$(1,2)$ | 0.90 / 0.89 | 0.84 / 0.85 | 0.82 / 0.85 | 0.77 / 0.80 |
+| Beta$(2,0.75)$ | 0.99 / 0.98 | 1.01 / 1.01 | 1.11 / 1.10 | 1.07 / 1.06 |
+| Beta$(2,1)$ | 0.94 / 0.92 | 0.87 / 0.88 | 0.76 / 0.80 | 0.78 / 0.81 |
+| Beta$(2,2)$ | 0.86 / 0.68 | 0.95 / 0.88 | 1.00 / 1.00 | 1.00 / 1.00 |
+| Uniform | 0.94 / 0.96 | 0.86 / 0.91 | 0.76 / 0.81 | 0.72 / 0.78 |
+
+#### Two-sided, degree 1
+
+| family | $n=25$ | $n=100$ | $n=1000$ | $n=2000$ |
+|---|---:|---:|---:|---:|
+| Beta$(0.75,0.75)$ | 0.90 / 0.90 | 0.86 / 0.88 | 0.95 / 0.95 | 0.98 / 0.98 |
+| Beta$(0.75,2)$ | 0.86 / 0.86 | 0.90 / 0.90 | 0.97 / 0.97 | 0.99 / 0.99 |
+| Beta$(1,2)$ | 0.79 / 0.74 | 0.71 / 0.73 | 0.72 / 0.75 | 0.69 / 0.72 |
+| Beta$(2,0.75)$ | 0.87 / 0.87 | 0.88 / 0.89 | 0.96 / 0.96 | 0.99 / 0.99 |
+| Beta$(2,1)$ | 0.84 / 0.81 | 0.72 / 0.75 | 0.69 / 0.74 | 0.71 / 0.73 |
+| Beta$(2,2)$ | 0.93 / 0.78 | 0.97 / 0.90 | 1.00 / 1.00 | 1.00 / 1.00 |
+| Uniform | 0.79 / 0.76 | 0.68 / 0.70 | 0.62 / 0.76 | 0.61 / 0.73 |
+
+#### Two-sided, degree 2
+
+| family | $n=25$ | $n=100$ | $n=1000$ | $n=2000$ |
+|---|---:|---:|---:|---:|
+| Beta$(0.75,0.75)$ | 0.90 / 0.91 | 0.94 / 0.95 | 1.04 / 1.04 | 1.07 / 1.07 |
+| Beta$(0.75,2)$ | 0.87 / 0.88 | 0.95 / 0.96 | 1.09 / 1.08 | 1.08 / 1.07 |
+| Beta$(1,2)$ | 0.74 / 0.77 | 0.82 / 0.83 | 0.81 / 0.83 | 0.77 / 0.79 |
+| Beta$(2,0.75)$ | 0.89 / 0.89 | 0.95 / 0.96 | 1.11 / 1.10 | 1.07 / 1.07 |
+| Beta$(2,1)$ | 0.89 / 0.89 | 0.80 / 0.82 | 0.76 / 0.80 | 0.77 / 0.81 |
+| Beta$(2,2)$ | 1.04 / 0.96 | 1.04 / 0.96 | 1.00 / 1.00 | 1.00 / 1.00 |
+| Uniform | 0.85 / 0.85 | 0.84 / 0.88 | 0.80 / 0.85 | 0.77 / 0.81 |
+
+
+
+## Bias, variance, classification, and scale
+
+The repair primarily reduces variance, with deliberate boundary bias. Against
+pre-PR at scale 1, aggregate global variance ratios are 0.38--0.62 one-sided
+for degree 0, 0.39--0.56 for degree 1, and 0.20--0.24 for degree 2. Bounded
+ratios are 0.96--0.98, 0.84--0.87, and 0.87--0.96, respectively. This matches
+the visual goal: substantially calmer one-sided endpoints and moderate
+smoothing on bounded support.
+
+Classifier behavior becomes appropriately conservative with sample size. For
+one-sided finite targets, selection rates at $n=100/2000$ are 58%/78%
+(Exponential), 75%/82% (half-normal), and 56%/76% (Lomax). False finite
+selection falls from 22% to 4% for Gamma$(0.75)$ and from 38% to 0% for
+Gamma$(2)$. For Uniform, either endpoint is selected in 90%/98% and both in
+51%/71% of samples. Beta$(2,2)$ falls from 4.5% both-endpoint selection at
+$n=100$ to 0% at $n\geq1000$.
+
+For Exponential samples at $n=100$, the 5%--95% range of the estimated value
+at the endpoint is 0.64--1.34 when the expert is selected, compared with
+0.01--49.6 when the classifier retains the raw bulk. Thus the expert delivers
+the intended visual stabilization when selected; remaining extreme endpoint
+values are primarily classifier false negatives.
+
+
+
+
+## Coarse runtime
+
+Median end-to-end R fit times for the default degree-2 estimator are below.
+
+| case | $n$ | pre-PR | candidate | slowdown |
+|---|---:|---:|---:|---:|
+| one active endpoint | 100 | 0.48 ms | 0.79 ms | 1.66x |
+| one active endpoint | 1,000 | 1.52 ms | 1.90 ms | 1.25x |
+| one active endpoint | 10,000 | 12.47 ms | 13.23 ms | 1.06x |
+| two active endpoints | 100 | 0.51 ms | 0.83 ms | 1.63x |
+| two active endpoints | 1,000 | 1.46 ms | 1.77 ms | 1.22x |
+| two active endpoints | 10,000 | 10.97 ms | 11.40 ms | 1.04x |
+| one-sided fallback | 100 / 1,000 / 10,000 | -- | -- | 1.08x / 1.03x / 1.00x |
+| two-sided fallback | 100 / 1,000 / 10,000 | -- | -- | 1.00x / 0.98x / 0.99x |
+
+## Less common paths and numerical checks
+
+- Public case weights are used consistently by the classifier, expert
+ estimand, bandwidth selector, effective sample size, and EDF approximation.
+ Tests cover weight rescaling and zero-weight-row invariance.
+- A manually supplied bandwidth controls the transformed bulk only. Boundary
+ experts retain their original-scale automatic bandwidth; the public
+ multiplier scales both. This unusual path has focused regression tests but
+ was not included in the full Monte Carlo sweep.
+- Setting `boundary_repair=false` recovers the transformed bulk and skips all
+ expert work.
+- Tests cover degrees 0--2, normalization, nonnegativity, monotone CDFs,
+ quantile round trips, support truncation, reflection, scales $10^{-4}$ to
+ $10^4$, weights, fixed bandwidths, and pure-bulk EDF parity.
+- The conditional EDF approximation differed from leave-one-out optimism by
+ about 0.08--0.30 EDF units for deterministic Uniform/Exponential examples;
+ the largest checked discrepancy was 0.84 for Beta$(1,2)$ at $n=25$.
diff --git a/diagnostic/README.md b/diagnostic/README.md
index d34d715..2d73909 100644
--- a/diagnostic/README.md
+++ b/diagnostic/README.md
@@ -1,31 +1,60 @@
-# Boundary diagnostics
+# Boundary-expert validation
-The diagnostic sweep evaluates deterministic samples from distributions with
-known densities and CDFs. It covers finite and one-sided supports, support
-scales from `1e-4` to `1e4`, polynomial degrees 0--2, grid sizes 100--1,000,
-and fixed or automatically selected bandwidths. Evaluation points approach the
-boundary on a logarithmic scale.
+This directory contains the retained validation path for the native estimator.
+The estimator itself, including the support transforms, endpoint classifier,
+boundary kernels, fusion, bandwidth handling, and EDF approximation, is
+documented in [`../docs/overview-continuous.dox`](../docs/overview-continuous.dox).
-Generate the raw CSV data and an optional summary and PDF with:
+## Final comparison
+
+`validate_native_boundary_experts.R` evaluates one installed implementation at
+a time. Run it with these method labels:
+
+| method | installed revision | estimator |
+|---|---|---|
+| `pre-pr` | parent `7faec53`, with C++ `b16d514` | production log transform on one-sided support and probit transform on two-sided support |
+| `bulk` | candidate | fourth-root/probit bulk with `boundary_repair = FALSE` |
+| `expert` | candidate | fourth-root/probit bulk with the selected boundary experts |
+
+Using the actual pre-PR revision is intentional: emulating its log transform in
+the candidate would miss other grid and numerical behavior of the production
+benchmark.
+
+The suite covers lower- and upper-bounded versions of five one-sided densities
+at scales `1e-4`, `1`, and `1e4`, seven two-sided densities, and sample sizes
+25, 100, 1,000, and 2,000 for bulk degrees zero, one, and two. It writes
+replication-level global and boundary ISE, EDF, and log-likelihood, plus
+integrated squared bias and variance summaries. Paired seeds are shared across
+methods and degrees.
+
+From the R-package root, install each requested revision into a separate
+library and run, for example:
```sh
-diagnostic/run_boundary_diagnostics.sh > boundary-diagnostics.csv
-Rscript diagnostic/plot_boundary_diagnostics.R boundary-diagnostics.csv
+Rscript inst/include/kde1d-cpp/diagnostic/validate_native_boundary_experts.R \
+ pre-pr diagnostic-results/pre-pr.csv 100 /path/to/pre-pr-library
+Rscript inst/include/kde1d-cpp/diagnostic/validate_native_boundary_experts.R \
+ bulk diagnostic-results/bulk.csv 100 /path/to/candidate-library
+Rscript inst/include/kde1d-cpp/diagnostic/validate_native_boundary_experts.R \
+ expert diagnostic-results/expert.csv 100 /path/to/candidate-library
```
-The executable is deliberately not registered with CTest. Its output is for
-comparing numerical behavior between revisions; stable invariants discovered
-through the sweep should be promoted to focused regression tests.
+Generated CSVs and figures are not versioned. The final aggregate comparison
+is generated with:
-The summary reports four complementary diagnostics:
+```sh
+Rscript inst/include/kde1d-cpp/diagnostic/summarize_native_boundary_validation.R \
+ diagnostic-results/pre-pr.csv diagnostic-results/bulk.csv \
+ diagnostic-results/expert.csv diagnostic-results/ratios
+```
-- PDF and CDF error against known distributions near the boundary;
-- equivariance when the same data are expressed at different scales;
-- agreement between left-bounded fits and reflected right-bounded fits;
-- sensitivity to the interpolation grid size.
+`benchmark_native_boundary_experts.R` performs coarse end-to-end R timings.
+`extract_native_boundary_examples.R` and
+`plot_native_boundary_examples.R` recreate the representative density montage
+using classifier-stratified median-ISE realizations selected by the summary
+script.
-The metrics deliberately remain separate. Relative density error is not useful
-where the true density vanishes, while CDF error can remain small despite a
-severely unstable pointwise density. Candidate repairs should therefore be
-checked against both finite nonzero boundary densities and densities that tend
-to zero.
+Final results are recorded in `NATIVE_RESULTS.md`; `PR_COMMENT.md` contains the
+copy-ready methodological description and report. Superseded R prototypes,
+candidate estimators, tuning outputs, and the standalone C++ diagnostic were
+removed; stable numerical properties are covered by unit and invariant tests.
diff --git a/diagnostic/benchmark_native_boundary_experts.R b/diagnostic/benchmark_native_boundary_experts.R
new file mode 100644
index 0000000..6458303
--- /dev/null
+++ b/diagnostic/benchmark_native_boundary_experts.R
@@ -0,0 +1,59 @@
+arguments <- commandArgs(trailingOnly = TRUE)
+if (length(arguments) != 3L) {
+ stop("usage: script method output.csv library-path")
+}
+.libPaths(c(arguments[[3L]], .libPaths()))
+library(kde1d)
+
+method <- arguments[[1L]]
+if (!method %in% c("pre-pr", "bulk", "expert")) {
+ stop("method must be one of: pre-pr, bulk, expert")
+}
+
+fit_kde1d <- function(observations, xmin, xmax) {
+ fit_arguments <- list(x = observations, xmin = xmin, xmax = xmax, deg = 2)
+ if (method != "pre-pr") {
+ fit_arguments$boundary_repair <- method == "expert"
+ }
+ do.call(kde1d, fit_arguments)
+}
+
+make_observations <- function(case, n) {
+ probabilities <- seq(0.5 / n, 1 - 0.5 / n, length.out = n)
+ switch(
+ case,
+ one_active = -log1p(-probabilities),
+ one_fallback = sqrt(probabilities),
+ two_active = probabilities,
+ two_fallback = qbeta(probabilities, 2, 2)
+ )
+}
+
+rows <- list()
+for (case in c("one_active", "one_fallback", "two_active", "two_fallback")) {
+ for (n in c(100L, 1000L, 10000L)) {
+ observations <- make_observations(case, n)
+ one_sided <- startsWith(case, "one_")
+ xmin <- 0
+ xmax <- if (one_sided) NaN else 1
+ iterations <- if (n == 100L) 500L else if (n == 1000L) 200L else 30L
+ fit_kde1d(observations, xmin, xmax)
+ for (batch in seq_len(7L)) {
+ elapsed <- system.time({
+ for (iteration in seq_len(iterations)) {
+ fit_kde1d(observations, xmin, xmax)
+ }
+ })[["elapsed"]]
+ rows[[length(rows) + 1L]] <- data.frame(
+ method = method,
+ case = case,
+ sample_size = n,
+ batch = batch,
+ milliseconds = 1000 * elapsed / iterations
+ )
+ }
+ }
+}
+
+dir.create(dirname(arguments[[2L]]), recursive = TRUE, showWarnings = FALSE)
+write.csv(do.call(rbind, rows), arguments[[2L]], row.names = FALSE)
diff --git a/diagnostic/boundary-examples.png b/diagnostic/boundary-examples.png
new file mode 100644
index 0000000..281cf42
Binary files /dev/null and b/diagnostic/boundary-examples.png differ
diff --git a/diagnostic/boundary_diagnostic.cpp b/diagnostic/boundary_diagnostic.cpp
deleted file mode 100644
index 27cc53a..0000000
--- a/diagnostic/boundary_diagnostic.cpp
+++ /dev/null
@@ -1,155 +0,0 @@
-#include
-#include
-#include
-#include
-#include
-#include
-#include
-#include
-
-namespace {
-
-struct Scenario
-{
- std::string name;
- Eigen::VectorXd observations;
- double xmin;
- double xmax;
- double scale;
-};
-
-Eigen::VectorXd
-probabilities(Eigen::Index size)
-{
- return Eigen::VectorXd::LinSpaced(
- size, 0.5 / static_cast(size),
- 1.0 - 0.5 / static_cast(size));
-}
-
-Scenario
-make_scenario(const std::string& name, Eigen::Index size, double scale)
-{
- Eigen::VectorXd probability = probabilities(size);
- if (name == "uniform") {
- return { name, probability * scale, 0.0, scale, scale };
- } else if (name == "beta-2-1") {
- return { name, probability.array().sqrt() * scale, 0.0, scale, scale };
- } else if (name == "beta-1-2") {
- return { name,
- (1.0 - (1.0 - probability.array()).sqrt()) * scale,
- 0.0,
- scale,
- scale };
- } else if (name == "exponential-left") {
- return { name,
- (-(1.0 - probability.array()).log()) * scale,
- 0.0,
- NAN,
- scale };
- }
- return { name,
- (1.0 - probability.array()).log() * scale,
- NAN,
- 0.0,
- scale };
-}
-
-double
-true_pdf(const std::string& name, double position, double scale)
-{
- if (name == "uniform")
- return 1.0 / scale;
- if (name == "beta-2-1")
- return 2.0 * position / scale;
- if (name == "beta-1-2")
- return 2.0 * (1.0 - position) / scale;
- return std::exp(-position) / scale;
-}
-
-double
-true_cdf(const std::string& name, double position)
-{
- if (name == "uniform")
- return position;
- if (name == "beta-2-1")
- return position * position;
- if (name == "beta-1-2")
- return 1.0 - (1.0 - position) * (1.0 - position);
- if (name == "exponential-left")
- return 1.0 - std::exp(-position);
- return std::exp(-position);
-}
-
-double
-evaluation_point(const Scenario& scenario, double position)
-{
- if (scenario.name == "exponential-right")
- return -position * scenario.scale;
- return position * scenario.scale;
-}
-
-void
-run_scenario(const Scenario& scenario,
- size_t degree,
- size_t grid_size,
- double bandwidth,
- const Eigen::VectorXd& positions)
-{
- kde1d::Kde1d fit(scenario.xmin,
- scenario.xmax,
- "continuous",
- 1.0,
- bandwidth,
- degree,
- grid_size);
- fit.fit(scenario.observations);
-
- for (Eigen::Index i = 0; i < positions.size(); ++i) {
- const double point = evaluation_point(scenario, positions(i));
- const double estimated_pdf =
- fit.pdf(Eigen::VectorXd::Constant(1, point))(0);
- const double estimated_cdf =
- fit.cdf(Eigen::VectorXd::Constant(1, point))(0);
- std::cout << scenario.name << ',' << scenario.scale << ',' << degree << ','
- << grid_size << ','
- << (std::isnan(bandwidth) ? "automatic"
- : std::to_string(bandwidth))
- << ',' << positions(i) << ',' << point << ',' << estimated_pdf
- << ',' << true_pdf(scenario.name, positions(i), scenario.scale)
- << ',' << estimated_cdf << ','
- << true_cdf(scenario.name, positions(i)) << '\n';
- }
-}
-
-} // namespace
-
-int
-main()
-{
- std::cout << std::setprecision(17);
- std::cout << "scenario,scale,degree,grid_size,bandwidth,position,x,"
- "estimated_pdf,true_pdf,estimated_cdf,true_cdf\n";
-
- Eigen::VectorXd positions(12);
- positions << 0.0, 1e-12, 1e-10, 1e-8, 1e-6, 1e-5, 1e-4, 1e-3,
- 1e-2, 1e-1, 0.5, 1.0;
- const std::array scenario_names = {
- "uniform", "beta-2-1", "beta-1-2", "exponential-left",
- "exponential-right"
- };
-
- for (const std::string& scenario_name : scenario_names) {
- for (double scale : { 1e-4, 1.0, 1e4 }) {
- Scenario scenario = make_scenario(scenario_name, 2000, scale);
- for (size_t degree = 0; degree < 3; ++degree) {
- for (size_t grid_size : { 100, 400, 1000 }) {
- for (double bandwidth :
- { 0.15, 0.3, std::numeric_limits::quiet_NaN() }) {
- run_scenario(
- scenario, degree, grid_size, bandwidth, positions);
- }
- }
- }
- }
- }
-}
diff --git a/diagnostic/extract_native_boundary_examples.R b/diagnostic/extract_native_boundary_examples.R
new file mode 100644
index 0000000..07ce3b1
--- /dev/null
+++ b/diagnostic/extract_native_boundary_examples.R
@@ -0,0 +1,116 @@
+arguments <- commandArgs(trailingOnly = TRUE)
+if (length(arguments) != 4L) {
+ stop("usage: script method output.csv library-path selections.csv")
+}
+.libPaths(c(arguments[[3L]], .libPaths()))
+library(kde1d)
+
+method <- arguments[[1L]]
+if (!method %in% c("pre-pr", "bulk", "expert")) {
+ stop("method must be one of: pre-pr, bulk, expert")
+}
+
+fit_kde1d <- function(observations, xmin, xmax) {
+ fit_arguments <- list(x = observations, xmin = xmin, xmax = xmax, deg = 2)
+ if (method != "pre-pr") {
+ fit_arguments$boundary_repair <- method == "expert"
+ }
+ do.call(kde1d, fit_arguments)
+}
+
+one_sided <- list(
+ exponential = list(random = rexp, density = dexp, quantile = qexp),
+ half_normal = list(
+ random = function(n) abs(rnorm(n)),
+ density = function(x) 2 * dnorm(x),
+ quantile = function(p) qnorm((p + 1) / 2)
+ ),
+ gamma_075 = list(
+ random = function(n) rgamma(n, 0.75),
+ density = function(x) dgamma(x, 0.75),
+ quantile = function(p) qgamma(p, 0.75)
+ ),
+ gamma_2 = list(
+ random = function(n) rgamma(n, 2),
+ density = function(x) dgamma(x, 2),
+ quantile = function(p) qgamma(p, 2)
+ )
+)
+two_sided <- list(
+ uniform = list(shape1 = 1, shape2 = 1),
+ beta_1_2 = list(shape1 = 1, shape2 = 2),
+ beta_075_2 = list(shape1 = 0.75, shape2 = 2),
+ beta_2_2 = list(shape1 = 2, shape2 = 2)
+)
+selections <- read.csv(arguments[[4L]], stringsAsFactors = FALSE)
+
+rows <- list()
+for (scenario_index in seq_along(one_sided)) {
+ scenario <- one_sided[[scenario_index]]
+ selection <- selections[
+ selections$support == "one-sided" &
+ selections$scenario == names(one_sided)[[scenario_index]],
+ ]
+ sample_size_index <- match(selection$sample_size, c(25L, 100L, 1000L, 2000L))
+ set.seed(20260824L + 10000000L * scenario_index +
+ 1000L * sample_size_index + selection$replication)
+ observations <- scenario$random(selection$sample_size)
+ points <- sort(unique(c(
+ 1e-10,
+ scenario$quantile(10^seq(-8, -1, length.out = 120L)),
+ scenario$quantile(seq(0.1, 0.999, length.out = 300L))
+ )))
+ if (is.finite(scenario$density(0))) {
+ points <- c(0, points)
+ }
+ fit <- fit_kde1d(observations, 0, NaN)
+ rows[[length(rows) + 1L]] <- data.frame(
+ method = method,
+ support = "one-sided",
+ scenario = paste0(
+ names(one_sided)[[scenario_index]], "_n", selection$sample_size
+ ),
+ sample_size = selection$sample_size,
+ x = points,
+ truth = scenario$density(points),
+ estimate = dkde1d(points, fit)
+ )
+}
+
+points <- sort(unique(c(
+ 10^seq(-10, -2, length.out = 121L),
+ seq(0.01, 0.99, length.out = 301L),
+ 1 - 10^seq(-2, -10, length.out = 121L)
+)))
+for (scenario_index in seq_along(two_sided)) {
+ scenario <- two_sided[[scenario_index]]
+ selection <- selections[
+ selections$support == "two-sided" &
+ selections$scenario == names(two_sided)[[scenario_index]],
+ ]
+ sample_size_index <- match(selection$sample_size, c(25L, 100L, 1000L, 2000L))
+ set.seed(20260825L + 100000L * scenario_index +
+ 1000L * sample_size_index + selection$replication)
+ observations <- rbeta(
+ selection$sample_size, scenario$shape1, scenario$shape2
+ )
+ scenario_points <- points
+ if (scenario$shape1 >= 1 && scenario$shape2 >= 1) {
+ scenario_points <- c(0, points, 1)
+ }
+ fit <- fit_kde1d(observations, 0, 1)
+ rows[[length(rows) + 1L]] <- data.frame(
+ method = method,
+ support = "two-sided",
+ scenario = paste0(
+ names(two_sided)[[scenario_index]], "_n", selection$sample_size
+ ),
+ sample_size = selection$sample_size,
+ x = scenario_points,
+ truth = dbeta(scenario_points, scenario$shape1, scenario$shape2),
+ estimate = dkde1d(scenario_points, fit)
+ )
+}
+
+dir.create(dirname(arguments[[2L]]), recursive = TRUE, showWarnings = FALSE)
+write.csv(do.call(rbind, rows), arguments[[2L]], row.names = FALSE)
diff --git a/diagnostic/plot_boundary_diagnostics.R b/diagnostic/plot_boundary_diagnostics.R
deleted file mode 100644
index f7b4517..0000000
--- a/diagnostic/plot_boundary_diagnostics.R
+++ /dev/null
@@ -1,101 +0,0 @@
-args <- commandArgs(trailingOnly = TRUE)
-if (length(args) != 1L) {
- stop("usage: Rscript diagnostic/plot_boundary_diagnostics.R diagnostics.csv")
-}
-
-results <- read.csv(args[[1]])
-results$pdf_ratio <- results$estimated_pdf / results$true_pdf
-results$pdf_ratio[!is.finite(results$pdf_ratio) | results$pdf_ratio <= 0] <-
- NA_real_
-results$cdf_error <- abs(results$estimated_cdf - results$true_cdf)
-
-tail_results <- subset(results, position <= 1e-2 & true_pdf > 0)
-summary <- aggregate(
- cbind(log_pdf_error = abs(log(pdf_ratio)), cdf_error = cdf_error) ~
- scenario + scale + degree + grid_size + bandwidth,
- data = tail_results,
- FUN = function(x) max(x, na.rm = TRUE)
-)
-summary <- summary[order(summary$log_pdf_error, summary$cdf_error), ]
-cat("Worst tail errors\n")
-print(head(summary[order(-summary$log_pdf_error), ], 20L), row.names = FALSE)
-
-results$scaled_pdf <- results$estimated_pdf * results$scale
-scale_summary <- aggregate(
- cbind(scaled_pdf = scaled_pdf, estimated_cdf = estimated_cdf) ~
- scenario + degree + grid_size + bandwidth + position,
- data = subset(results, true_pdf > 0),
- FUN = function(x) max(x) - min(x)
-)
-scale_summary$discrepancy <- pmax(
- scale_summary$scaled_pdf,
- scale_summary$estimated_cdf
-)
-cat("\nLargest scale-equivariance discrepancies\n")
-print(
- head(scale_summary[order(-scale_summary$discrepancy), ], 20L),
- row.names = FALSE
-)
-
-left <- subset(results, scenario == "exponential-left")
-right <- subset(results, scenario == "exponential-right")
-reflection <- merge(
- left,
- right,
- by = c("scale", "degree", "grid_size", "bandwidth", "position"),
- suffixes = c("_left", "_right")
-)
-reflection$pdf_error <- abs(
- reflection$scaled_pdf_left - reflection$scaled_pdf_right
-)
-reflection$cdf_error <- abs(
- reflection$estimated_cdf_left + reflection$estimated_cdf_right - 1
-)
-cat("\nReflection-equivariance maxima\n")
-cat("PDF:", max(reflection$pdf_error), "CDF:", max(reflection$cdf_error), "\n")
-
-grid_summary <- aggregate(
- cbind(scaled_pdf = scaled_pdf, estimated_cdf = estimated_cdf) ~
- scenario + scale + degree + bandwidth + position,
- data = subset(results, true_pdf > 0),
- FUN = function(x) max(x) - min(x)
-)
-grid_summary$discrepancy <- pmax(
- grid_summary$scaled_pdf,
- grid_summary$estimated_cdf
-)
-cat("\nLargest grid-size discrepancies\n")
-print(
- head(grid_summary[order(-grid_summary$discrepancy), ], 20L),
- row.names = FALSE
-)
-
-pdf(sub("[.]csv$", ".pdf", args[[1]]))
-for (scenario_name in unique(results$scenario)) {
- selected <- subset(
- results,
- scenario == scenario_name & scale == 1 & grid_size == 400 &
- bandwidth == "0.300000" & position > 0 & position <= 0.1
- )
- matplot(
- selected$position[selected$degree == 0],
- sapply(0:2, function(polynomial_degree) {
- selected$estimated_pdf[selected$degree == polynomial_degree]
- }),
- type = "l",
- log = "x",
- lty = 1,
- xlab = "distance from boundary / support scale",
- ylab = "density",
- main = scenario_name
- )
- lines(
- selected$position[selected$degree == 0],
- selected$true_pdf[selected$degree == 0],
- lty = 2,
- lwd = 2
- )
- legend("topright", c("degree 0", "degree 1", "degree 2", "truth"),
- col = c(1:3, 1), lty = c(1, 1, 1, 2), bty = "n")
-}
-dev.off()
diff --git a/diagnostic/plot_native_boundary_examples.R b/diagnostic/plot_native_boundary_examples.R
new file mode 100644
index 0000000..25b8182
--- /dev/null
+++ b/diagnostic/plot_native_boundary_examples.R
@@ -0,0 +1,25 @@
+arguments <- commandArgs(trailingOnly = TRUE)
+if (length(arguments) != 4L) {
+ stop("usage: script pre-pr.csv bulk.csv expert.csv output.png")
+}
+library(ggplot2)
+
+estimates <- do.call(rbind, lapply(arguments[1:3], read.csv))
+estimates$method <- factor(
+ estimates$method,
+ levels = c("pre-pr", "bulk", "expert"),
+ labels = c("pre-PR log/probit", "new bulk", "new bulk + experts")
+)
+truth <- unique(estimates[c("support", "scenario", "sample_size", "x", "truth")])
+
+plot <- ggplot(estimates, aes(x, estimate, color = method)) +
+ geom_line(linewidth = 0.55) +
+ geom_line(data = truth, aes(x = x, y = truth), inherit.aes = FALSE,
+ color = "black", linewidth = 0.65, linetype = 2) +
+ facet_wrap(~scenario, scales = "free", ncol = 2) +
+ scale_color_manual(values = c("#777777", "#377eb8", "#e41a1c")) +
+ labs(x = NULL, y = "density", color = NULL) +
+ theme_bw(base_size = 9) +
+ theme(legend.position = "bottom")
+
+ggsave(arguments[[4L]], plot, width = 8, height = 10, dpi = 160)
diff --git a/diagnostic/run_boundary_diagnostics.sh b/diagnostic/run_boundary_diagnostics.sh
deleted file mode 100755
index a0ad5bb..0000000
--- a/diagnostic/run_boundary_diagnostics.sh
+++ /dev/null
@@ -1,20 +0,0 @@
-#!/bin/sh
-set -eu
-
-diagnostic_source_dir=$(CDPATH= cd -- "$(dirname -- "$0")/.." && pwd)
-diagnostic_build_dir=${1:-"$diagnostic_source_dir/build-diagnostic"}
-
-cmake -S "$diagnostic_source_dir" -B "$diagnostic_build_dir" \
- -DCMAKE_BUILD_TYPE=Release \
- -DBUILD_TESTING=OFF \
- -DBUILD_DIAGNOSTICS=ON 1>&2
-cmake --build "$diagnostic_build_dir" \
- --target kde1d-boundary-diagnostic \
- --config Release 1>&2
-
-diagnostic_executable="$diagnostic_build_dir/bin/kde1d-boundary-diagnostic"
-if [ ! -x "$diagnostic_executable" ]; then
- diagnostic_executable="$diagnostic_build_dir/bin/Release/kde1d-boundary-diagnostic.exe"
-fi
-
-exec "$diagnostic_executable"
diff --git a/diagnostic/summarize_native_boundary_validation.R b/diagnostic/summarize_native_boundary_validation.R
new file mode 100644
index 0000000..ba24502
--- /dev/null
+++ b/diagnostic/summarize_native_boundary_validation.R
@@ -0,0 +1,228 @@
+arguments <- commandArgs(trailingOnly = TRUE)
+if (length(arguments) != 4L) {
+ stop("usage: script pre-pr.csv bulk.csv expert.csv output-prefix")
+}
+
+read_results <- function(path) {
+ results <- read.csv(path, stringsAsFactors = FALSE)
+ results$direction[is.na(results$direction)] <- "none"
+ results
+}
+
+pre_pr <- read_results(arguments[[1L]])
+bulk <- read_results(arguments[[2L]])
+expert <- read_results(arguments[[3L]])
+output_prefix <- arguments[[4L]]
+keys <- c(
+ "support", "direction", "scale", "scenario", "degree", "sample_size",
+ "replication"
+)
+
+paired_ratios <- function(reference, comparison) {
+ paired <- merge(
+ expert,
+ reference,
+ by = keys,
+ suffixes = c("_expert", "_reference")
+ )
+ groups <- split(
+ seq_len(nrow(paired)),
+ interaction(paired[c(
+ "support", "direction", "scale", "scenario", "degree", "sample_size"
+ )], drop = TRUE)
+ )
+ rows <- lapply(groups, function(index) {
+ lapply(c("global_ise", "boundary_ise"), function(metric) {
+ candidate <- paired[[paste0(metric, "_expert")]][index]
+ baseline <- paired[[paste0(metric, "_reference")]][index]
+ ratio <- mean(candidate) / mean(baseline)
+ data.frame(
+ comparison = comparison,
+ metric = metric,
+ paired[index[[1L]], c(
+ "support", "direction", "scale", "scenario", "degree",
+ "sample_size"
+ )],
+ ratio = ratio,
+ mcse = sd(candidate - ratio * baseline) /
+ (sqrt(length(index)) * mean(baseline))
+ )
+ })
+ })
+ do.call(rbind, unlist(rows, recursive = FALSE))
+}
+
+cell_ratios <- rbind(
+ paired_ratios(pre_pr, "pre-pr"),
+ paired_ratios(bulk, "bulk")
+)
+
+# Average scale/reflection cells within a family while retaining their paired
+# covariance. Different families use independent simulation streams.
+family_paired_ratios <- function(reference, comparison) {
+ paired <- merge(
+ expert,
+ reference,
+ by = keys,
+ suffixes = c("_expert", "_reference")
+ )
+ cell_names <- c(
+ "support", "direction", "scale", "scenario", "degree", "sample_size"
+ )
+ family_names <- c("support", "scenario", "degree", "sample_size")
+ rows <- lapply(c("global_ise", "boundary_ise"), function(metric) {
+ paired$cell <- interaction(paired[cell_names], drop = TRUE)
+ paired$influence <- NA_real_
+ paired$cell_ratio <- NA_real_
+ for (index in split(seq_len(nrow(paired)), paired$cell)) {
+ candidate <- paired[[paste0(metric, "_expert")]][index]
+ baseline <- paired[[paste0(metric, "_reference")]][index]
+ ratio <- mean(candidate) / mean(baseline)
+ paired$cell_ratio[index] <- ratio
+ paired$influence[index] <-
+ (candidate - ratio * baseline) / mean(baseline)
+ }
+ family_groups <- split(
+ seq_len(nrow(paired)), interaction(paired[family_names], drop = TRUE)
+ )
+ lapply(family_groups, function(index) {
+ first_cell_rows <- index[!duplicated(paired$cell[index])]
+ replication_groups <- split(index, paired$replication[index])
+ replication_influence <- vapply(
+ replication_groups,
+ function(replication_index) mean(paired$influence[replication_index]),
+ numeric(1)
+ )
+ data.frame(
+ comparison = comparison,
+ metric = metric,
+ paired[index[[1L]], family_names, drop = FALSE],
+ ratio = mean(paired$cell_ratio[first_cell_rows]),
+ mcse = sd(replication_influence) / sqrt(length(replication_influence)),
+ cells = length(first_cell_rows)
+ )
+ })
+ })
+ do.call(rbind, unlist(rows, recursive = FALSE))
+}
+
+family_ratios <- rbind(
+ family_paired_ratios(pre_pr, "pre-pr"),
+ family_paired_ratios(bulk, "bulk")
+)
+support_names <- c(
+ "comparison", "metric", "support", "degree", "sample_size"
+)
+support_groups <- split(
+ seq_len(nrow(family_ratios)),
+ interaction(family_ratios[support_names], drop = TRUE)
+)
+support_ratios <- do.call(rbind, lapply(support_groups, function(index) {
+ data.frame(
+ family_ratios[index[[1L]], support_names, drop = FALSE],
+ ratio = mean(family_ratios$ratio[index]),
+ mcse = sqrt(sum(family_ratios$mcse[index]^2)) / length(index),
+ cells = sum(family_ratios$cells[index])
+ )
+}))
+
+classification <- expert[
+ expert$degree == 0 &
+ ((expert$support == "one-sided" & expert$scale == 1 &
+ expert$direction == "lower") | expert$support == "two-sided"),
+]
+activation_groups <- split(
+ seq_len(nrow(classification)),
+ interaction(classification[c("support", "scenario", "sample_size")],
+ drop = TRUE)
+)
+activation <- do.call(rbind, lapply(activation_groups, function(index) {
+ data.frame(
+ classification[index[[1L]], c("support", "scenario", "sample_size")],
+ lower_rate = mean(classification$repair_lower[index]),
+ upper_rate = mean(classification$repair_upper[index]),
+ any_rate = mean(
+ classification$repair_lower[index] | classification$repair_upper[index]
+ ),
+ both_rate = mean(
+ classification$repair_lower[index] & classification$repair_upper[index]
+ )
+ )
+}))
+
+read_summaries <- function(path) {
+ summary <- read.csv(path, stringsAsFactors = FALSE)
+ summary$direction[is.na(summary$direction)] <- "none"
+ summary
+}
+summary_path <- function(path) sub("\\.csv$", "-summary.csv", path)
+pre_pr_summary <- read_summaries(summary_path(arguments[[1L]]))
+bulk_summary <- read_summaries(summary_path(arguments[[2L]]))
+expert_summary <- read_summaries(summary_path(arguments[[3L]]))
+summary_keys <- keys[keys != "replication"]
+
+component_ratios <- function(reference, comparison) {
+ paired <- merge(
+ expert_summary,
+ reference,
+ by = summary_keys,
+ suffixes = c("_expert", "_reference")
+ )
+ metrics <- c(
+ "global_squared_bias", "global_variance",
+ "boundary_squared_bias", "boundary_variance"
+ )
+ do.call(rbind, lapply(metrics, function(metric) {
+ data.frame(
+ comparison = comparison,
+ metric = metric,
+ paired[summary_keys],
+ ratio = paired[[paste0(metric, "_expert")]] /
+ paired[[paste0(metric, "_reference")]]
+ )
+ }))
+}
+
+components <- rbind(
+ component_ratios(pre_pr_summary, "pre-pr"),
+ component_ratios(bulk_summary, "bulk")
+)
+
+example_targets <- data.frame(
+ support = c(rep("one-sided", 4), rep("two-sided", 4)),
+ scenario = c(
+ "exponential", "half_normal", "gamma_075", "gamma_2",
+ "uniform", "beta_1_2", "beta_075_2", "beta_2_2"
+ ),
+ sample_size = c(100L, 25L, 1000L, 1000L, 100L, 1000L, 1000L, 1000L),
+ repair_lower = c(TRUE, TRUE, FALSE, FALSE, TRUE, TRUE, FALSE, FALSE),
+ repair_upper = c(FALSE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE)
+)
+example_rows <- lapply(seq_len(nrow(example_targets)), function(i) {
+ target <- example_targets[i, ]
+ candidates <- expert[
+ expert$support == target$support &
+ expert$scenario == target$scenario &
+ expert$sample_size == target$sample_size &
+ expert$degree == 2 &
+ expert$scale == 1 &
+ expert$repair_lower == target$repair_lower &
+ expert$repair_upper == target$repair_upper &
+ expert$direction ==
+ if (target$support == "one-sided") "lower" else "none",
+ ]
+ median_ise <- median(candidates$global_ise)
+ selected <- candidates[
+ which.min(abs(candidates$global_ise - median_ise)),
+ ]
+ data.frame(target, replication = selected$replication)
+})
+examples <- do.call(rbind, example_rows)
+
+dir.create(dirname(output_prefix), recursive = TRUE, showWarnings = FALSE)
+write.csv(cell_ratios, paste0(output_prefix, "-cells.csv"), row.names = FALSE)
+write.csv(family_ratios, paste0(output_prefix, "-families.csv"), row.names = FALSE)
+write.csv(support_ratios, paste0(output_prefix, "-support.csv"), row.names = FALSE)
+write.csv(activation, paste0(output_prefix, "-activation.csv"), row.names = FALSE)
+write.csv(components, paste0(output_prefix, "-components.csv"), row.names = FALSE)
+write.csv(examples, paste0(output_prefix, "-examples.csv"), row.names = FALSE)
diff --git a/diagnostic/validate_native_boundary_experts.R b/diagnostic/validate_native_boundary_experts.R
new file mode 100644
index 0000000..22c2e1d
--- /dev/null
+++ b/diagnostic/validate_native_boundary_experts.R
@@ -0,0 +1,263 @@
+# Evaluate one native kde1d method. Install the pre-PR revision for `pre-pr`;
+# install the candidate revision for `bulk` and `expert`. The pre-PR one-sided
+# fit is the production log-transform estimator, not a local emulation of it.
+
+arguments <- commandArgs(trailingOnly = TRUE)
+if (length(arguments) < 2L || length(arguments) > 4L) {
+ stop("usage: script method output.csv [replications] [library-path]")
+}
+if (length(arguments) == 4L) {
+ .libPaths(c(arguments[[4L]], .libPaths()))
+}
+library(kde1d)
+
+method <- arguments[[1L]]
+output_file <- arguments[[2L]]
+if (!method %in% c("pre-pr", "bulk", "expert")) {
+ stop("method must be one of: pre-pr, bulk, expert")
+}
+replications <- if (length(arguments) >= 3L) {
+ as.integer(arguments[[3L]])
+} else {
+ 100L
+}
+if (!is.finite(replications) || replications < 2L) {
+ stop("replications must be at least two")
+}
+
+integrate_values <- function(points, values) {
+ sum(diff(points) * (values[-length(values)] + values[-1L]) / 2)
+}
+
+fit_kde1d <- function(observations, degree,
+ xmin = NA_real_, xmax = NA_real_) {
+ arguments <- list(x = observations, xmin = xmin, xmax = xmax, deg = degree)
+ if (method != "pre-pr") {
+ arguments$boundary_repair <- method == "expert"
+ }
+ do.call(kde1d, arguments)
+}
+
+classify_finite_endpoint <- function(distances) {
+ distances <- sort(distances)
+ k <- min(length(distances) - 1L, ceiling(2 * sqrt(length(distances))))
+ reference_distance <- distances[[k + 1L]]
+ if (!is.finite(reference_distance) || reference_distance <= 0) {
+ return(FALSE)
+ }
+ denominator <- sum(log(
+ reference_distance /
+ pmax(distances[seq_len(k)], .Machine$double.eps * reference_distance)
+ ))
+ if (!is.finite(denominator) || denominator <= 0) {
+ return(FALSE)
+ }
+ beta <- k / denominator
+ beta >= 0.9 && beta * (1 - qnorm(0.95) / sqrt(k)) <= 1
+}
+
+one_sided <- list(
+ exponential = list(
+ random = function(n) rexp(n),
+ density = dexp,
+ quantile = qexp
+ ),
+ half_normal = list(
+ random = function(n) abs(rnorm(n)),
+ density = function(x) 2 * dnorm(x),
+ quantile = function(p) qnorm((p + 1) / 2)
+ ),
+ gamma_2 = list(
+ random = function(n) rgamma(n, 2),
+ density = function(x) dgamma(x, 2),
+ quantile = function(p) qgamma(p, 2)
+ ),
+ gamma_075 = list(
+ random = function(n) rgamma(n, 0.75),
+ density = function(x) dgamma(x, 0.75),
+ quantile = function(p) qgamma(p, 0.75)
+ ),
+ lomax_2 = list(
+ random = function(n) (1 - runif(n))^(-0.5) - 1,
+ density = function(x) 2 / (1 + x)^3,
+ quantile = function(p) (1 - p)^(-0.5) - 1
+ )
+)
+two_sided <- list(
+ uniform = c(1, 1),
+ beta_1_2 = c(1, 2),
+ beta_2_1 = c(2, 1),
+ beta_075_2 = c(0.75, 2),
+ beta_2_075 = c(2, 0.75),
+ beta_075_075 = c(0.75, 0.75),
+ beta_2_2 = c(2, 2)
+)
+sample_sizes <- c(25L, 100L, 1000L, 2000L)
+degrees <- 0:2
+one_sided_scales <- c(1e-4, 1, 1e4)
+one_sided_directions <- c("lower", "upper")
+result_rows <- list()
+summary_rows <- list()
+
+summarize_estimates <- function(points, estimates, truth, boundary) {
+ mean_density <- rowMeans(estimates)
+ c(
+ global_squared_bias = integrate_values(
+ points, (mean_density - truth)^2
+ ),
+ global_variance = integrate_values(
+ points, rowMeans((estimates - mean_density)^2)
+ ),
+ boundary_squared_bias = integrate_values(
+ points[boundary], (mean_density[boundary] - truth[boundary])^2
+ ),
+ boundary_variance = integrate_values(
+ points[boundary],
+ rowMeans((estimates[boundary, , drop = FALSE] -
+ mean_density[boundary])^2)
+ )
+ )
+}
+
+for (scenario_index in seq_along(one_sided)) {
+ scenario <- one_sided[[scenario_index]]
+ base_points <- sort(unique(c(
+ 1e-10,
+ scenario$quantile(10^seq(-8, -1, length.out = 100L)),
+ scenario$quantile(seq(0.1, 0.9999, length.out = 401L))
+ )))
+ for (scale_index in seq_along(one_sided_scales)) {
+ scale <- one_sided_scales[[scale_index]]
+ scaled_points <- scale * base_points
+ scaled_truth <- scenario$density(base_points) / scale
+ for (direction_index in seq_along(one_sided_directions)) {
+ direction <- one_sided_directions[[direction_index]]
+ if (direction == "lower") {
+ points <- scaled_points
+ truth <- scaled_truth
+ boundary <- points <= scale * scenario$quantile(0.1)
+ } else {
+ points <- -rev(scaled_points)
+ truth <- rev(scaled_truth)
+ boundary <- points >= -scale * scenario$quantile(0.1)
+ }
+ for (degree in degrees) {
+ for (sample_size in sample_sizes) {
+ estimates <- matrix(NA_real_, length(points), replications)
+ for (replication in seq_len(replications)) {
+ # Reuse the base sample across scales and reflections so departures
+ # from scale/reflection equivariance are paired directly.
+ set.seed(20260824L + 10000000L * scenario_index +
+ 1000L * match(sample_size, sample_sizes) + replication)
+ observations <- scale * scenario$random(sample_size)
+ repair_endpoint <- classify_finite_endpoint(observations)
+ if (direction == "lower") {
+ fit <- fit_kde1d(observations, degree, xmin = 0)
+ } else {
+ fit <- fit_kde1d(-observations, degree, xmax = 0)
+ }
+ estimates[, replication] <- dkde1d(points, fit)
+ result_rows[[length(result_rows) + 1L]] <- data.frame(
+ support = "one-sided",
+ direction = direction,
+ scale = scale,
+ scenario = names(one_sided)[[scenario_index]],
+ degree = degree,
+ sample_size = sample_size,
+ replication = replication,
+ method = method,
+ repair_lower = direction == "lower" && repair_endpoint,
+ repair_upper = direction == "upper" && repair_endpoint,
+ global_ise = integrate_values(
+ points, (estimates[, replication] - truth)^2
+ ),
+ boundary_ise = integrate_values(
+ points[boundary],
+ (estimates[boundary, replication] - truth[boundary])^2
+ ),
+ edf = fit$edf,
+ loglik = fit$loglik
+ )
+ }
+ summary_rows[[length(summary_rows) + 1L]] <- data.frame(
+ support = "one-sided",
+ direction = direction,
+ scale = scale,
+ scenario = names(one_sided)[[scenario_index]],
+ degree = degree,
+ sample_size = sample_size,
+ method = method,
+ as.list(summarize_estimates(points, estimates, truth, boundary))
+ )
+ }
+ }
+ }
+ }
+}
+
+points <- sort(unique(c(
+ 10^seq(-10, -2, length.out = 101L),
+ seq(0.01, 0.99, length.out = 301L),
+ 1 - 10^seq(-2, -10, length.out = 101L)
+)))
+boundary <- points <= 0.1 | points >= 0.9
+for (scenario_index in seq_along(two_sided)) {
+ shapes <- two_sided[[scenario_index]]
+ truth <- dbeta(points, shapes[[1L]], shapes[[2L]])
+ for (degree in degrees) {
+ for (sample_size in sample_sizes) {
+ estimates <- matrix(NA_real_, length(points), replications)
+ for (replication in seq_len(replications)) {
+ set.seed(20260825L + 100000L * scenario_index +
+ 1000L * match(sample_size, sample_sizes) + replication)
+ observations <- rbeta(sample_size, shapes[[1L]], shapes[[2L]])
+ fit <- fit_kde1d(
+ observations,
+ degree,
+ xmin = 0,
+ xmax = 1
+ )
+ estimates[, replication] <- dkde1d(points, fit)
+ result_rows[[length(result_rows) + 1L]] <- data.frame(
+ support = "two-sided",
+ direction = NA_character_,
+ scale = 1,
+ scenario = names(two_sided)[[scenario_index]],
+ degree = degree,
+ sample_size = sample_size,
+ replication = replication,
+ method = method,
+ repair_lower = classify_finite_endpoint(observations),
+ repair_upper = classify_finite_endpoint(1 - observations),
+ global_ise = integrate_values(
+ points, (estimates[, replication] - truth)^2
+ ),
+ boundary_ise = integrate_values(
+ points[boundary],
+ (estimates[boundary, replication] - truth[boundary])^2
+ ),
+ edf = fit$edf,
+ loglik = fit$loglik
+ )
+ }
+ summary_rows[[length(summary_rows) + 1L]] <- data.frame(
+ support = "two-sided",
+ direction = NA_character_,
+ scale = 1,
+ scenario = names(two_sided)[[scenario_index]],
+ degree = degree,
+ sample_size = sample_size,
+ method = method,
+ as.list(summarize_estimates(points, estimates, truth, boundary))
+ )
+ }
+ }
+}
+
+dir.create(dirname(output_file), recursive = TRUE, showWarnings = FALSE)
+write.csv(do.call(rbind, result_rows), output_file, row.names = FALSE)
+write.csv(
+ do.call(rbind, summary_rows),
+ sub("\\.csv$", "-summary.csv", output_file),
+ row.names = FALSE
+)
diff --git a/docs/overview-continuous.dox b/docs/overview-continuous.dox
index 28c1ed3..4d3dc0e 100644
--- a/docs/overview-continuous.dox
+++ b/docs/overview-continuous.dox
@@ -1,4 +1,114 @@
-/** @page overview-continuous Continuous Distributions
+/** @page overview-continuous Continuous distributions
+Continuous fits combine a transformed local-likelihood bulk estimator with an
+optional finite-endpoint repair. The public `degree`, `bandwidth`, and
+`multiplier` arguments retain their usual meaning for the bulk fit. Boundary
+repair is enabled by default and can be disabled with `boundary_repair=false`.
+
+@section continuous-bulk Bulk estimator and support transforms
+
+With no finite support endpoint, data are fitted on their original scale. If
+both endpoints @f$acheck_xmin_xmax(xmin, xmax);
if (multiplier <= 0.0) {
@@ -276,20 +357,24 @@ inline Kde1d::Kde1d(const interp::InterpolationGrid& grid,
//! selection).
//! @param degree degree of the local polynomial.
//! @param grid_size number of grid points for the interpolation grid.
+//! @param boundary_repair whether finite endpoints are eligible for boundary
+//! experts. A data-driven endpoint classifier can still retain the bulk fit.
inline Kde1d::Kde1d(double xmin,
double xmax,
std::string type,
double multiplier,
double bandwidth,
size_t degree,
- size_t grid_size)
+ size_t grid_size,
+ bool boundary_repair)
: Kde1d(xmin,
xmax,
this->as_enum(type),
multiplier,
bandwidth,
degree,
- grid_size)
+ grid_size,
+ boundary_repair)
{
}
@@ -327,6 +412,14 @@ Kde1d::fit(const Eigen::VectorXd& x, const Eigen::VectorXd& weights)
Eigen::VectorXd w = weights;
tools::remove_nans(xx, w);
+ const bool use_boundary_repair = boundary_repair_ &&
+ type_ == VarType::continuous &&
+ xx.size() >= 16 &&
+ (!std::isnan(xmin_) || !std::isnan(xmax_));
+ Eigen::VectorXd boundary_observations;
+ if (use_boundary_repair)
+ boundary_observations = xx;
+
if (w.size() > 0) {
w /= w.mean();
}
@@ -363,16 +456,16 @@ Kde1d::fit(const Eigen::VectorXd& x, const Eigen::VectorXd& weights)
else
distances = xx.array() - xmin_;
if (w.size() == 0) {
- boundary_offset_ = 1e-5 * stats::median(distances);
+ boundary_scale_ = stats::median(distances);
} else {
- boundary_offset_ =
- 1e-5 * stats::quantile(
- distances, Eigen::VectorXd::Constant(1, 0.5), w)(0);
+ boundary_scale_ = stats::quantile(
+ distances, Eigen::VectorXd::Constant(1, 0.5), w)(0);
}
- if (!(boundary_offset_ > 0.0))
- boundary_offset_ = 1e-5 * distances.maxCoeff();
- if (!(boundary_offset_ > 0.0))
- boundary_offset_ = 1e-5;
+ if (!(boundary_scale_ > 0.0))
+ boundary_scale_ = distances.maxCoeff();
+ if (!(boundary_scale_ > 0.0))
+ boundary_scale_ = 1.0;
+ boundary_offset_ = 1e-5 * boundary_scale_;
}
xx = boundary_transform(xx);
@@ -390,9 +483,15 @@ Kde1d::fit(const Eigen::VectorXd& x, const Eigen::VectorXd& weights)
// order grid points from left to right
grid_points = finalize_grid(grid_points);
+ Eigen::VectorXd influences = fitted.col(1);
+ if (std::isnan(xmin_) && !std::isnan(xmax_))
+ influences.reverseInPlace();
+
// construct interpolation grid
// (3 iterations for normalization to a proper density)
grid_ = interp::InterpolationGrid(grid_points, values, 3);
+ if (use_boundary_repair)
+ repair_boundaries(boundary_observations, grid_points, influences, w);
// calculate log-likelihood of final estimate
xx = boundary_transform(xx, true);
@@ -413,9 +512,7 @@ Kde1d::fit(const Eigen::VectorXd& x, const Eigen::VectorXd& weights)
}
// calculate effective degrees of freedom
- Eigen::VectorXd influences = fitted.col(1).cwiseMin(3.0).cwiseMax(0);
- if (std::isnan(xmin_) && !std::isnan(xmax_))
- influences.reverseInPlace();
+ influences = influences.cwiseMin(3.0).cwiseMax(0);
interp::InterpolationGrid infl_grid(grid_points, influences, 0);
influences = infl_grid.interpolate(xx).array();
edf_ = influences.sum() + static_cast(prob0_ > 0);
@@ -779,11 +876,17 @@ Kde1d::boundary_transform(const Eigen::VectorXd& x, bool inverse)
x_new = (x.array() - xmin_ + 5e-5 * rng) / (1.0001 * rng);
x_new = stats::qnorm(x_new);
} else if (!std::isnan(xmin_)) {
- // left boundary -> log transform
- x_new = (boundary_offset_ + x.array() - xmin_).log();
+ // left boundary -> regularized fourth-root transform
+ x_new = 4.0 * (((boundary_offset_ + x.array() - xmin_) /
+ boundary_scale_)
+ .pow(0.25) -
+ std::pow(boundary_offset_ / boundary_scale_, 0.25));
} else if (!std::isnan(xmax_)) {
- // right boundary -> negative log transform
- x_new = (boundary_offset_ + xmax_ - x.array()).log();
+ // right boundary -> reflected regularized fourth-root transform
+ x_new = 4.0 * (((boundary_offset_ + xmax_ - x.array()) /
+ boundary_scale_)
+ .pow(0.25) -
+ std::pow(boundary_offset_ / boundary_scale_, 0.25));
} else {
// no boundary -> no transform
}
@@ -793,11 +896,19 @@ Kde1d::boundary_transform(const Eigen::VectorXd& x, bool inverse)
auto rng = xmax_ - xmin_;
x_new = stats::pnorm(x).array() * 1.0001 * rng + xmin_ - 5e-5 * rng;
} else if (!std::isnan(xmin_)) {
- // left boundary -> log transform
- x_new = x.array().exp() + xmin_ - boundary_offset_;
+ // left boundary -> inverse regularized fourth-root transform
+ x_new = boundary_scale_ *
+ (x.array() / 4.0 +
+ std::pow(boundary_offset_ / boundary_scale_, 0.25))
+ .pow(4) +
+ xmin_ - boundary_offset_;
} else if (!std::isnan(xmax_)) {
- // right boundary -> negative log transform
- x_new = -(x.array().exp() - xmax_ - boundary_offset_);
+ // right boundary -> inverse reflected regularized fourth-root transform
+ x_new = xmax_ + boundary_offset_ -
+ boundary_scale_ *
+ (x.array() / 4.0 +
+ std::pow(boundary_offset_ / boundary_scale_, 0.25))
+ .pow(4);
} else {
// no boundary -> no transform
}
@@ -824,14 +935,18 @@ Kde1d::boundary_correct(const Eigen::VectorXd& x, const Eigen::VectorXd& fhat)
auto rng = xmax_ - xmin_;
corr_term = (x.array() - xmin_ + 5e-5 * rng) / (xmax_ - xmin_ + 1e-4 * rng);
corr_term = stats::dnorm(stats::qnorm(corr_term));
- corr_term /= (xmax_ - xmin_ + 1e-4 * rng);
+ corr_term *= (xmax_ - xmin_ + 1e-4 * rng);
corr_term = 1.0 / corr_term.array();
} else if (!std::isnan(xmin_)) {
- // left boundary -> log transform
- corr_term = 1.0 / (boundary_offset_ + x.array() - xmin_);
+ // left boundary -> fourth-root-transform Jacobian
+ corr_term = ((boundary_offset_ + x.array() - xmin_) / boundary_scale_)
+ .pow(-0.75) /
+ boundary_scale_;
} else if (!std::isnan(xmax_)) {
- // right boundary -> negative log transform
- corr_term = 1.0 / (boundary_offset_ + xmax_ - x.array());
+ // right boundary -> reflected fourth-root-transform Jacobian
+ corr_term = ((boundary_offset_ + xmax_ - x.array()) / boundary_scale_)
+ .pow(-0.75) /
+ boundary_scale_;
} else {
// no boundary -> no transform
corr_term.fill(1.0);
@@ -897,11 +1012,403 @@ Kde1d::finalize_grid(Eigen::VectorXd& grid_points)
return grid_points;
}
+//! Classifies a support endpoint from its ordered distances @f$d@f$ to the
+//! endpoint. For @f$k = \min(n-1, \lceil 2\sqrt n \rceil)@f$, the lower-tail
+//! index is
+//! @f[
+//! \widehat\beta = \frac{k}{\sum_{i=1}^k
+//! \log\{d_{(k+1)} / d_{(i)}\}}.
+//! @f]
+//! It is treated as finite when @f$\widehat\beta \geq 0.9@f$ and its one-sided
+//! 95% lower confidence bound does not exceed one. The weighted threshold is
+//! 0.975 to avoid unstable repair of exploding targets. Ambiguous and
+//! numerically degenerate cases retain the bulk fit.
+//! With nonconstant case weights, their effective sample size replaces @f$n@f$
+//! and the lower tail is defined by cumulative normalized weight.
+//! @param dist ordered distances of the observations from the support endpoint.
+//! @param weights correspondingly ordered case weights with mean one.
+//! @return whether the endpoint is confidently classified as finite.
+inline bool
+Kde1d::is_finite_endpoint(const Eigen::VectorXd& dist,
+ const Eigen::VectorXd& weights) const
+{
+ const double n_eff = weights.sum() * weights.sum() / weights.squaredNorm();
+ if (!(n_eff > 1.0) || !std::isfinite(n_eff))
+ return false;
+
+ // Scale weights to effective counts and locate the (k+1)th weighted order
+ // statistic, where k = min(n_eff, ceil(2 sqrt(n_eff))).
+ const Eigen::VectorXd counts = weights * n_eff / weights.sum();
+ const double target =
+ std::min(std::nextafter(n_eff, 0.0), std::ceil(2.0 * std::sqrt(n_eff)));
+ Eigen::Index k1 = 0;
+ double cumulative = counts(0);
+ while (cumulative <= target && k1 + 1 < dist.size())
+ cumulative += counts(++k1);
+ const double dist_k1 = dist(k1);
+ if (!(dist_k1 > 0.0) || !std::isfinite(dist_k1))
+ return false;
+
+ // tail_mass and denom define the weighted lower-tail index beta.
+ double tail_mass = 0.0;
+ double denom = 0.0;
+ const double dist_min = std::numeric_limits::epsilon() * dist_k1;
+ for (Eigen::Index i = 0; i < k1; ++i) {
+ tail_mass += counts(i);
+ denom += counts(i) * std::log(dist_k1 / std::max(dist(i), dist_min));
+ }
+ if (!(tail_mass > 0.0) || !(denom > 0.0) || !std::isfinite(denom))
+ return false;
+
+ const double beta = tail_mass / denom;
+ // 1.64485 is the 95% standard-normal quantile used by the R selector.
+ const double beta_min =
+ weights.minCoeff() == weights.maxCoeff() ? 0.9 : 0.975;
+ return beta >= beta_min &&
+ beta * (1.0 - 1.6448536269514722 / std::sqrt(tail_mass)) <= 1.0;
+}
+
+//! Constructs sorted distances and case weights for one support endpoint, then
+//! classifies its local density behavior.
+//! @param x observations on the original scale.
+//! @param weights normalized case weights with mean one.
+//! @param side endpoint from which distances are measured.
+//! @return sorted endpoint data and its finite-endpoint decision.
+inline Kde1d::EndpointData
+Kde1d::prepare_endpoint(const Eigen::VectorXd& x,
+ const Eigen::VectorXd& weights,
+ EndpointSide side) const
+{
+ EndpointData endpoint;
+ if (side == EndpointSide::lower)
+ endpoint.dist = x.array() - xmin_;
+ else
+ endpoint.dist = xmax_ - x.array();
+ endpoint.weights.resize(x.size());
+ const Eigen::VectorXi order = tools::get_order(endpoint.dist);
+ const Eigen::VectorXd unsorted_dist = endpoint.dist;
+ for (Eigen::Index i = 0; i < endpoint.dist.size(); ++i) {
+ endpoint.dist(i) = unsorted_dist(order(i));
+ endpoint.weights(i) = weights(order(i));
+ }
+ endpoint.finite = is_finite_endpoint(endpoint.dist, endpoint.weights);
+ return endpoint;
+}
+
+//! Selects the shared original-scale bandwidth for a boundary expert. The
+//! closest `fraction` of positive-weight observations are used by row count.
+//! @param endpoint sorted endpoint distances and case weights.
+//! @param fraction fraction of positive-weight observations used for selection.
+//! @return selected boundary bandwidth after applying the public multiplier.
+inline double
+Kde1d::select_boundary_bandwidth(const EndpointData& endpoint,
+ double fraction) const
+{
+ const Eigen::Index n_positive = (endpoint.weights.array() > 0.0).count();
+ Eigen::VectorXd dist(n_positive);
+ Eigen::VectorXd weights(n_positive);
+ Eigen::Index j = 0;
+ for (Eigen::Index i = 0; i < endpoint.dist.size(); ++i) {
+ if (endpoint.weights(i) > 0.0) {
+ dist(j) = endpoint.dist(i);
+ weights(j++) = endpoint.weights(i);
+ }
+ }
+ const Eigen::Index n_bw = std::min(
+ dist.size(),
+ std::max(
+ 4,
+ static_cast(
+ std::ceil(fraction * static_cast(dist.size())))));
+ bandwidth::PluginBandwidthSelector selector(dist.head(n_bw),
+ weights.head(n_bw));
+ return selector.select_bandwidth(2) * multiplier_;
+}
+
+//! Computes the ordinary Gaussian KDE and its first derivative on the
+//! regular FFT grid used by a boundary component. The influence weights are
+//! the average normalized case weights in the same bins.
+//! @param dist ordered endpoint distances.
+//! @param weights correspondingly ordered normalized case weights.
+//! @param bandwidth boundary-component bandwidth.
+//! @param eval_upper largest endpoint distance that will be evaluated.
+//! @return FFT convolutions, influence weights, and FFT-grid upper limit.
+inline Kde1d::BoundaryConvolutions
+Kde1d::compute_boundary_convolutions(const Eigen::VectorXd& dist,
+ const Eigen::VectorXd& weights,
+ double bandwidth,
+ double eval_upper) const
+{
+ BoundaryConvolutions convolutions;
+ convolutions.upper = eval_upper + 6.0 * bandwidth;
+ const double* last = std::upper_bound(
+ dist.data(), dist.data() + dist.size(), convolutions.upper);
+ const Eigen::Index n_fft = last - dist.data();
+ const double fft_weight = weights.head(n_fft).sum();
+ if (n_fft == 0 || !(fft_weight > 0.0))
+ return convolutions;
+
+ // This resolution kept the worst paired ISE perturbation below 0.7%.
+ const size_t num_bins = 256;
+ fft::KdeFFT kde_fft(dist.head(n_fft),
+ bandwidth,
+ 0.0,
+ convolutions.upper,
+ weights.head(n_fft),
+ num_bins);
+ const double scale = fft_weight / weights.sum();
+ convolutions.f0 = scale * kde_fft.kde_drv(0);
+ convolutions.f1 = scale * kde_fft.kde_drv(1);
+ convolutions.influence_weight = Eigen::VectorXd::Ones(num_bins + 1);
+ if (weights.minCoeff() != weights.maxCoeff()) {
+ const Eigen::VectorXd count = tools::linbin(dist.head(n_fft),
+ 0.0,
+ convolutions.upper,
+ num_bins,
+ Eigen::VectorXd::Ones(n_fft));
+ convolutions.influence_weight =
+ (count.array() == 0.0)
+ .select(Eigen::VectorXd::Zero(count.size()),
+ weights.head(n_fft).mean() *
+ kde_fft.get_bin_counts().cwiseQuotient(count));
+ }
+ return convolutions;
+}
+
+//! Computes the first inverse-moment row for the local-linear Gaussian
+//! boundary kernel at standardized endpoint distance @f$a@f$.
+//! @param a endpoint distance divided by bandwidth.
+//! @return coefficients @f$c_0,c_1@f$.
+inline Kde1d::BoundaryKernelCoefficients
+Kde1d::boundary_kernel_coefficients(double a) const
+{
+ const double phi = K0_ * std::exp(-0.5 * a * a);
+ const double mu0 = 0.5 * std::erfc(-a * 0.7071067811865475);
+ const double mu1 = -phi;
+ const double mu2 = mu0 - a * phi;
+ const double denom = mu0 * mu2 - mu1 * mu1;
+ return { mu2 / denom, -mu1 / denom };
+}
+
+//! Fits a Gaussian local-linear equivalent kernel from sample distances
+//! @f$d@f$ to endpoint-distance grid @f$t@f$. For bandwidth @f$h@f$,
+//! truncated Gaussian moments @f$\mu_r(a)@f$, and
+//! @f$M(a)=[\mu_{j+k}(a)]_{j,k=0}^1@f$,
+//! @f[
+//! \begin{align*}
+//! a &= t/h, & u_i &= (t-d_i)/h, \\[2pt]
+//! \widehat f_E(t;h)
+//! &= \frac{1}{nh}\sum_{i=1}^n w_i\phi(u_i)\,
+//! e_0^\mathsf{T}M(a)^{-1}(1,u_i)^\mathsf{T}.
+//! \end{align*}
+//! @f]
+//! Writing @f$f^{(1)}@f$ for the first derivative of the ordinary Gaussian KDE
+//! and @f$(c_0,c_1)=e_0^\mathsf{T}M(a)^{-1}@f$, the same estimator is
+//! @f[
+//! \begin{align*}
+//! \widehat f_E(t;h) &= c_0f(t)-hc_1f^{(1)}(t).
+//! \end{align*}
+//! @f]
+//! This identity evaluates the two convolutions by FFT. The returned influence
+//! numerator is the diagonal contribution with normalized case weights
+//! @f$w_i@f$.
+//! @param dist ordered distances of the observations from the support endpoint.
+//! @param eval_dist endpoint distances at which to evaluate the component.
+//! @param bandwidth bandwidth of the equivalent kernel.
+//! @param weights correspondingly ordered case weights with mean one.
+//! @return the density and influence numerator on `eval_dist`.
+inline Kde1d::BoundaryComponent
+Kde1d::fit_boundary_component(const Eigen::VectorXd& dist,
+ const Eigen::VectorXd& eval_dist,
+ double bandwidth,
+ const Eigen::VectorXd& weights) const
+{
+ const double h = bandwidth;
+ if (!(h > 0.0) || !std::isfinite(h))
+ return {};
+
+ BoundaryComponent fit{ Eigen::VectorXd::Zero(eval_dist.size()),
+ Eigen::VectorXd::Zero(eval_dist.size()) };
+ if (eval_dist.size() == 0)
+ return fit;
+
+ const BoundaryConvolutions convolutions = compute_boundary_convolutions(
+ dist, weights, h, eval_dist(eval_dist.size() - 1));
+ if (convolutions.f0.size() == 0)
+ return fit;
+
+ const size_t num_bins = convolutions.f0.size() - 1;
+ for (Eigen::Index j = 0; j < eval_dist.size(); ++j) {
+ const double a = eval_dist(j) / h;
+ const BoundaryKernelCoefficients c = boundary_kernel_coefficients(a);
+ const double position =
+ eval_dist(j) * static_cast(num_bins) / convolutions.upper;
+ const size_t bin = std::min(static_cast(position), num_bins - 1);
+ const double fraction = position - static_cast(bin);
+ auto interpolate = [&](const Eigen::VectorXd& values) {
+ return (1.0 - fraction) * values(bin) + fraction * values(bin + 1);
+ };
+ fit.density(j) = c.c0 * interpolate(convolutions.f0) -
+ h * c.c1 * interpolate(convolutions.f1);
+ fit.influence_num(j) =
+ K0_ * c.c0 * interpolate(convolutions.influence_weight) /
+ (static_cast(dist.size()) * h);
+ }
+
+ // Remove unstable local-polynomial values.
+ for (Eigen::Index j = 0; j < fit.density.size(); ++j) {
+ if (!(fit.density(j) > 0.0) || !std::isfinite(fit.density(j)) ||
+ !std::isfinite(fit.influence_num(j))) {
+ fit.density(j) = 0.0;
+ fit.influence_num(j) = 0.0;
+ }
+ }
+ return fit;
+}
+
+//! Fuses the bulk density and influence numerators with available lower and
+//! upper boundary components. Upper components are stored by increasing
+//! distance from the upper endpoint and are therefore indexed in reverse.
+//! @param grid increasing evaluation grid on the original scale.
+//! @param w_lower lower-endpoint fusion weights on `grid`.
+//! @param w_upper upper-endpoint fusion weights on `grid`.
+//! @param lower lower boundary component; empty means unavailable.
+//! @param upper upper boundary component; empty means unavailable.
+//! @param influences on input, bulk influences; on output, fused influences.
+inline void
+Kde1d::fuse_boundary_components(const Eigen::VectorXd& grid,
+ const Eigen::VectorXd& w_lower,
+ const Eigen::VectorXd& w_upper,
+ const BoundaryComponent& lower,
+ const BoundaryComponent& upper,
+ Eigen::VectorXd& influences)
+{
+ const Eigen::VectorXd f_bulk = grid_.get_values();
+ Eigen::VectorXd f(grid.size());
+ Eigen::VectorXd influence_num(grid.size());
+ for (Eigen::Index j = 0; j < grid.size(); ++j) {
+ const bool use_lower =
+ w_lower(j) > 0.0 && j < lower.density.size();
+ const Eigen::Index upper_j = grid.size() - 1 - j;
+ const bool use_upper =
+ w_upper(j) > 0.0 && upper_j < upper.density.size();
+ const double w_bulk =
+ 1.0 - (use_lower ? w_lower(j) : 0.0) -
+ (use_upper ? w_upper(j) : 0.0);
+ f(j) = w_bulk * f_bulk(j);
+ influence_num(j) = w_bulk * f_bulk(j) * influences(j);
+ if (use_lower) {
+ f(j) += w_lower(j) * lower.density(j);
+ influence_num(j) += w_lower(j) * lower.influence_num(j);
+ }
+ if (use_upper) {
+ f(j) += w_upper(j) * upper.density(upper_j);
+ influence_num(j) += w_upper(j) * upper.influence_num(upper_j);
+ }
+ influences(j) = f(j) > 0.0 ? influence_num(j) / f(j) : 0.0;
+ }
+ grid_ = interp::InterpolationGrid(grid, f.cwiseMax(0.0), 3);
+}
+
+//! Fuses the transformed bulk density @f$\widehat f_B@f$ with endpoint
+//! densities @f$\widehat f_L@f$ and @f$\widehat f_U@f$. Each endpoint density
+//! uses a local-linear kernel with an automatically selected bandwidth, based on
+//! all distances on bounded support and the closest
+//! 75% of positive-weight observations on one-sided support. Smooth weights
+//! @f$w_L@f$ and @f$w_U@f$ are functions of the weighted bulk CDF and shrink
+//! with effective sample size @f$n_e@f$:
+//! @f[
+//! \begin{align*}
+//! q &= \min(1/4,n_e^{-1/2}), & z(p) &= \min(1,p/q), \\[2pt]
+//! g(p) &= \{1-z(p)^2\}^2, \\[2pt]
+//! w_L &= g(F_B), & w_U &= g(1-F_B), \\[2pt]
+//! \widehat f
+//! &= w_L\widehat f_L+(1-w_L-w_U)\widehat f_B+w_U\widehat f_U, \\[2pt]
+//! \nu &= w_L\nu_L+(1-w_L-w_U)\nu_B+w_U\nu_U, \\[2pt]
+//! \mathrm{influence} &= \nu/\widehat f.
+//! \end{align*}
+//! @f]
+//! Here @f$\nu@f$ is the influence numerator, i.e., the diagonal kernel
+//! contribution on the density scale. Only endpoints classified as finite are
+//! repaired.
+//! @param x observations on the original scale.
+//! @param grid increasing evaluation grid on the original scale.
+//! @param influences on input, bulk-fit influences on `grid`; on output,
+//! influences of the fused estimate.
+//! @param weights normalized case weights; empty means equal weights.
+inline void
+Kde1d::repair_boundaries(const Eigen::VectorXd& x,
+ const Eigen::VectorXd& grid,
+ Eigen::VectorXd& influences,
+ const Eigen::VectorXd& weights)
+{
+ const Eigen::VectorXd case_weights =
+ weights.size() > 0 ? weights : Eigen::VectorXd::Ones(x.size());
+ const double n_eff =
+ case_weights.sum() * case_weights.sum() / case_weights.squaredNorm();
+ if (n_eff < 16.0)
+ return;
+
+ EndpointData lower_endpoint;
+ EndpointData upper_endpoint;
+ if (!std::isnan(xmin_))
+ lower_endpoint = prepare_endpoint(x, case_weights, EndpointSide::lower);
+ if (!std::isnan(xmax_))
+ upper_endpoint = prepare_endpoint(x, case_weights, EndpointSide::upper);
+ if (!lower_endpoint.finite && !upper_endpoint.finite)
+ return;
+
+ // q is the shrinking probability width of each endpoint weight.
+ const double q = std::min(0.25, 1.0 / std::sqrt(n_eff));
+ auto endpoint_weight = [&](double probability) {
+ const double z = std::min(1.0, probability / q);
+ return (1.0 - z * z) * (1.0 - z * z);
+ };
+
+ // Only evaluate endpoint fits where their fusion weights are nonzero.
+ const Eigen::VectorXd bulk_cdf = grid_.integrate(grid, true);
+ Eigen::VectorXd w_l = Eigen::VectorXd::Zero(grid.size());
+ Eigen::VectorXd w_u = Eigen::VectorXd::Zero(grid.size());
+ for (Eigen::Index j = 0; j < grid.size(); ++j) {
+ if (lower_endpoint.finite)
+ w_l(j) = endpoint_weight(bulk_cdf(j));
+ if (upper_endpoint.finite)
+ w_u(j) = endpoint_weight(1.0 - bulk_cdf(j));
+ }
+ const Eigen::Index n_lower = (w_l.array() > 0.0).count();
+ const Eigen::Index n_upper = (w_u.array() > 0.0).count();
+
+ // Reflection invariance gives both endpoints the same boundary bandwidth.
+ const double bandwidth_fraction =
+ (!std::isnan(xmin_) && !std::isnan(xmax_)) ? 1.0 : 0.75;
+ const EndpointData& bandwidth_endpoint =
+ lower_endpoint.finite ? lower_endpoint : upper_endpoint;
+ const double h =
+ select_boundary_bandwidth(bandwidth_endpoint, bandwidth_fraction);
+
+ BoundaryComponent lower;
+ BoundaryComponent upper;
+ if (lower_endpoint.finite && n_lower > 0) {
+ lower = fit_boundary_component(lower_endpoint.dist,
+ grid.head(n_lower).array() - xmin_,
+ h,
+ lower_endpoint.weights);
+ }
+ if (upper_endpoint.finite && n_upper > 0) {
+ Eigen::VectorXd upper_eval_dist =
+ (xmax_ - grid.tail(n_upper).array()).reverse();
+ upper = fit_boundary_component(
+ upper_endpoint.dist, upper_eval_dist, h, upper_endpoint.weights);
+ }
+ if (lower.density.size() == 0 && upper.density.size() == 0)
+ return;
+ fuse_boundary_components(grid, w_l, w_u, lower, upper, influences);
+}
+
// Bandwidth for Kernel Density Estimation
//' @param x vector of observations
//' @param bandwidth bandwidth parameter, NA for automatic selection.
-//' @param multiplier bandwidth multiplieriplier.
-//' @param discrete whether a jittered estimate is computed.
+//' @param multiplier bandwidth multiplier.
//' @param weights vector of weights for each observation (can be empty).
//' @param degree polynomial degree.
//' @return the selected bandwidth
diff --git a/include/kde1d/kdefft.hpp b/include/kde1d/kdefft.hpp
index 1434486..8485027 100644
--- a/include/kde1d/kdefft.hpp
+++ b/include/kde1d/kdefft.hpp
@@ -38,7 +38,7 @@ class KdeFFT
//! @param bandwidth the bandwidth parameter.
//! @param lower lower bound of the grid.
//! @param upper bound of the grid.
-//! @param weigths optional vector of weights for each observation.
+//! @param weights optional vector of weights for each observation.
//! @param num_bins number of bins for the FFT grid.
inline KdeFFT::KdeFFT(const Eigen::VectorXd& x,
double bandwidth,
diff --git a/test/test.cpp b/test/test.cpp
index a3c76e2..f1a3ab6 100644
--- a/test/test.cpp
+++ b/test/test.cpp
@@ -119,10 +119,13 @@ TEST_CASE("boundary grids resolve the transformed support", "[boundary-grid]")
TEST_CASE("one-boundary fits are reflection equivariant", "[boundary-reflection]")
{
Eigen::VectorXd observations = Eigen::VectorXd::LinSpaced(500, 0.01, 5.0);
- Kde1d left_bounded(0.0, NAN, "continuous", 1.0, 0.5);
+ Kde1d left_bounded(0.0, NAN, "continuous");
left_bounded.fit(observations);
+ Kde1d manual(0.0, NAN, "continuous", 1.0, left_bounded.get_bandwidth());
+ manual.fit(observations);
+ CHECK(left_bounded.get_values().isApprox(manual.get_values(), 1e-12));
- Kde1d right_bounded(NAN, 0.0, "continuous", 1.0, 0.5);
+ Kde1d right_bounded(NAN, 0.0, "continuous");
right_bounded.fit(-observations);
CHECK(left_bounded.get_grid_points().isApprox(
@@ -133,6 +136,214 @@ TEST_CASE("one-boundary fits are reflection equivariant", "[boundary-reflection]
CHECK(left_bounded.get_edf() == Approx(right_bounded.get_edf()));
}
+TEST_CASE("one-sided finite endpoints use the boundary expert",
+ "[boundary-expert]")
+{
+ const Eigen::VectorXd probabilities =
+ Eigen::VectorXd::LinSpaced(200, 0.5 / 200.0, 1.0 - 0.5 / 200.0);
+ const Eigen::VectorXd observations = (-(1.0 - probabilities.array()).log());
+ Kde1d fit(0.0, NAN, "continuous");
+ fit.fit(observations);
+
+ Kde1d manual(0.0, NAN, "continuous", 1.0, fit.get_bandwidth());
+ manual.fit(observations);
+ CHECK(fit.get_values().isApprox(manual.get_values(), 1e-12));
+
+ Eigen::VectorXd expected_density(6);
+ expected_density << 0.9797951, 0.9798032, 0.9780643, 0.9448344, 0.7316556,
+ 0.0077380;
+ Eigen::VectorXd selected_density(6);
+ selected_density << fit.get_values()(0), fit.get_values()(4),
+ fit.get_values()(24), fit.get_values()(49), fit.get_values()(99),
+ fit.get_values()(199);
+ CHECK(selected_density.isApprox(expected_density, 2e-4));
+ CHECK(std::isfinite(fit.get_edf()));
+
+ Kde1d reflected(NAN, 0.0, "continuous");
+ reflected.fit(-observations);
+ CHECK(fit.get_grid_points().isApprox(-reflected.get_grid_points().reverse(),
+ 1e-12));
+ CHECK(fit.get_values().isApprox(reflected.get_values().reverse(), 1e-12));
+ CHECK(fit.get_edf() == Approx(reflected.get_edf()).epsilon(1e-12));
+}
+
+TEST_CASE("boundary experts can be disabled", "[boundary-expert]")
+{
+ const Eigen::VectorXd probabilities =
+ Eigen::VectorXd::LinSpaced(200, 0.5 / 200.0, 1.0 - 0.5 / 200.0);
+ const Eigen::VectorXd observations = (-(1.0 - probabilities.array()).log());
+
+ Kde1d default_fit(0.0, NAN, "continuous");
+ default_fit.fit(observations);
+ Kde1d enabled(0.0, NAN, "continuous", 1.0, NAN, 2, 400, true);
+ enabled.fit(observations);
+ Kde1d disabled(0.0, NAN, "continuous", 1.0, NAN, 2, 400, false);
+ disabled.fit(observations);
+ Kde1d disabled_manual(
+ 0.0, NAN, "continuous", 1.0, disabled.get_bandwidth(), 2, 400, false);
+ disabled_manual.fit(observations);
+
+ CHECK(default_fit.get_boundary_repair());
+ CHECK(enabled.get_boundary_repair());
+ CHECK_FALSE(disabled.get_boundary_repair());
+ CHECK(default_fit.get_values().isApprox(enabled.get_values(), 1e-12));
+ CHECK_FALSE(default_fit.get_values().isApprox(disabled.get_values(), 1e-6));
+ CHECK(disabled.get_values().isApprox(disabled_manual.get_values(), 1e-12));
+ CHECK(disabled.str().find("boundary_repair=false") != std::string::npos);
+}
+
+TEST_CASE("one-sided vanishing endpoints retain the bulk fit",
+ "[boundary-expert]")
+{
+ const Eigen::VectorXd probabilities =
+ Eigen::VectorXd::LinSpaced(200, 0.5 / 200.0, 1.0 - 0.5 / 200.0);
+ const Eigen::VectorXd observations = probabilities.array().sqrt();
+ Kde1d fit(0.0, NAN, "continuous");
+ fit.fit(observations);
+
+ Kde1d bulk(0.0, NAN, "continuous", 1.0, fit.get_bandwidth());
+ bulk.fit(observations);
+ CHECK(fit.get_grid_points().isApprox(bulk.get_grid_points(), 1e-12));
+ CHECK(fit.get_values().isApprox(bulk.get_values(), 1e-12));
+ CHECK(fit.get_edf() == Approx(bulk.get_edf()).epsilon(1e-12));
+}
+
+TEST_CASE("two-sided finite endpoints use the boundary experts",
+ "[boundary-expert]")
+{
+ const Eigen::VectorXd observations =
+ Eigen::VectorXd::LinSpaced(200, 0.5 / 200.0, 1.0 - 0.5 / 200.0);
+ Kde1d fit(0.0, 1.0, "continuous");
+ fit.fit(observations);
+
+ Kde1d manual(0.0, 1.0, "continuous", 1.0, fit.get_bandwidth());
+ manual.fit(observations);
+ CHECK(fit.get_values().isApprox(manual.get_values(), 1e-12));
+
+ Eigen::VectorXd expected_density(6);
+ expected_density << 0.9993158, 0.9992803, 0.9985546, 1.0001208, 0.9985473,
+ 0.9993158;
+ Eigen::VectorXd selected_density(6);
+ selected_density << fit.get_values()(0), fit.get_values()(49),
+ fit.get_values()(99), fit.get_values()(199), fit.get_values()(299),
+ fit.get_values()(400);
+ CHECK(selected_density.isApprox(expected_density, 2e-4));
+ CHECK(std::isfinite(fit.get_edf()));
+
+ Kde1d scaled(-3.0, 7.0, "continuous");
+ scaled.fit(-3.0 + 10.0 * observations.array());
+ CHECK(((scaled.get_grid_points().array() + 3.0) / 10.0)
+ .matrix()
+ .isApprox(fit.get_grid_points(), 1e-12));
+ CHECK((10.0 * scaled.get_values()).isApprox(fit.get_values(), 1e-10));
+ CHECK(scaled.get_edf() == Approx(fit.get_edf()).epsilon(1e-10));
+}
+
+TEST_CASE("two-sided endpoints are classified independently",
+ "[boundary-expert]")
+{
+ const Eigen::VectorXd probabilities =
+ Eigen::VectorXd::LinSpaced(200, 0.5 / 200.0, 1.0 - 0.5 / 200.0);
+ const Eigen::VectorXd observations =
+ 1.0 - (1.0 - probabilities.array()).sqrt();
+ Kde1d fit(0.0, 1.0, "continuous");
+ fit.fit(observations);
+
+ Eigen::VectorXd expected_density(6);
+ expected_density << 1.9985571, 1.9950070, 1.9824915, 0.9946390, 0.0562128,
+ 0.0000020;
+ Eigen::VectorXd selected_density(6);
+ selected_density << fit.get_values()(0), fit.get_values()(49),
+ fit.get_values()(99), fit.get_values()(199), fit.get_values()(299),
+ fit.get_values()(400);
+ CHECK(selected_density.isApprox(expected_density, 2e-4));
+
+ Kde1d reflected(0.0, 1.0, "continuous");
+ reflected.fit(1.0 - observations.array());
+ CHECK(fit.get_values().isApprox(reflected.get_values().reverse(), 1e-12));
+ CHECK(fit.get_edf() == Approx(reflected.get_edf()).epsilon(1e-12));
+}
+
+TEST_CASE("two-sided exploding endpoints retain the bulk fit",
+ "[boundary-expert]")
+{
+ const Eigen::VectorXd probabilities =
+ Eigen::VectorXd::LinSpaced(200, 0.5 / 200.0, 1.0 - 0.5 / 200.0);
+ const Eigen::VectorXd observations =
+ (0.5 * std::acos(-1.0) * probabilities.array()).sin().square();
+ Kde1d fit(0.0, 1.0, "continuous");
+ fit.fit(observations);
+
+ Kde1d bulk(0.0, 1.0, "continuous", 1.0, fit.get_bandwidth());
+ bulk.fit(observations);
+ CHECK(fit.get_values().isApprox(bulk.get_values(), 1e-12));
+ CHECK(fit.get_edf() == Approx(bulk.get_edf()).epsilon(1e-12));
+}
+
+TEST_CASE("boundary experts support weights and manual bandwidths",
+ "[boundary-expert]")
+{
+ const Eigen::VectorXd probabilities =
+ Eigen::VectorXd::LinSpaced(200, 0.5 / 200.0, 1.0 - 0.5 / 200.0);
+ const Eigen::VectorXd observations = (-(1.0 - probabilities.array()).log());
+ const Eigen::VectorXd weights =
+ Eigen::VectorXd::LinSpaced(observations.size(), 0.5, 1.5);
+
+ Kde1d weighted(0.0, NAN, "continuous");
+ weighted.fit(observations, weights);
+ Kde1d weighted_manual(0.0, NAN, "continuous", 1.0, weighted.get_bandwidth());
+ weighted_manual.fit(observations, weights);
+ CHECK(weighted.get_values().isApprox(weighted_manual.get_values(), 1e-12));
+ CHECK(weighted.get_edf() == Approx(weighted_manual.get_edf()).epsilon(1e-12));
+
+ Kde1d rescaled_weights(0.0, NAN, "continuous");
+ rescaled_weights.fit(observations, 7.0 * weights);
+ CHECK(weighted.get_values().isApprox(rescaled_weights.get_values(), 1e-12));
+ CHECK(weighted.get_edf() ==
+ Approx(rescaled_weights.get_edf()).epsilon(1e-12));
+ CHECK(weighted.get_loglik() ==
+ Approx(rescaled_weights.get_loglik()).epsilon(1e-12));
+
+ Eigen::VectorXd zero_weights = weights;
+ for (Eigen::Index i = 1; i < zero_weights.size(); i += 5)
+ zero_weights(i) = 0.0;
+ Eigen::VectorXd retained_observations(160);
+ Eigen::VectorXd retained_weights(160);
+ Eigen::Index retained_index = 0;
+ for (Eigen::Index i = 0; i < zero_weights.size(); ++i) {
+ if (zero_weights(i) > 0.0) {
+ retained_observations(retained_index) = observations(i);
+ retained_weights(retained_index++) = zero_weights(i);
+ }
+ }
+ Kde1d with_zeros(0.0, NAN, "continuous");
+ with_zeros.fit(observations, zero_weights);
+ Kde1d without_zeros(0.0, NAN, "continuous");
+ without_zeros.fit(retained_observations, retained_weights);
+ CHECK(with_zeros.get_values().isApprox(without_zeros.get_values(), 1e-12));
+ CHECK(with_zeros.get_edf() == Approx(without_zeros.get_edf()).epsilon(1e-12));
+
+ Kde1d reflected(NAN, 0.0, "continuous");
+ reflected.fit(-observations, weights);
+ CHECK(
+ weighted.get_values().isApprox(reflected.get_values().reverse(), 1e-12));
+ CHECK(weighted.get_edf() == Approx(reflected.get_edf()).epsilon(1e-12));
+
+ Kde1d linear(0.0, NAN, "continuous", 1.0, NAN, 1);
+ linear.fit(observations);
+ Kde1d linear_bulk(0.0, NAN, "continuous", 1.0, linear.get_bandwidth(), 1);
+ linear_bulk.fit(observations);
+ CHECK(linear.get_values().isApprox(linear_bulk.get_values(), 1e-12));
+ CHECK(linear.get_edf() == Approx(linear_bulk.get_edf()).epsilon(1e-12));
+
+ Kde1d small(0.0, NAN, "continuous");
+ small.fit(observations.head(15));
+ Kde1d small_bulk(0.0, NAN, "continuous", 1.0, small.get_bandwidth());
+ small_bulk.fit(observations.head(15));
+ CHECK(small.get_values().isApprox(small_bulk.get_values(), 1e-12));
+ CHECK(small.get_edf() == Approx(small_bulk.get_edf()).epsilon(1e-12));
+}
+
TEST_CASE("finite support truncates density", "[finite-support]")
{
Eigen::VectorXd observations = Eigen::VectorXd::LinSpaced(200, 0.0, 1.0);
@@ -439,6 +650,7 @@ TEST_CASE("bandwidth selection on refit", "[bandwidth][refit]")
CHECK(std::isnan(from_grid.get_bandwidth()));
CHECK(from_grid.get_multiplier() == Approx(1.0));
CHECK(from_grid.get_degree() == 2);
+ CHECK(from_grid.get_boundary_repair());
}
}
diff --git a/test/test_numerical_invariants.cpp b/test/test_numerical_invariants.cpp
index c726e9d..a29345f 100644
--- a/test/test_numerical_invariants.cpp
+++ b/test/test_numerical_invariants.cpp
@@ -85,20 +85,23 @@ TEST_CASE("finite-support fits are scale equivariant",
Eigen::VectorXd points(10);
points << 0.0, 1e-10, 1e-8, 1e-6, 1e-4, 1e-2, 0.1, 0.5, 0.9, 1.0;
- kde1d::Kde1d reference(0.0, 1.0, "continuous", 1.0, 0.3, 2, 400);
- reference.fit(observations);
-
- for (double scale : { 1e-4, 1e4 }) {
- INFO("scale=" << scale);
- kde1d::Kde1d scaled(
- 0.0, scale, "continuous", 1.0, 0.3, 2, 400);
- scaled.fit(observations * scale);
-
- CHECK((scaled.get_grid_points() / scale)
- .isApprox(reference.get_grid_points(), 1e-11));
- CHECK((scaled.pdf(points * scale) * scale)
- .isApprox(reference.pdf(points), 1e-10));
- CHECK(scaled.cdf(points * scale).isApprox(reference.cdf(points), 1e-10));
+ for (size_t degree = 0; degree < 3; ++degree) {
+ kde1d::Kde1d reference(
+ 0.0, 1.0, "continuous", 1.0, 0.3, degree, 400);
+ reference.fit(observations);
+
+ for (double scale : { 1e-4, 1e4 }) {
+ INFO("degree=" << degree << ", scale=" << scale);
+ kde1d::Kde1d scaled(
+ 0.0, scale, "continuous", 1.0, 0.3, degree, 400);
+ scaled.fit(observations * scale);
+
+ CHECK((scaled.get_grid_points() / scale)
+ .isApprox(reference.get_grid_points(), 1e-11));
+ CHECK((scaled.pdf(points * scale) * scale)
+ .isApprox(reference.pdf(points), 1e-10));
+ CHECK(scaled.cdf(points * scale).isApprox(reference.cdf(points), 1e-10));
+ }
}
}
@@ -111,23 +114,42 @@ TEST_CASE("one-sided fits are scale equivariant",
Eigen::VectorXd points(10);
points << 0.0, 1e-10, 1e-8, 1e-6, 1e-4, 1e-2, 0.1, 0.5, 1.0, 5.0;
- kde1d::Kde1d reference(0.0, NAN, "continuous", 1.0, 0.3, 2, 400);
- reference.fit(observations);
-
- for (double scale : { 1e-4, 1e4 }) {
- INFO("scale=" << scale);
- kde1d::Kde1d scaled(
- 0.0, NAN, "continuous", 1.0, 0.3, 2, 400);
- scaled.fit(observations * scale);
-
- CHECK((scaled.get_grid_points() / scale)
- .isApprox(reference.get_grid_points(), 1e-11));
- CHECK((scaled.pdf(points * scale) * scale)
- .isApprox(reference.pdf(points), 1e-10));
- CHECK(scaled.cdf(points * scale).isApprox(reference.cdf(points), 1e-10));
+ for (size_t degree = 0; degree < 3; ++degree) {
+ kde1d::Kde1d reference(
+ 0.0, NAN, "continuous", 1.0, 0.3, degree, 400);
+ reference.fit(observations);
+
+ for (double scale : { 1e-4, 1e4 }) {
+ INFO("degree=" << degree << ", scale=" << scale);
+ kde1d::Kde1d scaled(
+ 0.0, NAN, "continuous", 1.0, 0.3, degree, 400);
+ scaled.fit(observations * scale);
+
+ CHECK((scaled.get_grid_points() / scale)
+ .isApprox(reference.get_grid_points(), 1e-11));
+ CHECK((scaled.pdf(points * scale) * scale)
+ .isApprox(reference.pdf(points), 1e-10));
+ CHECK(scaled.cdf(points * scale).isApprox(reference.cdf(points), 1e-10));
+ }
}
}
+TEST_CASE("one-sided power grids do not explode on exponential data",
+ "[numerical-invariants][boundary-grid]")
+{
+ Eigen::VectorXd probabilities =
+ Eigen::VectorXd::LinSpaced(500, 0.5 / 500.0, 1.0 - 0.5 / 500.0);
+ Eigen::VectorXd observations = (-(1.0 - probabilities.array()).log());
+ kde1d::Kde1d fit(0.0, NAN, "continuous");
+ fit.fit(observations);
+
+ Eigen::VectorXd boundary_density =
+ fit.pdf(Eigen::VectorXd::Constant(1, 0.0));
+ CHECK(fit.get_grid_points().maxCoeff() < 10.0 * observations.maxCoeff());
+ CHECK(std::isfinite(boundary_density(0)));
+ CHECK(boundary_density(0) > 0.0);
+}
+
TEST_CASE("interpolation preserves nonuniform-grid reference values",
"[interpolation][parity]")
{