ox, oy;
+ if ((int) hx.size() == 1) {
+ for (int k = 0; k < 3; ++k) { ox.push_back(hx[0]); oy.push_back(hy[0]); }
+ } else {
+ ox = hx; oy = hy;
+ }
+ for (size_t i = 0; i < ox.size(); ++i) {
+ ox[i] = std::min(std::max(ox[i], minx), maxx);
+ oy[i] = std::min(std::max(oy[i], miny), maxy);
+ }
+
+ return DataFrame::create(_["x"] = wrap(ox), _["y"] = wrap(oy));
+}
diff --git a/tools/NNS/src/partial_moments.cpp b/tools/NNS/src/partial_moments.cpp
index 2bc723e..40eb2d2 100644
--- a/tools/NNS/src/partial_moments.cpp
+++ b/tools/NNS/src/partial_moments.cpp
@@ -1027,4 +1027,4 @@ List PMMatrix_CPv(
Named("cov.matrix") = covMat
)
);
-}
+}
\ No newline at end of file
diff --git a/tools/NNS/tests/testthat/Rplots.pdf b/tools/NNS/tests/testthat/Rplots.pdf
index 4646286..d3be77e 100644
Binary files a/tools/NNS/tests/testthat/Rplots.pdf and b/tools/NNS/tests/testthat/Rplots.pdf differ
diff --git a/tools/NNS/tests/testthat/test_NNS_reg_infinite.R b/tools/NNS/tests/testthat/test_NNS_reg_infinite.R
new file mode 100644
index 0000000..d75eb22
--- /dev/null
+++ b/tools/NNS/tests/testthat/test_NNS_reg_infinite.R
@@ -0,0 +1,12 @@
+test_that("native partition complete-case handling keeps infinities", {
+ part <- NNS.part(
+ x = c(Inf, Inf, NA_real_, NaN),
+ y = c(1, 2, 3, 4),
+ order = 1,
+ obs.req = 0,
+ noise.reduction = "median"
+ )
+
+ expect_true(any(is.infinite(part$regression.points$x)))
+ expect_false(any(is.na(part$regression.points$x)))
+})
diff --git a/tools/NNS/vignettes/NNSvignette_01_Overview.R b/tools/NNS/vignettes/NNSvignette_01_Overview.R
new file mode 100644
index 0000000..9fb1efc
--- /dev/null
+++ b/tools/NNS/vignettes/NNSvignette_01_Overview.R
@@ -0,0 +1,205 @@
+## ----setup, message=FALSE-----------------------------------------------------
+# Prereqs (uncomment if needed):
+# install.packages("NNS")
+# install.packages(c("data.table","xts","zoo","Rfast"))
+
+library(NNS)
+library(data.table)
+
+## ----include=FALSE, message=FALSE---------------------------------------------
+data.table::setDTthreads(1L)
+options(mc.cores = 1)
+RcppParallel::setThreadOptions(numThreads = 1)
+Sys.setenv("OMP_THREAD_LIMIT" = 1)
+
+## -----------------------------------------------------------------------------
+set.seed(42)
+
+# Normal sample
+y <- rnorm(3000)
+mu <- mean(y)
+L2 <- LPM(2, mu, y); U2 <- UPM(2, mu, y)
+cat(sprintf("LPM2 + UPM2 = %.6f vs var(y)=%.6f\n", (L2+U2)*(length(y) / (length(y) - 1)), var(y)))
+
+# Empirical CDF via LPM.ratio(0, t, x)
+for (t in c(-1,0,1)) {
+ cdf_lpm <- LPM.ratio(0, t, y)
+ cat(sprintf("CDF at t=%+.1f : LPM.ratio=%.4f | empirical=%.4f\n", t, cdf_lpm, mean(y<=t)))
+}
+
+# Asymmetry on a skewed distribution
+z <- rexp(3000)-1; mu_z <- mean(z)
+cat(sprintf("Skewed z: LPM2=%.4f, UPM2=%.4f (expect imbalance)\n", LPM(2,mu_z,z), UPM(2,mu_z,z)))
+
+## -----------------------------------------------------------------------------
+M <- NNS.moments(y)
+M
+
+## -----------------------------------------------------------------------------
+set.seed(23)
+multimodal <- c(rnorm(1500,-2,.5), rnorm(1500,2,.5))
+NNS.mode(multimodal,multi = TRUE)
+
+## -----------------------------------------------------------------------------
+qgrid <- LPM.VaR(seq(0.05,0.95,.1),0,z) # equivalent to quantile(z,probs = seq(0.05,0.95,by=0.1))
+CDF_tbl <- data.table(threshold = as.numeric(qgrid), CDF = LPM.ratio(0,qgrid,z))
+CDF_tbl
+
+## -----------------------------------------------------------------------------
+set.seed(1)
+x <- runif(2000,-1,1)
+y <- x^2 + rnorm(2000, sd=.05)
+cat(sprintf("Pearson r = %.4f\n", cor(x,y)))
+cat(sprintf("NNS.dep = %.4f\n", NNS.dep(x,y)$Dependence))
+
+X <- data.frame(a=x, b=y, c=x*y + rnorm(2000, sd=.05))
+pm <- PM.matrix(1, 1, target = "means", variable=X, pop_adj=TRUE)
+pm
+
+cop <- NNS.copula(X, continuous=TRUE, plot=FALSE)
+cop
+
+## ----eval=FALSE---------------------------------------------------------------
+# # Data
+# set.seed(123); x = rnorm(100); y = rnorm(100); z = expand.grid(x, y)
+#
+# # Plot
+# rgl::plot3d(z[,1], z[,2], Co.LPM(0, z[,1], z[,2], z[,1], z[,2]), col = "red")
+#
+# # Uniform values
+# u_x = LPM.ratio(0, x, x); u_y = LPM.ratio(0, y, y); z = expand.grid(u_x, u_y)
+#
+# # Plot
+# rgl::plot3d(z[,1], z[,2], Co.LPM(0, z[,1], z[,2], z[,1], z[,2]), col = "blue")
+
+## -----------------------------------------------------------------------------
+A <- rnorm(100, mean = 0, sd = 1)
+B <- rnorm(100, mean = 0, sd = 5)
+C <- rnorm(100, mean = 10, sd = 1)
+D <- rnorm(100, mean = 10, sd = 10)
+
+X <- data.frame(A, B, C, D)
+
+# Linear scaling
+lin_norm <- NNS.norm(X, linear = TRUE, chart.type=NULL, location=NULL)
+
+## -----------------------------------------------------------------------------
+px <- 100 + cumsum(rnorm(260, sd = 1))
+rn <- NNS.rescale(px, a=100, b=0.03, method="riskneutral", T=1, type="Terminal")
+c( target = 100*exp(0.03*1), mean_rn = mean(rn) )
+
+## -----------------------------------------------------------------------------
+ctrl <- rnorm(200, 0, 1)
+trt <- rnorm(180, 0.35, 1.2)
+NNS.ANOVA(control=ctrl, treatment=trt, means.only=FALSE, plot=FALSE)
+
+A <- list(g1=rnorm(150,0.0,1.1), g2=rnorm(150,0.2,1.0), g3=rnorm(150,-0.1,0.9))
+NNS.ANOVA(control=A, means.only=TRUE, plot=FALSE)
+
+## ----stochsuperiority, echo=TRUE----------------------------------------------
+set.seed(123)
+x = rnorm(1000, mean = 0, sd = 1)
+y = rnorm(1000, mean = 1, sd = 1)
+
+NNS.SS(x, y)
+
+## ----stochsuperiorityci, echo=TRUE, eval=FALSE--------------------------------
+# NNS.SS(x, y, confidence.interval = TRUE, reps = 999, ci = 0.95)[1:5]
+#
+# $p_gt
+# [1] 0.233915
+#
+# $p_tie
+# [1] 0
+#
+# $p_star
+# [1] 0.233915
+#
+# $lower
+# [1] 0.2105631
+#
+# $upper
+# [1] 0.2537789
+
+## ----stochsuperioritydiscrete, echo=TRUE--------------------------------------
+set.seed(123)
+x = sample(1:5, 100, replace = TRUE)
+y = sample(1:5, 100, replace = TRUE)
+
+NNS.SS(x, y)
+
+## ----fig.width=7, fig.height=5, fig.align='center'----------------------------
+# Example 1: Nonlinear regression
+set.seed(123)
+x_train <- runif(1000, -2, 2)
+y_train <- sin(pi * x_train) + rnorm(1000, sd = 0.2)
+
+x_test <- seq(-2, 2, length.out = 100)
+
+NNS.reg(x = x_train, y = y_train, order = NULL, point.est = x_test)
+
+## ----eval = FALSE-------------------------------------------------------------
+# # Simple train/test for boosting & stacking
+# test.set = 141:150
+#
+# boost <- NNS.boost(IVs.train = iris[-test.set, 1:4],
+# DV.train = iris[-test.set, 5],
+# IVs.test = iris[test.set, 1:4],
+# epochs = 10, learner.trials = 10,
+# status = FALSE, balance = TRUE,
+# type = "CLASS", folds = 5)
+#
+#
+# mean(boost$results == as.numeric(iris[test.set,5]))
+# # [1] 1
+#
+#
+# boost$feature.weights; boost$feature.frequency
+#
+# stacked <- NNS.stack(IVs.train = iris[-test.set, 1:4],
+# DV.train = iris[-test.set, 5],
+# IVs.test = iris[test.set, 1:4],
+# type = "CLASS", balance = TRUE,
+# ncores = 1, folds = 1)
+# mean(stacked$stack == as.numeric(iris[test.set,5]))
+# # [1] 1
+
+## -----------------------------------------------------------------------------
+NNS.caus(mtcars$hp, mtcars$mpg) # hp -> mpg
+NNS.caus(mtcars$mpg, mtcars$hp) # hp -> mpg
+
+## ----fig.width=7, fig.align='center'------------------------------------------
+# Univariate nonlinear ARMA
+z <- as.numeric(scale(sin(1:480/8) + rnorm(480, sd=.35)))
+
+# Seasonality detection (prints a summary)
+seasonal_period <- NNS.seas(z, plot = FALSE)
+head(seasonal_period$all.periods)
+
+# Validate seasonal periods
+NNS.ARMA.optim(z, h = 48, seasonal.factor = seasonal_period$periods, plot = TRUE, ncores = 1)
+
+## -----------------------------------------------------------------------------
+x_ts <- cumsum(rnorm(350, sd=.7))
+mb <- NNS.meboot(x_ts, reps=5, rho = 1)
+dim(mb["replicates", ]$replicates)
+
+## -----------------------------------------------------------------------------
+mc <- NNS.MC(x_ts, reps=5, lower_rho=-1, upper_rho=1, by=.5, exp=1)
+length(mc$ensemble); names(mc$replicates)
+
+head(mc$replicates$`rho = 0`)
+
+## -----------------------------------------------------------------------------
+RA <- rnorm(240, 0.005, 0.03)
+RB <- rnorm(240, 0.003, 0.02)
+RC <- rnorm(240, 0.006, 0.04)
+
+NNS.FSD.uni(RA, RB)
+NNS.SSD.uni(RA, RB)
+NNS.TSD.uni(RA, RB)
+
+Rmat <- cbind(A=RA, B=RB, C=RC)
+try(NNS.SD.cluster(Rmat, degree = 1))
+try(NNS.SD.efficient.set(Rmat, degree = 1))
+
diff --git a/tools/NNS/vignettes/NNSvignette_01_Overview.html b/tools/NNS/vignettes/NNSvignette_01_Overview.html
new file mode 100644
index 0000000..56f9d4a
--- /dev/null
+++ b/tools/NNS/vignettes/NNSvignette_01_Overview.html
@@ -0,0 +1,1362 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+Getting Started with NNS: Overview
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+Getting Started with NNS: Overview
+Fred Viole
+
+
+
+# Prereqs (uncomment if needed):
+# install.packages("NNS")
+# install.packages(c("data.table","xts","zoo","Rfast"))
+
+library(NNS)
+library(data.table)
+
Orientation
+
Goal. A complete, hands‑on curriculum for Nonlinear
+Nonparametric Statistics (NNS) using partial moments.
+Each section blends narrative intuition, precise math, and executable
+code.
+
Structure. 1. Foundations — partial moments &
+variance decomposition 2. Descriptive & distributional tools 3.
+Dependence & nonlinear association 4. Normalization & Rescaling
+5. Hypothesis testing, ANOVA & Stochastic Superiority 6. Regression,
+boosting, stacking & causality 7. Time series & forecasting 8.
+Simulation (max‑entropy) & Monte Carlo 9. Portfolio & stochastic
+dominance
+
Notation. For a random variable \(X\) and threshold/target \(t\), the population \(n\)‑th partial moments are
+defined as:
+
\[
+\operatorname{LPM}(n,t,X)
+= \int_{-\infty}^{t} (t-x)^{n} \, dF_X(x),
+\qquad
+\operatorname{UPM}(n,t,X)
+= \int_{t}^{\infty} (x-t)^{n} \, dF_X(x).
+\]
+
The empirical estimators replace \(F_X\) with the empirical CDF \(\hat F_n\) (or, equivalently, use indicator
+functions):
+
\[
+\widehat{\operatorname{LPM}}_n(t;X) = \frac{1}{n} \sum_{i=1}^n (t-x_i)^n
+\, \mathbf{1}_{\{x_i \le t\}},
+\qquad
+\widehat{\operatorname{UPM}}_n(t;X) = \frac{1}{n} \sum_{i=1}^n (x_i-t)^n
+\, \mathbf{1}_{\{x_i > t\}}.
+\]
+
These correspond to integrals over the measurable subsets \(\{X \le t\}\) and \(\{X > t\}\) in a \(\sigma\)‑algebra; the empirical sums are
+discrete analogues of Lebesgue integrals.
+
+
+
1. Foundations — Partial Moments & Variance Decomposition
+
+
1.1 Why partial moments
+
+- Classical variance treats upside and downside symmetrically. Partial
+moments separate them, allowing asymmetric risk/reward
+analysis around a chosen target \(t\)
+(often the mean or a benchmark).
+- At \(t=\mu_X\): \[
+\operatorname{Var}(X) = \operatorname{UPM}(2,\mu_X,X) +
+\operatorname{LPM}(2,\mu_X,X)\quad\text{(exact empirical identity)}.
+\] This is not the same as splitting conditional
+variances around a threshold; partial moments use a global
+reference, preserving the between‑group contribution.
+
+
+
+
1.2 Core functions and headers
+
+LPM(degree, target, variable)UPM(degree, target, variable)
+
+
+
1.3 Code: variance decomposition & CDF
+
set.seed(42)
+
+# Normal sample
+y <- rnorm(3000)
+mu <- mean(y)
+L2 <- LPM(2, mu, y); U2 <- UPM(2, mu, y)
+cat(sprintf("LPM2 + UPM2 = %.6f vs var(y)=%.6f\n", (L2+U2)*(length(y) / (length(y) - 1)), var(y)))
+
## LPM2 + UPM2 = 1.011889 vs var(y)=1.011889
+
# Empirical CDF via LPM.ratio(0, t, x)
+for (t in c(-1,0,1)) {
+ cdf_lpm <- LPM.ratio(0, t, y)
+ cat(sprintf("CDF at t=%+.1f : LPM.ratio=%.4f | empirical=%.4f\n", t, cdf_lpm, mean(y<=t)))
+}
+
## CDF at t=-1.0 : LPM.ratio=0.1633 | empirical=0.1633
+## CDF at t=+0.0 : LPM.ratio=0.5043 | empirical=0.5043
+## CDF at t=+1.0 : LPM.ratio=0.8480 | empirical=0.8480
+
# Asymmetry on a skewed distribution
+z <- rexp(3000)-1; mu_z <- mean(z)
+cat(sprintf("Skewed z: LPM2=%.4f, UPM2=%.4f (expect imbalance)\n", LPM(2,mu_z,z), UPM(2,mu_z,z)))
+
## Skewed z: LPM2=0.2780, UPM2=0.7682 (expect imbalance)
+
Interpretation. The equality
+LPM2 + UPM2 == var(x) (Bessel adjustment used) holds
+because deviations are measured against the global mean.
+LPM.ratio(0, t, x) constructs an empirical CDF directly
+from partial‑moment counts.
+
+
+
+
3. Dependence & Nonlinear Association
+
+
3.1 Why move beyond Pearson \(r\)
+
Pearson captures linear monotone relationships. Many structures
+(U‑shapes, saturation, asymmetric tails) produce near‑zero \(r\) despite strong dependence.
+Partial‑moment dependence metrics respond to such structure.
+
Headers.
+
+Co.LPM(degree_lpm, x, y, target_x, target_y, degree_y)
+/ Co.UPM(...) (co‑partial moments)PM.matrix(LPM_degree, UPM_degree, target=NULL, variable, pop_adj=TRUE)NNS.dep(x, y) (scalar dependence coefficient)NNS.copula(X, target=NULL, continuous=TRUE, plot=FALSE, independence.overlay=FALSE)
+
+
+
3.2 Code: nonlinear dependence
+
set.seed(1)
+x <- runif(2000,-1,1)
+y <- x^2 + rnorm(2000, sd=.05)
+cat(sprintf("Pearson r = %.4f\n", cor(x,y)))
+
## Pearson r = 0.0006
+
cat(sprintf("NNS.dep = %.4f\n", NNS.dep(x,y)$Dependence))
+
## NNS.dep = 0.7097
+
X <- data.frame(a=x, b=y, c=x*y + rnorm(2000, sd=.05))
+pm <- PM.matrix(1, 1, target = "means", variable=X, pop_adj=TRUE)
+pm
+
## $cupm
+## a b c
+## a 0.17384174 0.05668152 0.10450858
+## b 0.05668152 0.05566363 0.04414923
+## c 0.10450858 0.04414923 0.07529373
+##
+## $dupm
+## a b c
+## a 0.0000000000 0.05675501 0.0005598221
+## b 0.0143108307 0.00000000 0.0036839026
+## c 0.0004239566 0.04430691 0.0000000000
+##
+## $dlpm
+## a b c
+## a 0.0000000000 0.014310831 0.0004239566
+## b 0.0567550147 0.000000000 0.0443069142
+## c 0.0005598221 0.003683903 0.0000000000
+##
+## $clpm
+## a b c
+## a 0.16803827 0.014485430 0.102709867
+## b 0.01448543 0.037120650 0.003051617
+## c 0.10270987 0.003051617 0.074865823
+##
+## $cov.matrix
+## a b c
+## a 0.3418800141 0.0001011068 0.206234664
+## b 0.0001011068 0.0927842833 -0.000789973
+## c 0.2062346637 -0.0007899730 0.150159552
+
cop <- NNS.copula(X, continuous=TRUE, plot=FALSE)
+cop
+
## [1] 0.5692785
+
+
+
3.3 Code: copula
+
# Data
+set.seed(123); x = rnorm(100); y = rnorm(100); z = expand.grid(x, y)
+
+# Plot
+rgl::plot3d(z[,1], z[,2], Co.LPM(0, z[,1], z[,2], z[,1], z[,2]), col = "red")
+
+# Uniform values
+u_x = LPM.ratio(0, x, x); u_y = LPM.ratio(0, y, y); z = expand.grid(u_x, u_y)
+
+# Plot
+rgl::plot3d(z[,1], z[,2], Co.LPM(0, z[,1], z[,2], z[,1], z[,2]), col = "blue")
+
Interpretation. NNS.dep remains high
+for curved relationships; PM.matrix collects co‑partial
+moments across variables; NNS.copula summarizes
+higher‑dimensional dependence using partial‑moment ratios. Copulas are
+returned and evaluated via Co.LPM functions.
+
+
+
+
4. Normalization and Rescaling
+
NNS provides two main tools for scaling data while preserving rank
+structure and distributional shape. Both operate via deterministic
+affine transformations.
+
+
4.1 Normalization
+
NNS.norm() rescales variables to a common magnitude
+while preserving distributional structure. The method can be
+linear (all variables forced to have the same mean) or
+nonlinear (using dependence weights to produce a more
+nuanced scaling). In the nonlinear case, the degree of association
+between variables influences the final normalized values.
+
Header.
+
+NNS.norm(x, linear=TRUE, chart.type = NULL)
+
A <- rnorm(100, mean = 0, sd = 1)
+B <- rnorm(100, mean = 0, sd = 5)
+C <- rnorm(100, mean = 10, sd = 1)
+D <- rnorm(100, mean = 10, sd = 10)
+
+X <- data.frame(A, B, C, D)
+
+# Linear scaling
+lin_norm <- NNS.norm(X, linear = TRUE, chart.type=NULL, location=NULL)
+
Interpretation. NNS.norm() brings
+variables to a common scale without distorting their distributional
+shape. Linear mode equalizes means; nonlinear mode additionally weights
+each variable by its dependence with others, so more correlated
+variables exert greater influence on the final scaling.
+
+
+
4.2 Risk‑neutral rescale (pricing context)
+
NNS.rescale() performs one‑dimensional affine
+transformations.
+
Header.
+
+NNS.rescale(x, a, b, method=c("minmax","riskneutral"), T=NULL, type=c("Terminal","Discounted"))
+
px <- 100 + cumsum(rnorm(260, sd = 1))
+rn <- NNS.rescale(px, a=100, b=0.03, method="riskneutral", T=1, type="Terminal")
+c( target = 100*exp(0.03*1), mean_rn = mean(rn) )
+
## target mean_rn
+## 103.0455 103.0455
+
Interpretation. riskneutral shifts the
+mean to match \(S_0 e^{rT}\) (Terminal)
+or \(S_0\) (Discounted), preserving
+distributional shape.
+
+
+
+
5. Hypothesis Testing, ANOVA & Stochastic Superiority
+
+
5.1 Concept
+
Instead of distributional assumptions, compare groups via
+LPM‑based CDFs. Output is a degree of
+certainty (not a p‑value) for equality of populations or means.
+
Header.
+
+NNS.ANOVA(control, treatment, means.only=FALSE, medians=FALSE, confidence.interval=.95, tails=c("Both","left","right"), pairwise=FALSE, plot=TRUE, robust=FALSE)NNS.SS(x, y, ...)
+
+
+
5.2 Code: two‑sample & multi‑group
+
ctrl <- rnorm(200, 0, 1)
+trt <- rnorm(180, 0.35, 1.2)
+NNS.ANOVA(control=ctrl, treatment=trt, means.only=FALSE, plot=FALSE)
+
## $Control
+## [1] 0.05568255
+##
+## $Treatment
+## [1] 0.2771257
+##
+## $Grand_Statistic
+## [1] 0.1605767
+##
+## $Control_CDF
+## [1] 0.5670595
+##
+## $Treatment_CDF
+## [1] 0.4385169
+##
+## $Certainty
+## [1] 0.6905098
+##
+## $Effect_Size_LB
+## 2.5%
+## -0.07055716
+##
+## $Effect_Size_UB
+## 97.5%
+## 0.5317766
+##
+## $Confidence_Level
+## [1] 0.95
+
A <- list(g1=rnorm(150,0.0,1.1), g2=rnorm(150,0.2,1.0), g3=rnorm(150,-0.1,0.9))
+NNS.ANOVA(control=A, means.only=TRUE, plot=FALSE)
+
## Certainty
+## 0.6876008
+
Math sketch. For each quantile/threshold \(t\), compare CDFs built from
+LPM.ratio(0, t, •) (possibly with one‑sided tails).
+Aggregate across \(t\) to a certainty
+score.
+
+
+
5.3 Stochastic Superiority
+
Stochastic superiority asks a different question than equality of
+means or equality of distributions. Rather than testing whether two
+samples came from the same population, or whether they share the same
+mean or median, stochastic superiority measures the probability that a
+random draw from one distribution exceeds a random draw from
+another.
+
For two random variables \(X\) and
+\(Y\), the stochastic superiority
+probability is:
+
\[
+P(X > Y)
+\]
+
and with ties accounted for, the tie-adjusted stochastic superiority
+measure is:
+
\[
+P^* = P(X > Y) + \frac{1}{2} P(X = Y)
+\]
+
A value of \(P^* = 0.5\) indicates
+no directional advantage, values above \(0.5\) favor \(X\), and values below \(0.5\) favor \(Y\).
+
This differs from stochastic dominance. Stochastic superiority is a
+pairwise exceedance probability, while stochastic dominance requires one
+distribution to be preferred to another over the entire shared
+support.
+
Below is an example comparing two distributions with unequal
+means.
+
set.seed(123)
+x = rnorm(1000, mean = 0, sd = 1)
+y = rnorm(1000, mean = 1, sd = 1)
+
+NNS.SS(x, y)
+
## $p_gt
+## [1] 0.233915
+##
+## $p_tie
+## [1] 0
+##
+## $p_star
+## [1] 0.233915
+
Since \(y\) was generated with a
+higher mean, the stochastic superiority probability for \(x\) relative to \(y\) should be less than \(0.5\), indicating that a draw from \(x\) is less likely to exceed a draw from
+\(y\).
+
We can also obtain confidence intervals for the tie-adjusted
+superiority probability using maximum entropy bootstrap replicates.
+
NNS.SS(x, y, confidence.interval = TRUE, reps = 999, ci = 0.95)[1:5]
+
+$p_gt
+[1] 0.233915
+
+$p_tie
+[1] 0
+
+$p_star
+[1] 0.233915
+
+$lower
+[1] 0.2105631
+
+$upper
+[1] 0.2537789
+
This provides an interpretable effect size for directional comparison
+between two distributions without requiring identical distributions or
+equal variances.
+
For discrete variables, ties may occur with positive probability, and
+the reported p_tie and p_star values reflect
+that adjustment explicitly.
+
set.seed(123)
+x = sample(1:5, 100, replace = TRUE)
+y = sample(1:5, 100, replace = TRUE)
+
+NNS.SS(x, y)
+
## $p_gt
+## [1] 0.3982
+##
+## $p_tie
+## [1] 0.1992
+##
+## $p_star
+## [1] 0.4978
+
+
+
+
6. Regression, Boosting, Stacking & Causality
+
+
6.1 Philosophy
+
NNS.reg learns partitioned
+relationships using partial‑moment weights — linear where appropriate,
+nonlinear where needed — avoiding fragile global parametric forms.
+
Headers.
+
+NNS.reg(x, y, order=NULL, smooth=TRUE, ncores=1, ...) →
+$Fitted.xy, $Point.est, …NNS.boost(IVs.train, DV.train, IVs.test, epochs, learner.trials, status, balance, type, folds)NNS.stack(IVs.train, DV.train, IVs.test, type, balance, ncores, folds)NNS.caus(x, y) (directional causality score via
+conditional dependence)
+
+
+
6.2 Code: classification via regression + ensembles
+
# Example 1: Nonlinear regression
+set.seed(123)
+x_train <- runif(1000, -2, 2)
+y_train <- sin(pi * x_train) + rnorm(1000, sd = 0.2)
+
+x_test <- seq(-2, 2, length.out = 100)
+
+NNS.reg(x = x_train, y = y_train, order = NULL, point.est = x_test)
+
+
## $R2
+## [1] 0.9276761
+##
+## $SE
+## [1] 0.2015258
+##
+## $Prediction.Accuracy
+## NULL
+##
+## $equation
+## NULL
+##
+## $x.star
+## NULL
+##
+## $derivative
+## Coefficient X.Lower.Range X.Upper.Range
+## <num> <num> <num>
+## 1: 3.0485215 -1.998138604 -1.934370540
+## 2: 3.5169373 -1.934370540 -1.804387149
+## 3: 1.8605016 -1.804387149 -1.692769075
+## 4: 0.6783073 -1.692769075 -1.590915710
+## 5: 0.4272848 -1.590915710 -1.465816449
+## 6: -0.5144026 -1.465816449 -1.376464546
+## 7: -1.9381128 -1.376464546 -1.229726997
+## 8: -3.0106084 -1.229726997 -1.110428636
+## 9: -2.5210796 -1.110428636 -0.976623793
+## 10: -3.7347021 -0.976623793 -0.870193992
+## 11: -2.0861598 -0.870193992 -0.754706576
+## 12: -2.1796417 -0.754706576 -0.636846031
+## 13: -0.9300308 -0.636846031 -0.533099369
+## 14: 1.0359249 -0.533099369 -0.417818767
+## 15: 0.9115004 -0.417818767 -0.323764665
+## 16: 2.3250859 -0.323764665 -0.184330858
+## 17: 3.0769180 -0.184330858 -0.132632209
+## 18: 3.3162510 -0.132632209 -0.080933560
+## 19: 3.5323950 -0.080933560 -0.004108338
+## 20: 2.1862481 -0.004108338 0.121863569
+## 21: 3.4805229 0.121863569 0.216987038
+## 22: 1.8001452 0.216987038 0.336388996
+## 23: 0.2295375 0.336388996 0.516729182
+## 24: -0.5625172 0.516729182 0.668479078
+## 25: -2.5532272 0.668479078 0.830264570
+## 26: -2.4765129 0.830264570 0.988320504
+## 27: -3.1248612 0.988320504 1.083380900
+## 28: -2.9622550 1.083380900 1.218812429
+## 29: -1.5047059 1.218812429 1.279773569
+## 30: -1.5723118 1.279773569 1.445979675
+## 31: 0.1804598 1.445979675 1.571628940
+## 32: 0.8726461 1.571628940 1.689536565
+## 33: 3.4198918 1.689536565 1.860999223
+## 34: 1.5206901 1.860999223 1.997618112
+## Coefficient X.Lower.Range X.Upper.Range
+##
+## $Point.est
+## [1] 0.01571202 0.10442684 0.23470933 0.37680781 0.51890629 0.65039140
+## [7] 0.72556318 0.80073496 0.85698998 0.88439634 0.91180269 0.93033285
+## [13] 0.94759688 0.96486091 0.95248720 0.93170326 0.87827479 0.79996720
+## [19] 0.72165961 0.64335202 0.52449573 0.40285499 0.28121424 0.17901835
+## [25] 0.07715655 -0.02470525 -0.15949123 -0.31038828 -0.45880078 -0.54309006
+## [31] -0.62737935 -0.71234468 -0.80041101 -0.88847734 -0.96331853 -1.00089554
+## [37] -1.03847254 -1.02090671 -0.97905116 -0.93719561 -0.89956805 -0.86273975
+## [43] -0.79660165 -0.70265879 -0.60871593 -0.51288396 -0.38856404 -0.25667590
+## [49] -0.11829230 0.02443074 0.13442845 0.22276171 0.31109496 0.42473203
+## [55] 0.56535922 0.69718932 0.76992246 0.84265560 0.90432326 0.91359751
+## [61] 0.92287175 0.93214599 0.94142023 0.92794242 0.90521445 0.88248648
+## [67] 0.85975852 0.76020581 0.65704511 0.55388442 0.45072372 0.35051056
+## [73] 0.25044944 0.15038831 0.04930377 -0.07695325 -0.20321027 -0.32495818
+## [79] -0.44464525 -0.56433232 -0.66432673 -0.72512293 -0.78817431 -0.85170206
+## [85] -0.91522981 -0.97875756 -0.99186193 -0.98457062 -0.97727932 -0.95314667
+## [91] -0.91788824 -0.88262981 -0.77697785 -0.63880041 -0.50062296 -0.36244551
+## [97] -0.25805231 -0.19661029 -0.13516826 -0.07526941
+##
+## $pred.int
+## NULL
+##
+## $regression.points
+## x y
+## <num> <num>
+## 1: -1.998138604 -0.01307124
+## 2: -1.934370540 0.18132707
+## 3: -1.804387149 0.63847051
+## 4: -1.692769075 0.84613612
+## 5: -1.590915710 0.91522399
+## 6: -1.465816449 0.96867700
+## 7: -1.376464546 0.92271416
+## 8: -1.229726997 0.63832023
+## 9: -1.110428636 0.27915958
+## 10: -0.976623793 -0.05817308
+## 11: -0.870193992 -0.45565668
+## 12: -0.754706576 -0.69658188
+## 13: -0.636846031 -0.95347564
+## 14: -0.533099369 -1.04996323
+## 15: -0.417818767 -0.93054118
+## 16: -0.323764665 -0.84481083
+## 17: -0.184330858 -0.52061525
+## 18: -0.132632209 -0.36154275
+## 19: -0.080933560 -0.19009705
+## 20: -0.004108338 0.08127998
+## 21: 0.121863569 0.35668582
+## 22: 0.216987038 0.68776523
+## 23: 0.336388996 0.90270609
+## 24: 0.516729182 0.94410093
+## 25: 0.668479078 0.85873901
+## 26: 0.830264570 0.44566388
+## 27: 0.988320504 0.05423632
+## 28: 1.083380900 -0.24281423
+## 29: 1.218812429 -0.64399694
+## 30: 1.279773569 -0.73572553
+## 31: 1.445979675 -0.99705336
+## 32: 1.571628940 -0.97437872
+## 33: 1.689536565 -0.87148709
+## 34: 1.860999223 -0.28510334
+## 35: 1.997618112 -0.07734835
+## x y
+##
+## $Fitted.xy
+## x y y.hat NNS.ID gradient residuals
+## <num> <num> <num> <char> <num> <num>
+## 1: -0.8496899 -0.5752368 -0.4984314 q121122 -2.0861598 0.076805376
+## 2: 1.1532205 -0.6617217 -0.4496971 q221122 -2.9622550 0.212024652
+## 3: -0.3640923 -0.7048691 -0.8815695 q122122 0.9115004 -0.176700402
+## 4: 1.5320696 -0.8447168 -0.9815176 q222121 0.1804598 -0.136800802
+## 5: 1.7618691 -0.9820881 -0.6241175 q222212 3.4198918 0.357970569
+## ---
+## 996: 1.3184955 -0.7988901 -0.7966085 q221222 -1.5723118 0.002281548
+## 997: 0.5684553 1.1554781 0.9150041 q212122 -0.5625172 -0.240473993
+## 998: -0.4340050 -0.7748325 -0.9473089 q122121 1.0359249 -0.172476359
+## 999: 0.8383194 0.7041960 0.4257159 q212222 -2.4765129 -0.278480031
+## 1000: -1.5647037 0.9467853 0.9264240 q111222 0.4272848 -0.020361366
+## standard.errors
+## <num>
+## 1: 0.1769692
+## 2: 0.1783713
+## 3: 0.1905081
+## 4: 0.2044300
+## 5: 0.2636784
+## ---
+## 996: 0.1971693
+## 997: 0.2137362
+## 998: 0.1831159
+## 999: 0.2108312
+## 1000: 0.2078031
+
# Simple train/test for boosting & stacking
+test.set = 141:150
+
+boost <- NNS.boost(IVs.train = iris[-test.set, 1:4],
+ DV.train = iris[-test.set, 5],
+ IVs.test = iris[test.set, 1:4],
+ epochs = 10, learner.trials = 10,
+ status = FALSE, balance = TRUE,
+ type = "CLASS", folds = 5)
+
+
+mean(boost$results == as.numeric(iris[test.set,5]))
+# [1] 1
+
+
+boost$feature.weights; boost$feature.frequency
+
+stacked <- NNS.stack(IVs.train = iris[-test.set, 1:4],
+ DV.train = iris[-test.set, 5],
+ IVs.test = iris[test.set, 1:4],
+ type = "CLASS", balance = TRUE,
+ ncores = 1, folds = 1)
+mean(stacked$stack == as.numeric(iris[test.set,5]))
+# [1] 1
+
+
+
6.3 Code: directional causality
+
NNS.caus(mtcars$hp, mtcars$mpg) # hp -> mpg
+
## Causation.x.given.y Causation.y.given.x C(x--->y)
+## 0.2607148 0.3863580 0.3933374
+
NNS.caus(mtcars$mpg, mtcars$hp) # hp -> mpg
+
## Causation.x.given.y Causation.y.given.x C(y--->x)
+## 0.3863580 0.2607148 0.3933374
+
Interpretation. Examine asymmetry in scores to infer
+direction. The method conditions partial‑moment dependence on candidate
+drivers.
+
+
+
+
7. Time Series & Forecasting
+
Headers.
+
+NNS.ARMANNS.ARMA.optimNNS.seasNNS.VAR
+
# Univariate nonlinear ARMA
+z <- as.numeric(scale(sin(1:480/8) + rnorm(480, sd=.35)))
+
+# Seasonality detection (prints a summary)
+seasonal_period <- NNS.seas(z, plot = FALSE)
+head(seasonal_period$all.periods)
+
## Period Coefficient.of.Variation Variable.Coefficient.of.Variation
+## 1 200 0.4267885 8.540159e+16
+## 2 96 0.4425880 8.540159e+16
+## 3 49 0.4615546 8.540159e+16
+## 4 198 0.4812956 8.540159e+16
+## 5 199 0.4885608 8.540159e+16
+## 6 146 0.4901054 8.540159e+16
+
# Validate seasonal periods
+NNS.ARMA.optim(z, h = 48, seasonal.factor = seasonal_period$periods, plot = TRUE, ncores = 1)
+
## [1] "CURRNET METHOD: lin"
+## [1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
+## [1] "NNS.ARMA(... method = 'lin' , seasonal.factor = c( 51 ) ...)"
+## [1] "CURRENT lin OBJECTIVE FUNCTION = 0.398327414917885"
+## [1] "BEST method = 'lin', seasonal.factor = c( 51 )"
+## [1] "BEST lin OBJECTIVE FUNCTION = 0.398327414917885"
+## [1] "CURRNET METHOD: nonlin"
+## [1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
+## [1] "NNS.ARMA(... method = 'nonlin' , seasonal.factor = c( 51 ) ...)"
+## [1] "CURRENT nonlin OBJECTIVE FUNCTION = 2.75408671013046"
+## [1] "BEST method = 'nonlin' PATH MEMBER = c( 51 )"
+## [1] "BEST nonlin OBJECTIVE FUNCTION = 2.75408671013046"
+## [1] "CURRNET METHOD: both"
+## [1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
+## [1] "NNS.ARMA(... method = 'both' , seasonal.factor = c( 51 ) ...)"
+## [1] "CURRENT both OBJECTIVE FUNCTION = 0.778172239627562"
+## [1] "BEST method = 'both' PATH MEMBER = c( 51 )"
+## [1] "BEST both OBJECTIVE FUNCTION = 0.778172239627562"
+
+
## $periods
+## [1] 51
+##
+## $weights
+## NULL
+##
+## $obj.fn
+## [1] 0.3983274
+##
+## $method
+## [1] "lin"
+##
+## $shrink
+## [1] FALSE
+##
+## $nns.regress
+## [1] FALSE
+##
+## $bias.shift
+## [1] 0.01738357
+##
+## $errors
+## [1] -0.4754897523 -0.4609730867 -0.3018876142 0.0439513384 -0.1128600832
+## [6] 0.9193835234 0.0160010547 -0.7516578805 -0.8195384972 -0.1274709629
+## [11] -0.0093477175 0.1480424491 0.0345888303 0.0009331215 -0.2819915138
+## [16] -0.3474821395 -1.2543849202 -0.5442948705 -0.0049610072 -0.4702036102
+## [21] 0.1846614137 1.6541950586 -0.2046795992 0.9691745476 1.1460606178
+## [26] -0.5141738440 -1.3562787956 0.3853973272 -0.3364881552 -0.5604890777
+## [31] -0.3175309175 -0.1677932189 -0.1511705981 0.4541183441 -0.1377055180
+## [36] 0.4279932502 1.3576081283 -0.0645315976 1.1430476887 0.1399600873
+## [41] 0.0874395694 -0.3703494531 0.3046994756 0.2057574931 -0.7602832912
+## [46] 0.6902933417 0.2238850985 -0.2775974238 -0.8250763050 0.5817787408
+## [51] -0.8733350647 0.2906911996 0.1863210948 -0.2484232855 0.1444232735
+## [56] 1.1655644133 0.0821221969 -0.2813315730 -0.7959981329 -0.3601165470
+## [61] -0.4617740020 -0.0593491905 0.0143389607 0.1016580238 0.0300275332
+## [66] -1.7237406556 -0.0930802461 -0.9348574200 -0.7189682901 0.0700766333
+## [71] -0.3547205444 0.2333233909 0.6840012123 -0.0779445509 0.8409902584
+## [76] 0.0130711684 0.8074727217 -1.1462424589 0.0926963526 -0.4674150054
+## [81] 0.1308248298 -0.6493713604 0.0713668583 0.4889233461 0.4293197750
+## [86] -0.4397639878 0.4261287370 0.7556075116 0.6016698079 -0.1086690282
+## [91] 0.6426872057 -0.4175612763 0.0250816728 0.9344147185 0.5444153587
+## [96] -0.8746369897
+##
+## $results
+## [1] -0.495145166 -0.629911085 -0.423647703 -1.217211533 -1.313660334
+## [6] -1.507558621 -1.512809568 -1.244102492 -0.765300445 -2.402307464
+## [11] -1.325990243 -0.928756118 -1.819067479 -0.855732188 -1.152690586
+## [16] -1.039006594 -0.562011496 -1.103503510 -0.685097085 -0.727417601
+## [21] -0.044018500 -0.030435409 0.002633325 -0.314902491 0.232587264
+## [26] 1.030889038 0.556722546 0.680351082 1.101193382 0.941245213
+## [31] 1.648626820 1.225992916 1.806473740 0.964372963 1.627354696
+## [36] 0.460925955 1.318674310 1.692295367 0.854538440 0.768654797
+## [41] 0.739228654 1.582319086 0.402156303 0.902802567 0.718513288
+## [46] 0.086635865 0.193748286 0.283357285
+##
+## $lower.pred.int
+## [1] -1.67077922 -1.80554514 -1.59928176 -2.39284559 -2.48929439 -2.68319268
+## [7] -2.68844363 -2.41973655 -1.94093450 -3.57794152 -2.50162430 -2.10439018
+## [13] -2.99470154 -2.03136624 -2.32832464 -2.21464065 -1.73764555 -2.27913757
+## [19] -1.86073114 -1.90305166 -1.21965256 -1.20606947 -1.17300073 -1.49053655
+## [25] -0.94304679 -0.14474502 -0.61891151 -0.49528298 -0.07444068 -0.23438884
+## [31] 0.47299276 0.05035886 0.63083968 -0.21126109 0.45172064 -0.71470810
+## [37] 0.14304025 0.51666131 -0.32109562 -0.40697926 -0.43640540 0.40668503
+## [43] -0.77347775 -0.27283149 -0.45712077 -1.08899819 -0.98188577 -0.89227677
+##
+## $upper.pred.int
+## [1] 0.68048889 0.54572297 0.75198635 -0.04157748 -0.13802628 -0.33192456
+## [7] -0.33717551 -0.06846843 0.41033361 -1.22667341 -0.15035619 0.24687794
+## [13] -0.64343342 0.31990187 0.02294347 0.13662746 0.61362256 0.07213055
+## [19] 0.49053697 0.44821646 1.13161556 1.14519865 1.17826738 0.86073157
+## [25] 1.40822132 2.20652310 1.73235660 1.85598514 2.27682744 2.11687927
+## [31] 2.82426088 2.40162697 2.98210780 2.14000702 2.80298875 1.63656001
+## [37] 2.49430837 2.86792942 2.03017250 1.94428885 1.91486271 2.75795314
+## [43] 1.57779036 2.07843662 1.89414735 1.26226992 1.36938234 1.45899134
+
Notes. NNS seasonality uses coefficient of variation
+instead of ACF/PACFs, and NNS ARMA blends multiple seasonal periods into
+the linear or nonlinear regression forecasts.
+
+
+
8. Simulation & Bootstrap & Risk‑Neutral Rescaling
+
+
8.1 Maximum entropy bootstrap (shape‑preserving)
+
Header.
+
+NNS.meboot(x, reps=999, rho=NULL, type="spearman", drift=TRUE, ...)
+
x_ts <- cumsum(rnorm(350, sd=.7))
+mb <- NNS.meboot(x_ts, reps=5, rho = 1)
+dim(mb["replicates", ]$replicates)
+
## [1] 350 5
+
+
+
8.2 Monte Carlo over the full correlation space
+
Header.
+
+NNS.MC(x, reps=30, lower_rho=-1, upper_rho=1, by=.01, exp=1, type="spearman", ...)
+
mc <- NNS.MC(x_ts, reps=5, lower_rho=-1, upper_rho=1, by=.5, exp=1)
+length(mc$ensemble); names(mc$replicates)
+
## [1] 350
+
## [1] "rho = 1" "rho = 0.5" "rho = 0" "rho = -0.5" "rho = -1"
+
head(mc$replicates$`rho = 0`)
+
## Replicate 1 Replicate 2 Replicate 3 Replicate 4 Replicate 5
+## [1,] 8.561720 11.097841 12.140974 3.478574 16.25845
+## [2,] 4.989649 9.142348 6.298598 2.573488 11.23749
+## [3,] 5.489892 11.635826 9.151404 4.146175 13.61840
+## [4,] 7.175210 13.194315 11.614209 5.906763 19.23707
+## [5,] 8.443500 12.157572 13.263425 4.369562 13.40513
+## [6,] 7.386515 10.979258 11.705842 2.410838 15.31133
+
+
+
+
9. Portfolio & Stochastic Dominance
+
Stochastic dominance orders uncertain prospects for broad classes of
+risk‑averse utilities; partial moments supply practical, nonparametric
+estimators.
+
Headers.
+
+NNS.FSD.uni(x, y)NNS.SSD.uni(x, y)NNS.TSD.uni(x, y)NNS.SD.cluster(R)NNS.SD.efficient.set(R)
+
RA <- rnorm(240, 0.005, 0.03)
+RB <- rnorm(240, 0.003, 0.02)
+RC <- rnorm(240, 0.006, 0.04)
+
+NNS.FSD.uni(RA, RB)
+
## [1] 0
+
+
## [1] 0
+
+
## [1] 0
+
Rmat <- cbind(A=RA, B=RB, C=RC)
+try(NNS.SD.cluster(Rmat, degree = 1))
+
## $Clusters
+## $Clusters$Cluster_1
+## [1] "C" "A" "B"
+
try(NNS.SD.efficient.set(Rmat, degree = 1))
+
## Checking 1 of 2Checking 2 of 2
+
## [1] "C" "A" "B"
+
+
+
Appendix A — Measure‑theoretic sketch (why partial moments are
+rigorous)
+
Let \((\Omega, \mathcal{F},
+\mathbb{P})\) be a probability space, \(X: \Omega\to\mathbb{R}\) measurable. For
+any fixed \(t\in\mathbb{R}\), the sets
+\(\{X\le t\}\) and \(\{X>t\}\) are in \(\mathcal{F}\) because they are preimages of
+Borel sets. The population partial moments are
+
\[
+\operatorname{LPM}(k,t,X) = \int_{-\infty}^{t} (t-x)^k\, dF_X(x),
+\qquad
+\operatorname{UPM}(k,t,X) = \int_{t}^{\infty} (x-t)^k\, dF_X(x).
+\]
+
The empirical versions correspond to replacing \(F_X\) with the empirical measure \(\mathbb{P}_n\) (or CDF \(\hat F_n\)):
+
\[
+\widehat{\operatorname{LPM}}_k(t;X) = \int_{(-\infty,t]} (t-x)^k\,
+d\mathbb{P}_n(x),
+\qquad
+\widehat{\operatorname{UPM}}_k(t;X) = \int_{(t,\infty)} (x-t)^k\,
+d\mathbb{P}_n(x).
+\]
+
Centering at \(t=\mu_X\) yields the
+variance decomposition identity in Section 1.
+
+
+
Appendix B — Quick Reference (Grouped by Topic)
+
+
+
1. Partial Moments & Ratios
+
+LPM(degree, target, variable) — lower partial moment of
+order degree at target.UPM(degree, target, variable) — upper partial moment of
+order degree at target.LPM.ratio(degree, target, variable);
+UPM.ratio(...) — normalized shares; degree=0
+gives CDF.LPM.VaR(p, degree, variable) — partial-moment quantile
+at probability p.Co.LPM(degree_lpm, x, y, target_x, target_y, degree_y)
+— co-lower partial moment between two variables.Co.UPM(degree_upm, x, y, target_x, target_y, degree_y)
+— co-upper partial moment between two variables.D.LPM(degree, target, variable) — divergent lower
+partial moment (away from target).D.UPM(degree, target, variable) — divergent upper
+partial moment (away from target).NNS.CDF(x, target = NULL, points = NULL, plot = TRUE/FALSE)
+— CDF from partial moments.NNS.moments(x) — mean/var/skew/kurtosis via partial
+moments.
+
+
+
2. Descriptive Statistics & Distributions
+
+NNS.mode(x, multi = FALSE) — nonparametric
+mode(s).PM.matrix(l_degree, u_degree, target, variable, pop_adj)
+— co-/divergent partial-moment matrices.NNS.gravity(x, w = NULL) — partial-moment weighted
+location (gravity center).
+
See NNS Vignette: Getting Started with NNS:
+Partial Moments
+
+
+
+
4. Normalization & Rescaling
+
+NNS.norm(x, linear=FALSE) — normalization retaining
+target moments.NNS.rescale(x, a, b, method=c("minmax","riskneutral"), T=NULL, type=c("Terminal","Discounted"))
+— risk-neutral or min–max rescaling.
+
See NNS Vignette: Getting Started
+with NNS: Normalization and Rescaling
+
+
+
5. Hypothesis Testing
+
+NNS.ANOVA(control, treatment, ...) — certainty of
+equality (distributions or means).NNS.SS(x, y, ...) — stochastic superiority between two
+variables.
+
See NNS Vignette: Getting Started with
+NNS: Comparing Distributions
+
+
+
6. Regression, Classification & Causality
+
+NNS.part(x, y, ...) — partition analysis for variable
+segmentation.NNS.reg(x, y, ...) — partition-based
+regression/classification ($Fitted.xy,
+$Point.est).NNS.boost(IVs, DV, ...),
+NNS.stack(IVs, DV, ...) — ensembles using
+NNS.reg base learners.NNS.caus(x, y) — directional causality score.
+
See NNS Vignette: Getting Started
+with NNS: Clustering and Regression
+
See NNS Vignette: Getting Started with NNS:
+Classification
+
+
+
7. Differentiation & Slope Measures
+
+dy.dx(x, y) — numerical derivative of y
+with respect to x via NNS.reg.dy.d_(x, Y, var) — partial derivative of multivariate
+Y w.r.t. var.NNS.diff(x, y) — derivative via secant
+projections.
+
+
+
8. Time Series & Forecasting
+
+NNS.ARMA(...), NNS.ARMA.optim(...) —
+nonlinear ARMA modeling.NNS.seas(...) — detect seasonality.NNS.VAR(...) — nonlinear VAR modeling.NNS.nowcast(x, h, ...) — near-term nonlinear
+forecast.
+
See NNS Vignette: Getting
+Started with NNS: Forecasting
+
+
+
+
10. Portfolio Analysis & Stochastic Dominance
+
+NNS.FSD.uni(x, y), NNS.SSD.uni(x, y),
+NNS.TSD.uni(x, y) — univariate stochastic dominance
+tests.NNS.SD.cluster(R), NNS.SD.efficient.set(R)
+— dominance-based portfolio sets.
+
For complete references, please see the Vignettes linked above and
+their specific referenced materials.
+
+
Getting Started with NNS: Partial
+Moments
+Fred Viole
+
+
+
+
+
Partial Moments
+
Why is it necessary to parse the variance with partial moments? The
+additional information generated from partial moments permits a level of
+analysis simply not possible with traditional summary statistics.
+
Below are some basic equivalences demonstrating partial moments role
+as the elements of variance.
+
+
Mean
+
library(NNS)
+set.seed(123) ; x = rnorm(100) ; y = rnorm(100)
+
+mean(x)
+
## [1] 0.09040591
+
UPM(1, 0, x) - LPM(1, 0, x)
+
## [1] 0.09040591
+
+
+
Variance
+
# Sample Variance (base R):
+var(x)
+
## [1] 0.8332328
+
# Sample Variance:
+(UPM(2, mean(x), x) + LPM(2, mean(x), x)) * (length(x) / (length(x) - 1))
+
## [1] 0.8332328
+
# Population Adjustment of Sample Variance (base R):
+var(x) * ((length(x) - 1) / length(x))
+
## [1] 0.8249005
+
# Population Variance:
+UPM(2, mean(x), x) + LPM(2, mean(x), x)
+
## [1] 0.8249005
+
# Variance is also the co-variance of itself:
+(Co.LPM(1, x, x, mean(x), mean(x)) + Co.UPM(1, x, x, mean(x), mean(x)) - D.LPM(1, 1, x, x, mean(x), mean(x)) - D.UPM(1, 1, x, x, mean(x), mean(x)))
+
## [1] 0.8249005
+
+
+
Standard Deviation
+
+
## [1] 0.9128159
+
((UPM(2, mean(x), x) + LPM(2, mean(x), x)) * (length(x) / (length(x) - 1))) ^ .5
+
## [1] 0.9128159
+
+
+
First 4 Moments
+
The first 4 moments are returned with the function
+NNS.moments. For sample statistics, set
+population = FALSE.
+
+
## $mean
+## [1] 0.09040591
+##
+## $variance
+## [1] 0.8249005
+##
+## $skewness
+## [1] 0.06049948
+##
+## $kurtosis
+## [1] -0.161053
+
NNS.moments(x, population = FALSE)
+
## $mean
+## [1] 0.09040591
+##
+## $variance
+## [1] 0.8332328
+##
+## $skewness
+## [1] 0.06235774
+##
+## $kurtosis
+## [1] -0.1069186
+
+
+
Statistical Mode of a Continuous Distribution
+
NNS.mode offers support for discrete valued
+distributions as well as recognizing multiple modes.
+
# Continuous
+NNS.mode(x)
+
## [1] -0.4132834
+
# Discrete and multiple modes
+NNS.mode(c(1, 2, 2, 3, 3, 4, 4, 5), discrete = TRUE, multi = TRUE)
+
## [1] 2 3 4
+
+
+
Covariance
+
+
## [1] -0.04372107
+
(Co.LPM(1, x, y, mean(x), mean(y)) + Co.UPM(1, x, y, mean(x), mean(y)) - D.LPM(1, 1, x, y, mean(x), mean(y)) - D.UPM(1, 1, x, y, mean(x), mean(y))) * (length(x) / (length(x) - 1))
+
## [1] -0.04372107
+
+
+
Covariance Elements and Covariance Matrix
+
The covariance matrix \((\Sigma)\)
+is equal to the sum of the co-partial moments matrices less the
+divergent partial moments matrices. \[ \Sigma
+= CLPM + CUPM - DLPM - DUPM \]
+
cov.mtx = PM.matrix(LPM_degree = 1, UPM_degree = 1, target = 'mean', variable = cbind(x, y), pop_adj = TRUE)
+cov.mtx
+
## $cupm
+## x y
+## x 0.4299250 0.1033601
+## y 0.1033601 0.5411626
+##
+## $dupm
+## x y
+## x 0.0000000 0.1469182
+## y 0.1560924 0.0000000
+##
+## $dlpm
+## x y
+## x 0.0000000 0.1560924
+## y 0.1469182 0.0000000
+##
+## $clpm
+## x y
+## x 0.4033078 0.1559295
+## y 0.1559295 0.3939005
+##
+## $cov.matrix
+## x y
+## x 0.83323283 -0.04372107
+## y -0.04372107 0.93506310
+
# Reassembled Covariance Matrix
+cov.mtx$clpm + cov.mtx$cupm - cov.mtx$dlpm - cov.mtx$dupm
+
## x y
+## x 0.83323283 -0.04372107
+## y -0.04372107 0.93506310
+
# Standard Covariance Matrix
+cov(cbind(x, y))
+
## x y
+## x 0.83323283 -0.04372107
+## y -0.04372107 0.93506310
+
+
+
Pearson Correlation
+
+
## [1] -0.04953215
+
cov.xy = (Co.LPM(1, x, y, mean(x), mean(y)) + Co.UPM(1, x, y, mean(x), mean(y)) - D.LPM(1, 1, x, y, mean(x), mean(y)) - D.UPM(1, 1, x, y, mean(x), mean(y))) * (length(x) / (length(x) - 1))
+sd.x = ((UPM(2, mean(x), x) + LPM(2, mean(x), x)) * (length(x) / (length(x) - 1))) ^ .5
+sd.y = ((UPM(2, mean(y), y) + LPM(2, mean(y) , y)) * (length(y) / (length(y) - 1))) ^ .5
+cov.xy / (sd.x * sd.y)
+
## [1] -0.04953215
+
+
+
CDFs (Discrete and Continuous)
+
P = ecdf(x)
+P(0) ; P(1)
+LPM(0, 0, x) ; LPM(0, 1, x)
+
+# Vectorized targets:
+LPM(0, c(0, 1), x)
+
+plot(ecdf(x))
+points(sort(x), LPM(0, sort(x), x), col = "red")
+legend("left", legend = c("ecdf", "LPM.CDF"), fill = c("black", "red"), border = NA, bty = "n")
+
+
# Joint CDF:
+Co.LPM(0, x, y, 0, 0)
+
+# Vectorized targets:
+Co.LPM(0, x, y, c(0, 1), c(0, 1))
+
+# Copula
+# Transform x and y so that they are uniform
+u_x = LPM.ratio(0, x, x)
+u_y = LPM.ratio(0, y, y)
+
+# Value of copula at c(.5, .5)
+Co.LPM(0, u_x, u_y, .5, .5)
+
+# Continuous CDF:
+NNS.CDF(x, 1)
+
+# CDF with target:
+NNS.CDF(x, 1, target = mean(x))
+
+
# Survival Function:
+NNS.CDF(x, 1, type = "survival")
+
+
+
+
Numerical Integration
+
Partial moments are asymptotic area approximations of \(f(x)\) akin to the familiar Trapezoidal and
+Simpson’s rules. More observations, more accuracy…
+
\[[UPM(1,0,f(x))-LPM(1,0,f(x))]\asymp\frac{[F(b)-F(a)]}{[b-a]}\]
+\[[UPM(1,0,f(x))-LPM(1,0,f(x))] *[b-a]
+\asymp[F(b)-F(a)]\]
+
x = seq(0, 1, .001) ; y = x ^ 2
+(UPM(1, 0, y) - LPM(1, 0, y)) * (1 - 0)
+
## [1] 0.3335
+
\[0.3333 * [1-0] = \int_{0}^{1} x^2
+dx\] For the total area, not just the definite integral, simply
+sum the partial moments and multiply by \([b -
+a]\): \[[UPM(1,0,f(x))+LPM(1,0,f(x))]
+*[b-a]\asymp\left\lvert{\int_{a}^{b} f(x)dx}\right\rvert\]
+
+
+
Bayes’ Theorem
+
For example, when ascertaining the probability of an increase in
+\(A\) given an increase in \(B\), the
+Co.UPM(degree_upm, x, y, target_x, target_y) target
+parameters are set to target_x = 0 and
+target_y = 0 and the
+UPM(degree, target, variable) target parameter is also set
+to target = 0.
+
\[P(A|B)=\frac{Co.UPM(0,A,B,0,0)}{UPM(0,0,B)}\]
+
+
+
References
+
If the user is so motivated, detailed arguments and proofs are
+provided within the following:
+
+
Getting Started with NNS: Correlation and
+Dependence
+Fred Viole
+
+
+
+library(NNS)
+library(data.table)
+require(knitr)
+require(rgl)
+
Correlation and Dependence
+
The limitations of linear correlation are well known. Often one uses
+correlation, when dependence is the intended measure for defining the
+relationship between variables. NNS dependence
+NNS.dep is a signal:noise measure robust
+to nonlinear signals.
+
Below are some examples comparing NNS correlation
+NNS.cor and
+NNS.dep with the standard Pearson’s
+correlation coefficient cor.
+
+
Linear Equivalence
+
Note the fact that all observations occupy the co-partial moment
+quadrants.
+
x = seq(0, 3, .01) ; y = 2 * x
+
+
+
## [1] 1
+
+
## $Correlation
+## [1] 1
+##
+## $Dependence
+## [1] 1
+
+
+
Nonlinear Relationship
+
Note the fact that all observations occupy the co-partial moment
+quadrants.
+
x = seq(0, 3, .01) ; y = x ^ 10
+
+
+
## [1] 0.6610183
+
+
## $Correlation
+## [1] 0.9595032
+##
+## $Dependence
+## [1] 0.9595032
+
+
+
Cyclic Relationship
+
Even the difficult inflection points, which span both the co- and
+divergent partial moment quadrants, are properly compensated for in
+NNS.dep.
+
x = seq(0, 12*pi, pi/100) ; y = sin(x)
+
+
+
## [1] -0.1297766
+
+
## $Correlation
+## [1] 0.202252
+##
+## $Dependence
+## [1] 0.8197963
+
+
+
Asymmetrical Analysis
+
The asymmetrical analysis is critical for further determining a
+causal path between variables which should be identifiable, i.e., it is
+asymmetrical in causes and effects.
+
The previous cyclic example visually highlights the asymmetry of
+dependence between the variables, which can be confirmed using
+NNS.dep(..., asym = TRUE).
+
+
## [1] -0.1297766
+
NNS.dep(x, y, asym = TRUE)
+
## $Correlation
+## [1] 0.202252
+##
+## $Dependence
+## [1] 0.8197963
+
+
## [1] -0.1297766
+
NNS.dep(y, x, asym = TRUE)
+
## $Correlation
+## [1] 0.07270847
+##
+## $Dependence
+## [1] 0.4086234
+
+
+
Dependence
+
Note the fact that all observations occupy only co- or divergent
+partial moment quadrants for a given subquadrant.
+
set.seed(123)
+df = data.frame(x = runif(10000, -1, 1), y = runif(10000, -1, 1))
+df = subset(df, (x ^ 2 + y ^ 2 <= 1 & x ^ 2 + y ^ 2 >= 0.95))
+
+
+
## $Correlation
+## [1] 0.05834412
+##
+## $Dependence
+## [1] 0.46764
+
+
+
p-values for NNS.dep()
+
p-values and confidence intervals can be obtained from sampling
+random permutations of \(y \rightarrow
+y_p\) and running NNS.dep(x,$y_p$)
+to compare against a null hypothesis of 0 correlation, or independence
+between \((x, y)\).
+
Simply set
+NNS.dep(..., p.value = TRUE, print.map = TRUE)
+to run 100 permutations and plot the results.
+
## p-values for [NNS.dep]
+set.seed(123)
+x = seq(-5, 5, .1); y = x^2 + rnorm(length(x))
+
+
NNS.dep(x, y, p.value = TRUE, print.map = TRUE)
+
+
## $Correlation
+## [1] 0.2957015
+##
+## $`Correlation p.value`
+## [1] 0.18
+##
+## $`Correlation 95% CIs`
+## 2.5% 97.5%
+## -0.1544429 0.4062421
+##
+## $Dependence
+## [1] 0.7932674
+##
+## $`Dependence p.value`
+## [1] 0
+##
+## $`Dependence 95% CIs`
+## 2.5% 97.5%
+## 0.5467152 0.6782456
+
+
Multivariate Dependence NNS.copula()
+
These partial moment insights permit us to extend the analysis to
+multivariate instances and deliver a dependence measure \((D)\) such that \(D \in [0,1]\). This level of analysis is
+simply impossible with Pearson or other rank based correlation methods,
+which are restricted to bivariate cases.
+
set.seed(123)
+x = rnorm(1000); y = rnorm(1000); z = rnorm(1000)
+NNS.copula(cbind(x, y, z), plot = TRUE, independence.overlay = TRUE)
+
## [1] 0.3278775
+
+
References
+
If the user is so motivated, detailed arguments and proofs are
+provided within the following:
+
+
Getting Started with NNS: Normalization and
+Rescaling
+Fred Viole
+
+
+
+library(NNS)
+library(data.table)
+require(knitr)
+require(rgl)
+
Overview
+
This vignette covers two related tools:
+
+NNS.norm() for cross‑variable normalization when
+comparing multiple series.NNS.rescale() for single‑vector rescaling with either
+min‑max or risk‑neutral targets.
+
Both functions perform deterministic affine transformations that
+preserve rank structure while modifying scale.
+
+
+
NNS.norm(): Normalize Multiple Variables
+
NNS.norm() rescales variables to a common magnitude
+while preserving distributional structure. The method can be
+linear (all variables forced to have the same mean) or
+nonlinear (using dependence weights to produce a more
+nuanced scaling). In the nonlinear case, the degree of association
+between variables influences the final normalized values.
+
+
Mathematical Structure
+
Let \(X\) be an \(n \times p\) matrix of variables.
+
+
Step 1: Compute Mean Vector
+
\[
+m_j = \text{mean}(X_{\cdot j})
+\]
+
If any \(m_j = 0\), it is replaced
+with \(10^{-10}\) to prevent division
+by zero.
+
+
+
+
Step 2: Construct Mean Ratio Matrix
+
\[
+RG_{ij} = \frac{m_i}{m_j}
+\]
+
In R this corresponds to:
+
+
+
+
+
Step 3: Dependence Weight Matrix
+
If linear = FALSE:
+
+- If number of variables \(p <
+10\): \[
+W = |\mathrm{cor}(X)|
+\]
+- Otherwise: \[
+W = |D| \quad \text{where } D = \text{NNS.dep}(X)\$Dependence
+\]
NNS.dep() returns a symmetric matrix of nonlinear
+dependence measures.
+
+
If linear = TRUE, the weighting effectively becomes:
+
\[
+W_{ij} = 1
+\]
+
+
+
+
Step 4: Scaling Factors
+
\[
+s_j = \frac{1}{p} \sum_{i=1}^{p} RG_{ij} W_{ij}
+\]
+
Each column is scaled:
+
\[
+X_{\cdot j}^{*} = s_j X_{\cdot j}
+\]
+
+
+
+
+
Linear Case Proof
+
If \(W_{ij} = 1\):
+
\[
+s_j = \frac{1}{p} \sum_{i=1}^{p} \frac{m_i}{m_j}
+= \frac{\bar{m}}{m_j}
+\]
+
Then:
+
\[
+\text{mean}(X_{\cdot j}^{*}) = s_j m_j = \bar{m}
+\]
+
All variables share the same mean.
+
+
+
+
Nonlinear Case Interpretation
+
\[
+\text{mean}(X_{\cdot j}^{*})
+=
+\frac{1}{p}
+\sum_{i=1}^{p}
+m_i W_{ij}
+\]
+
Thus, the normalized mean becomes a dependence‑weighted average of
+original means. Variables more strongly dependent with higher‑mean
+variables scale upward more.
+
+
+
+
Examples
+
+
Basic Multivariate Example
+
This holds for any distribution type and can be applied to vectors of
+different lengths.
+
set.seed(123)
+
+A <- rnorm(100, mean = 0, sd = 1)
+B <- rnorm(100, mean = 0, sd = 5)
+C <- rnorm(100, mean = 10, sd = 1)
+D <- rnorm(100, mean = 10, sd = 10)
+
+X <- data.frame(A, B, C, D)
+
+# Linear scaling
+lin_norm <- NNS.norm(X, linear = TRUE, chart.type = NULL)
+head(lin_norm)
+ A Normalized B Normalized C Normalized D Normalized
+[1,] -29.929719 31.889828 5.819152 1.4264014
+[2,] -12.291609 -11.531393 5.396317 1.2388239
+[3,] 83.235911 11.073887 4.643781 0.3078703
+[4,] 3.765188 15.601030 5.029380 -0.2630481
+[5,] 6.904039 42.717726 4.572611 2.8193657
+[6,] 91.585447 2.021274 4.543080 6.6681079
+
+# Verify means are equal
+apply(lin_norm, 2, function(x) c(mean = mean(x), sd = sd(x)))
+
+ A Normalized B Normalized C Normalized D Normalized
+mean 4.827727 4.827727 4.8277270 4.827727
+sd 48.744888 43.407590 0.4531172 5.203436
+
Now compare with nonlinear scaling:
+
nonlin_norm <- NNS.norm(X, linear = FALSE, chart.type = NULL)
+head(nonlin_norm)
+ A Normalized B Normalized C Normalized D Normalized
+[1,] -2.7834653 0.32807768 3.178568 0.7439872
+[2,] -1.1431202 -0.11863321 2.947605 0.6461499
+[3,] 7.7409438 0.11392645 2.536550 0.1605800
+[4,] 0.3501627 0.16050101 2.747174 -0.1372015
+[5,] 0.6420759 0.43947344 2.497676 1.4705341
+[6,] 8.5174510 0.02079456 2.481545 3.4779738
+
+apply(nonlin_norm, 2, function(x) c(mean = mean(x), sd = sd(x)))
+
+ A Normalized B Normalized C Normalized D Normalized
+mean 0.4489788 0.04966692 2.637026 2.518062
+sd 4.5332769 0.44657066 0.247504 2.714025
+
Note that the means differ and the standard deviations are smaller
+than in the linear case, reflecting the dependence structure.
+
+
Normalize list of unequal vector lengths
+
set.seed(123)
+vec1 <- rnorm(n = 10, mean = 0, sd = 1)
+vec2 <- rnorm(n = 5, mean = 5, sd = 5)
+vec3 <- rnorm(n = 8, mean = 10, sd = 10)
+
+vec_list <- list(vec1, vec2, vec3)
+
+NNS.norm(vec_list)
+
+$`x_1 Normalized`
+ [1] 13.074058 -3.004912 -11.745878 25.406891 -4.647966 -5.481229 6.225165 5.920719 6.113733 9.640242
+
+$`x_2 Normalized`
+[1] 2.875960212 0.008876158 1.230826150 5.855582361 10.779166523
+
+$`x_3 Normalized`
+[1] 4.0749062 2.2395840 0.4067264 0.7457562 15.6445780 5.1941416 2.3326665 2.5622994
+
+
+
+
+
Quantile Normalization Comparison
+
Quantile normalization forces distributions to be identical. This is
+literally the opposite intended effect of NNS.norm, which
+preserves individual distribution shapes while aligning ranges. The
+quantile normalized series become identical in distribution, while the
+NNS methods retain the original patterns.
+
+
+
+
+
Practical Applications
+
Normalization eliminates the need for multiple y‑axis charts and
+prevents their misuse. By placing variables on the same axes with shared
+ranges, we enable more relevant conditional probability analyses. This
+technique, combined with time normalization, is used in
+NNS.caus() to identify causal relationships between
+variables.
+
+
+
+
NNS.rescale(): Distribution Rescaling
+
NNS.rescale() performs one‑dimensional affine
+transformations.
+
Function signature:
+
NNS.rescale(x, a, b, method = "minmax", T = NULL, type = "Terminal")
+
+
+
1) Min-Max Scaling
+
If method = "minmax":
+
\[
+x^{*}
+=
+a
++
+(b - a)
+\frac{x - \min(x)}
+{\max(x) - \min(x)}
+\]
+
Properties:
+
+- Preserves order
+- Maps support to \([a,b]\)
+- Linear transformation
+
+
+
+
Example
+
raw_vals <- c(-2.5, 0.2, 1.1, 3.7, 5.0)
+
+scaled_minmax <- NNS.rescale(
+ x = raw_vals,
+ a = 5,
+ b = 10,
+ method = "minmax",
+ T = NULL,
+ type = "Terminal"
+)
+
+cbind(raw_vals, scaled_minmax)
+#> raw_vals scaled_minmax
+#> [1,] -2.5 5.000000
+#> [2,] 0.2 6.800000
+#> [3,] 1.1 7.400000
+#> [4,] 3.7 9.133333
+#> [5,] 5.0 10.000000
+range(scaled_minmax)
+#> [1] 5 10
+
+
+
+
+
2) Risk-Neutral Scaling
+
If method = "riskneutral":
+
Let:
+
+- \(S_0 = a\)
+- \(r = b\)
+- \(T\) = time horizon
+
+
+
Terminal Type
+
Target:
+
\[
+\mathbb{E}[S_T] = S_0 e^{rT}
+\]
+
Transformation form:
+
\[
+x^{*}
+=
+x
+\cdot
+\frac{S_0 e^{rT}}
+{\text{mean}(x)}
+\]
+
This enforces the required expectation.
+
+
+
+
Discounted Type
+
Target:
+
\[
+\mathbb{E}[e^{-rT} S_T] = S_0
+\]
+
Equivalent to:
+
\[
+\mathbb{E}[S_T] = S_0 e^{rT}
+\]
+
but the returned series is scaled so that its discounted mean equals
+\(S_0\). In practice, the function
+applies the same multiplicative factor as above, because:
+
\[
+\text{mean}(e^{-rT} x^{*}) = e^{-rT} \cdot \text{mean}(x^{*}) = e^{-rT}
+\cdot S_0 e^{rT} = S_0.
+\]
+
+
+
+
+
Risk-Neutral Example
+
set.seed(123)
+S0 <- 100
+r <- 0.05
+T <- 1
+
+# Simulate a price path
+prices <- S0 * exp(cumsum(rnorm(250, 0.0005, 0.02)))
+
+rn_terminal <- NNS.rescale(
+ x = prices,
+ a = S0,
+ b = r,
+ method = "riskneutral",
+ T = T,
+ type = "Terminal"
+)
+
+c(
+ mean_original = mean(prices),
+ mean_rescaled = mean(rn_terminal),
+ target = S0 * exp(r * T)
+)
+
+mean_original mean_rescaled target
+ 109.7019 105.1271 105.1271
+
+
+
+
Discounted Example
+
rn_discounted <- NNS.rescale(
+ x = prices,
+ a = S0,
+ b = r,
+ method = "riskneutral",
+ T = T,
+ type = "Discounted"
+)
+
+c(
+ mean_rescaled = mean(rn_discounted),
+ target_discounted_mean = S0
+)
+
+ mean_rescaled target_discounted_mean
+ 100 100
+
+
+
+
Conceptual Summary
+
+
NNS.norm()
+
+- Multivariate
+- Dependence‑aware scaling
+- Equalizes means only in linear mode
+- Preserves shape and order
+
+
+
+
NNS.rescale()
+
+- Univariate
+- Affine transformation
+- Either range‑targeted or expectation‑targeted
+- Preserves rank structure
+
+
Both functions maintain monotonicity and are therefore compatible
+with NNS copula and dependence modeling frameworks.
+
set.seed(123)
+
+x <- rnorm(1000, 5, 2)
+y <- rgamma(1000, 3, 1)
+
+# Combine variables
+X <- cbind(x, y)
+
+# NNS normalization
+X_norm_lin <- NNS.norm(X, linear = TRUE)
+X_norm_nonlin <- NNS.norm(X, linear = FALSE)
+
+# Standard min-max normalization
+minmax <- function(v) (v - min(v)) / (max(v) - min(v))
+X_minmax <- apply(X, 2, minmax)
+
+
+
+
+
References
+
If the user is so motivated, detailed arguments further examples are
+provided within the following:
+
+
Getting Started with NNS: Sampling and
+Simulation
+Fred Viole
+
+
+
+library(NNS)
+library(data.table)
+require(knitr)
+require(rgl)
NNS offers several novel sampling methods from any
+distribution, as well as simulating variables while maintaining their
+dependence.
+
+
Sampling
+
+
CDFs
+
Cumulative distribution functions (CDFs) represent the probability a
+variable \(X\) will take a value less
+than or equal to \(x\). \[F(x) = P(X \leq x)\]
+
+
Empirical CDF
+
The empirical CDF is a simple construct, provided in the base package
+of R. We can generate an empirical CDF with the ecdf
+function and create a function (P) to return the CDF of a
+given value of \(X\).
+
set.seed(123); x = rnorm(100)
+ecdf(x)
+
## Empirical CDF
+## Call: ecdf(x)
+## x[1:100] = -2.3092, -1.9666, -1.6867, ..., 2.169, 2.1873
+
+
## [1] 0.48
+
## [1] 0.83
+
+
+
Lower Partial Moment CDF
+(LPM.ratio)
+
The empirical CDF and Lower Partial Moment CDF
+(LPM.ratio) are identical when the degree
+term of the LPM.ratio is set to zero.
+
Degree 0 LPM: \[LPM(0,t,X)=\frac{1}{N}\sum_{n=1}^{N}[max(t-X_n),0]^0\]
+LPM.ratio is equivalent to the following form for any
+target \((t)\) and variable \(X\): \[LPM(0,t,X)=\frac{LPM(0,t,X)}{LPM(0,t,X)+UPM(0,t,X)}\]
+
Using the same targets from our ecdf example above (0,1)
+we can compare LPM.ratios.
+
LPM.ratio(degree = 0, target = 0, variable = x); LPM.ratio(degree = 0, target = 1, variable = x)
+
## [1] 0.48
+
## [1] 0.83
+
Calculating the probability for every target value in
+\(X\), we can plot both methods
+visualizing their identical results. ecdf function in black
+and LPM.ratio in red.
+
+
+
+
LPM.ratio degree > 0
+
By simply increasing the degree parameter to any
+positive real number, we can generate different CDFs of our initial
+distribution \(x\).
+
+
+
+
Generating PDFs with (LPM.VaR)
+
We can now generate distributions using the same insights and
+degree manipulation in the corresponding
+LPM.VaR function, a la value-at-risk,
+providing inverse CDF estimates.
+
The general form in the following plots is:
+
LPM.VaR(percentile = seq(0, 1, length.out = 100), degree = 0, x = x)
+
Any length percentile can be used to sample from the
+underlying distribution \(x\).
+
+
Viewing the first 10 samples from each of the degrees
+compared to our original \(X\).
+
degree.0.samples = LPM.VaR(percentile = seq(0, 1, length.out = 100), degree = 0, x = x)
+degree.0.25.samples = LPM.VaR(percentile = seq(0, 1, length.out = 100), degree = 0.25, x = x)
+degree.0.5.samples = LPM.VaR(percentile = seq(0, 1, length.out = 100), degree = 0.5, x = x)
+degree.1.samples = LPM.VaR(percentile = seq(0, 1, length.out = 100), degree = 1, x = x)
+degree.2.samples = LPM.VaR(percentile = seq(0, 1, length.out = 100), degree = 2, x = x)
+
+head(data.table::data.table(cbind("original x" = sort(x), degree.0.samples,
+ degree.0.25.samples,
+ degree.0.5.samples,
+ degree.1.samples,
+ degree.2.samples)), 10)
+
+ original x degree.0.samples degree.0.25.samples degree.0.5.samples
+ 1: -2.309169 -2.309169 -2.309097 -2.3090915
+ 2: -1.966617 -1.966617 -1.941190 -1.6935509
+ 3: -1.686693 -1.686693 -1.599486 -1.4541494
+ 4: -1.548753 -1.548753 -1.382553 -1.2462731
+ 5: -1.265396 -1.265396 -1.250823 -1.1453748
+ 6: -1.265061 -1.265061 -1.176436 -1.0745440
+ 7: -1.220718 -1.220718 -1.119655 -1.0252742
+ 8: -1.138137 -1.138137 -1.067793 -0.9868693
+ 9: -1.123109 -1.123109 -1.026429 -0.9322105
+ 10: -1.071791 -1.071791 -1.014276 -0.8710942
+ degree.1.samples degree.2.samples
+ 1: -2.3091021 -2.3091170
+ 2: -1.4744653 -1.1614908
+ 3: -1.2159961 -0.9709972
+ 4: -1.0823023 -0.8610192
+ 5: -0.9968028 -0.7810300
+ 6: -0.9290505 -0.7169770
+ 7: -0.8666886 -0.6631888
+ 8: -0.8090433 -0.6170691
+ 9: -0.7556644 -0.5765608
+ 10: -0.7069835 -0.5403318
+
+
+
+
Simulation
+
+
Bootstrapping (NNS.meboot)
+
NNS.meboot is based on the maximum
+entropy bootstrap, available in the R-package meboot. This
+procedure is specifically designed for time-series and avoids the IID
+assumption in traditional methods.
+
The ability to sample from specified correlations ensures the full
+spectrum of future paths is sampled from. Typical Monte Carlo samples
+are restricted to [-0.3, 0.3] correlations to the original data.
+
We will generate 1 replicate of \(X\) for each value of a sequence of \(\rho\) values (the \(ensemble\)), and then plot the results
+compared to our original \(X\) (black
+line). NNS.MC is a streamlined wrapper
+function for this functionality of
+NNS.meboot.
+
boots = NNS.MC(x, reps = 1, lower_rho = -1, upper_rho = 1, by = .5)$replicates
+reps = do.call(cbind, boots)
+
+
+matplot(reps, type = "l", col = rainbow(length(boots)))
+lines(x, type = "l", lwd = 3, ylim = c(min(reps), max(reps)))
+
+
Checking our replicate correlations:
+
sapply(boots, function(r) cor(r, x, method = "spearman"))
+
+ rho = 1 rho = 0.5 rho = 0 rho = -0.5 rho = -1
+ 0.99732373 0.51147915 0.01036904 -0.48720072 -0.98294629
+
More replicates and ensembles thereof can be generated for any number
+of \(\rho\) values.
+
+
target_drift Specification
+
We can also specify a target drift in our replicates with the
+target_drift parameter.
+
boots = NNS.MC(x, reps = 1, lower_rho = -1, upper_rho = 1, by = .5, target_drift = 0.05)$replicates
+reps = do.call(cbind, boots)
+
+plot(x, type = "l", lwd = 3, ylim = c(min(c(x, reps)), max(c(x, reps))))
+matplot(reps, type = "l", col = rainbow(length(boots)), add = TRUE)
+
+
Please see the full NNS.meboot and
+NNS.MC argument documentation.
+
+
+
+
Simulating a Multivariate Dependence Structure
+
Analogous to an empirical copula transformation, we can generate
+new data from the dependence structure of our
+original data via the following steps:
+
+- Determine the dependence structure:
+
+
This is accomplished using
+LPM.ratio(1, x, x) for continuous
+variables, and LPM.ratio(0, x, x) for
+discrete variables, which are the empirical CDFs of the marginal
+variables.
+
+- Generate or supply
new data:
+
+
new data does not have to be of the same distribution or
+dimension as the original data, nor does each dimension of
+new data have to share a distribution type.
+
+- Apply dependence structure to
+
new data:
+
+
We then utilize LPM.VaR to ascertain
+new data values corresponding to original data
+position mappings, and return a matrix of these transformed values with
+the same dimensions as new.data.
+
set.seed(123)
+x = rnorm(1000); y = rnorm(1000); z = rnorm(1000)
+
+# Add variable x to original data to avoid total independence (example only)
+original.data = cbind(x, y, z, x)
+
+# Determine dependence structure
+dep.structure = apply(original.data, 2, function(x) LPM.ratio(degree = 1, target = x, variable = x))
+
+# Generate new data with different mean, sd and length (or distribution type)
+new.data = sapply(1:ncol(original.data), function(x) rnorm(nrow(original.data)*2, mean = 10, sd = 20))
+
+# Apply dependence structure to new data
+new.dep.data = sapply(1:ncol(original.data), function(x) LPM.VaR(percentile = dep.structure[,x], degree = 1, x = new.data[,x]))
+
+
Compare Multivariate Dependence Structures
+
Similar dependence with radically different values, since we used
+\(N(10, 20)\) in place of our original
+\(N(0,1)\) observations.
+
NNS.copula(original.data)
+NNS.copula(new.dep.data)
+
+[1] 0.4743531
+[1] 0.4753264
+
head(original.data)
+head(new.dep.data)
+
+ x y z x
+[1,] -0.56047565 -0.99579872 -0.5116037 -0.56047565
+[2,] -0.23017749 -1.03995504 0.2369379 -0.23017749
+[3,] 1.55870831 -0.01798024 -0.5415892 1.55870831
+[4,] 0.07050839 -0.13217513 1.2192276 0.07050839
+[5,] 0.12928774 -2.54934277 0.1741359 0.12928774
+[6,] 1.71506499 1.04057346 -0.6152683 1.71506499
+ [,1] [,2] [,3] [,4]
+[1,] -2.028109 -10.498044 -0.2090467 -1.682949
+[2,] 4.608303 -11.390485 15.6213689 4.852534
+[3,] 39.478741 8.836581 -0.8508203 40.585505
+[4,] 10.683731 6.609255 36.0328589 10.877677
+[5,] 11.866922 -47.955235 14.3111350 12.064633
+[6,] 42.665726 29.639640 -2.4141874 43.797025
+
+
+
+
Alternative Using NNS.meboot
+
Alternatively, if we wish to keep the simulated values close to the
+original data, we can apply the NNS.meboot
+procedure to each of the variables.
+
We will generate 1 replicate (for brevity) of \(\rho = 0.95\) to our
+original.data, use their ensemble and note the
+multivariate dependence among our new.boot.dep.data.
+
# Apply bootstrap to each variable
+new.boot.dep.data = apply(original.data, 2, function(r) NNS.meboot(r, reps = 100, rho = .95))
+
+# Reformat into vectors
+boot.ensemble.vectors = lapply(new.boot.dep.data, function(z) unlist(z["ensemble",]))
+
+# Create matrix from vectors
+new.boot.dep.matrix = do.call(cbind, boot.ensemble.vectors)
+
Checking ensemble correlations with
+original.data:
+
for(i in 1:4) print(cor(new.boot.dep.matrix[,i], original.data[,i], method = "spearman"))
+
+[1] 0.9452863
+[1] 0.9499478
+[1] 0.945878
+[1] 0.9442845
+
+
Compare Multivariate Dependence Structures
+
Similar dependence with similar values.
+
NNS.copula(original.data)
+NNS.copula(new.boot.dep.matrix)
+
+[1] 0.4743531
+[1] 0.4517661
+
head(original.data)
+head(new.boot.dep.matrix)
+
+ x y z x
+[1,] -0.56047565 -0.99579872 -0.5116037 -0.56047565
+[2,] -0.23017749 -1.03995504 0.2369379 -0.23017749
+[3,] 1.55870831 -0.01798024 -0.5415892 1.55870831
+[4,] 0.07050839 -0.13217513 1.2192276 0.07050839
+[5,] 0.12928774 -2.54934277 0.1741359 0.12928774
+[6,] 1.71506499 1.04057346 -0.6152683 1.71506499
+ x y z x
+ensemble1 -0.4268047 -0.7794553 -0.6364458 -0.4642642
+ensemble2 -0.2965744 -1.0682197 0.3297265 -0.2531178
+ensemble3 1.3302149 0.3054734 -0.4014515 1.4914884
+ensemble4 0.2257378 0.3108846 1.0603892 0.1728540
+ensemble5 0.4716743 -3.3344967 -0.1917697 0.4309379
+ensemble6 1.3984978 1.1881374 -0.5295386 1.5326055
+
+
+
+
References
+
If the user is so motivated, detailed arguments and proofs are
+provided within the following:
+
+
Getting Started with NNS: Comparing
+Distributions
+Fred Viole
+
+
+
+library(NNS)
+library(data.table)
+require(knitr)
+require(rgl)
+
Comparing Distributions
+
NNS offers a multitude of ways to test
+if distributions came from the same population, or if they share the
+same mean or median. The underlying function for these tests is
+NNS.ANOVA().
+
The output from NNS.ANOVA() is a
+Certainty statistic, which compares CDFs of distributions
+from several shared quantiles and normalizes the similarity of these
+points to be within the interval \([0,1]\), with 1 representing identical
+distributions. For a complete analysis of Certainty to
+common p-values and the role of power, please see the References.
+
+
Test if Same Population
+
Below we run the analysis to whether automatic transmissions and
+manual transmissions have significantly different mpg
+distributions per the mtcars dataset.
+
The plot on the left shows the robust Certainty
+estimate, reflecting the distribution of Certainty
+estimates over 100 random permutations of both variables. The plot on
+the right illustrates the control and treatment variables, along with
+the grand mean among variables, and the confidence interval associated
+with the control mean.
+
mpg_auto_trans = mtcars[mtcars$am==1, "mpg"]
+mpg_man_trans = mtcars[mtcars$am==0, "mpg"]
+
+NNS.ANOVA(control = mpg_man_trans, treatment = mpg_auto_trans, robust = TRUE)
+
+
## $Control
+## [1] 17.14737
+##
+## $Treatment
+## [1] 24.39231
+##
+## $Grand_Statistic
+## [1] 20.09063
+##
+## $Control_CDF
+## [1] 0.8708501
+##
+## $Treatment_CDF
+## [1] 0.1294878
+##
+## $Certainty
+## [1] 0.02345583
+##
+## $`Effect_Size_LB.2.5%`
+## [1] 2.4708
+##
+## $`Effect_Size_UB.97.5%`
+## [1] 11.88554
+##
+## $Confidence_Level
+## [1] 0.95
+##
+## $`Robust Certainty Estimate`
+## [1] 0.01113453
+##
+## $`Lower 95% CI`
+## [1] 0
+##
+## $`Upper 95% CI`
+## [1] 0.1046841
+
The Certainty shows that these two distributions clearly
+do not come from the same population. This is verified with the
+Mann-Whitney-Wilcoxon test, which also does not assume a normality to
+the underlying data as a nonparametric test of identical
+distributions.
+
wilcox.test(mpg ~ am, data=mtcars)
+
##
+## Wilcoxon rank sum test with continuity correction
+##
+## data: mpg by am
+## W = 42, p-value = 0.001871
+## alternative hypothesis: true location shift is not equal to 0
+
+
+
Test if means are Equal
+
Here we provide the output from
+NNS.ANOVA() and t.test()
+functions on two Normal distribution samples, where we are pretty
+certain these two means are equal.
+
set.seed(123)
+x = rnorm(1000, mean = 0, sd = 1)
+y = rnorm(1000, mean = 0, sd = 2)
+
+NNS.ANOVA(control = x, treatment = y,
+ means.only = TRUE, robust = TRUE, plot = TRUE)
+
+
## $Control
+## [1] 0.01612787
+##
+## $Treatment
+## [1] 0.08493051
+##
+## $Grand_Statistic
+## [1] 0.05052919
+##
+## $Control_CDF
+## [1] 0.5218858
+##
+## $Treatment_CDF
+## [1] 0.4893919
+##
+## $Certainty
+## [1] 0.912545
+##
+## $`Effect_Size_LB.2.5%`
+## [1] -0.1215839
+##
+## $`Effect_Size_UB.97.5%`
+## [1] 0.2556542
+##
+## $Confidence_Level
+## [1] 0.95
+##
+## $`Robust Certainty Estimate`
+## [1] 0.9183685
+##
+## $`Lower 95% CI`
+## [1] 0.7311398
+##
+## $`Upper 95% CI`
+## [1] 0.9928339
+
+
##
+## Welch Two Sample t-test
+##
+## data: x and y
+## t = -0.96711, df = 1454.4, p-value = 0.3336
+## alternative hypothesis: true difference in means is not equal to 0
+## 95 percent confidence interval:
+## -0.20835512 0.07074984
+## sample estimates:
+## mean of x mean of y
+## 0.01612787 0.08493051
+
+
+
Test if means are Unequal
+
By altering the mean of the y variable, we can start to
+see the sensitivity of the results from the two methods, where both
+firmly reject the null hypothesis of identical means.
+
set.seed(123)
+x = rnorm(1000, mean = 0, sd = 1)
+y = rnorm(1000, mean = 1, sd = 1)
+
+NNS.ANOVA(control = x, treatment = y,
+ means.only = TRUE, robust = TRUE, plot = TRUE)
+
+
## $Control
+## [1] 0.01612787
+##
+## $Treatment
+## [1] 1.042465
+##
+## $Grand_Statistic
+## [1] 0.5292966
+##
+## $Control_CDF
+## [1] 0.7862176
+##
+## $Treatment_CDF
+## [1] 0.2197938
+##
+## $Certainty
+## [1] 0.1824463
+##
+## $`Effect_Size_LB.2.5%`
+## [1] 0.900412
+##
+## $`Effect_Size_UB.97.5%`
+## [1] 1.148409
+##
+## $Confidence_Level
+## [1] 0.95
+##
+## $`Robust Certainty Estimate`
+## [1] 0.1788691
+##
+## $`Lower 95% CI`
+## [1] 0.1484567
+##
+## $`Upper 95% CI`
+## [1] 0.2114944
+
+
##
+## Welch Two Sample t-test
+##
+## data: x and y
+## t = -22.933, df = 1997.4, p-value < 2.2e-16
+## alternative hypothesis: true difference in means is not equal to 0
+## 95 percent confidence interval:
+## -1.1141064 -0.9385684
+## sample estimates:
+## mean of x mean of y
+## 0.01612787 1.04246525
+
The effect size from NNS.ANOVA() is
+calculated from the confidence interval of the control mean and the
+specified y shift of 1 is within the provided lower and
+upper effect boundaries.
+
+
+
+
Stochastic Superiority
+
Stochastic superiority asks a different question than equality of
+means or equality of distributions. Rather than testing whether two
+samples came from the same population, or whether they share the same
+mean or median, stochastic superiority measures the probability that a
+random draw from one distribution exceeds a random draw from
+another.
+
For two random variables \(X\) and
+\(Y\), the stochastic superiority
+probability is:
+
\[
+P(X > Y)
+\]
+
and with ties accounted for, the tie-adjusted stochastic superiority
+measure is:
+
\[
+P^* = P(X > Y) + \frac{1}{2} P(X = Y)
+\]
+
A value of \(P^* = 0.5\) indicates
+no directional advantage, values above \(0.5\) favor \(X\), and values below \(0.5\) favor \(Y\).
+
This differs from stochastic dominance. Stochastic superiority is a
+pairwise exceedance probability, while stochastic dominance requires one
+distribution to be preferred to another over the entire shared
+support.
+
Below is an example using the same data generating process from the
+unequal means example.
+
set.seed(123)
+x = rnorm(1000, mean = 0, sd = 1)
+y = rnorm(1000, mean = 1, sd = 1)
+
+NNS.SS(x, y)
+
## $p_gt
+## [1] 0.233915
+##
+## $p_tie
+## [1] 0
+##
+## $p_star
+## [1] 0.233915
+
Since \(y\) was generated with a
+higher mean, the stochastic superiority probability for \(x\) relative to \(y\) should be less than \(0.5\), indicating that a draw from \(x\) is less likely to exceed a draw from
+\(y\).
+
We can also obtain confidence intervals for the tie-adjusted
+superiority probability using maximum entropy bootstrap replicates.
+
NNS.SS(x, y, confidence.interval = TRUE, reps = 999, ci = 0.95)[1:5]
+
+$p_gt
+[1] 0.233915
+
+$p_tie
+[1] 0
+
+$p_star
+[1] 0.233915
+
+$lower
+[1] 0.2105631
+
+$upper
+[1] 0.2537789
+
This provides an interpretable effect size for directional comparison
+between two distributions without requiring identical distributions or
+equal variances.
+
For discrete variables, ties may occur with positive probability, and
+the reported p_tie and p_star values reflect
+that adjustment explicitly.
+
set.seed(123)
+x = sample(1:5, 100, replace = TRUE)
+y = sample(1:5, 100, replace = TRUE)
+
+NNS.SS(x, y)
+
## $p_gt
+## [1] 0.3982
+##
+## $p_tie
+## [1] 0.1992
+##
+## $p_star
+## [1] 0.4978
+
+
Stochastic Dominance
+
Another method of comparing distributions involves a test for
+stochastic dominance. The first, second, and third degree stochastic
+dominance tests are available in NNS
+via:
+
+NNS.FSD()
NNS.SSD()
NNS.TSD()
+
set.seed(123)
+x = rnorm(1000, mean = 0, sd = 1)
+y = rnorm(1000, mean = 1, sd = 1)
+
+NNS.FSD(x, y)
+
+
## [1] "Y FSD X"
+
NNS.FSD() correctly identifies the
+shift in the y variable we specified when testing for
+unequal means.
+
+
Stochastic Dominant Efficient Sets
+
NNS also offers the ability to isolate
+a set of variables that do not have any dominated constituents with the
+NNS.SD.efficient.set() function.
+
x2, x4, x6, x8 all dominate their preceding
+distributions yet do not dominate one another, and are thus included in
+the first degree stochastic dominance efficient set.
+
set.seed(123)
+x1 = rnorm(1000)
+x2 = x1 + 1
+x3 = rnorm(1000)
+x4 = x3 + 1
+x5 = rnorm(1000)
+x6 = x5 + 1
+x7 = rnorm(1000)
+x8 = x7 + 1
+
+NNS.SD.efficient.set(cbind(x1, x2, x3, x4, x5, x6, x7, x8), degree = 1, status = FALSE)
+
## [1] "x4" "x2" "x8" "x6"
+
+
+
Stochastic Dominant Clusters
+
Further, we can assign clusters to non dominated constituents and
+represent the clustering in a dendrogram.
+
NNS.SD.cluster(cbind(x1, x2, x3, x4, x5, x6, x7, x8), degree = 1, dendrogram = TRUE)
+
+
## $Clusters
+## $Clusters$Cluster_1
+## [1] "x4" "x2" "x8" "x6"
+##
+## $Clusters$Cluster_2
+## [1] "x3" "x1" "x7" "x5"
+##
+##
+## $Dendrogram
+##
+## Call:
+## hclust(d = dist_matrix, method = "complete")
+##
+## Cluster method : complete
+## Number of objects: 8
+
+
+
References
+
If the user is so motivated, detailed arguments and proofs are
+provided within the following:
+
+
Getting Started with NNS: Clustering and
+Regression
+Fred Viole
+
+
+
+library(NNS)
+library(data.table)
+require(knitr)
+require(rgl)
+
Clustering and Regression
+
Below are some examples demonstrating unsupervised learning with NNS
+clustering and nonlinear regression using the resulting clusters. As
+always, for a more thorough description and definition, please view the
+References.
+
+
NNS Partitioning NNS.part()
+
NNS.part is both a partitional and
+hierarchical clustering method. NNS iteratively partitions
+the joint distribution into partial moment quadrants, and then assigns a
+quadrant identification (1:4) at each partition.
+
NNS.part returns a
+data.table of observations along with their final quadrant
+identification. It also returns the regression points, which are the
+quadrant means used in NNS.reg.
+
x = seq(-5, 5, .05); y = x ^ 3
+
+for(i in 1 : 4){NNS.part(x, y, order = i, Voronoi = TRUE, obs.req = 0)}
+
+
+
X-only Partitioning
+
NNS.part offers a partitioning based on
+\(x\) values only
+NNS.part(x, y, type = "XONLY", ...), using
+the entire bandwidth in its regression point derivation, and shares the
+same limit condition as partitioning via both \(x\) and \(y\) values.
+
for(i in 1 : 4){NNS.part(x, y, order = i, type = "XONLY", Voronoi = TRUE)}
+
+
Note the partition identifications are limited to 1’s and 2’s (left
+and right of the partition respectively), not the 4 values per the \(x\) and \(y\) partitioning.
+
## $order
+## [1] 4
+##
+## $dt
+## x y quadrant prior.quadrant
+## <num> <num> <char> <char>
+## 1: -5.00 -125.0000 q1111 q111
+## 2: -4.95 -121.2874 q1111 q111
+## 3: -4.90 -117.6490 q1111 q111
+## 4: -4.85 -114.0841 q1111 q111
+## 5: -4.80 -110.5920 q1111 q111
+## ---
+## 197: 4.80 110.5920 q2222 q222
+## 198: 4.85 114.0841 q2222 q222
+## 199: 4.90 117.6490 q2222 q222
+## 200: 4.95 121.2874 q2222 q222
+## 201: 5.00 125.0000 q2222 q222
+##
+## $regression.points
+## quadrant x y
+## <char> <num> <num>
+## 1: q111 -4.4523966 -89.31996002
+## 2: q112 -3.2250000 -31.51531806
+## 3: q121 -2.0023966 -7.46341667
+## 4: q122 -0.7590415 -0.51890098
+## 5: q211 0.3739355 0.08338409
+## 6: q212 1.3499632 2.26930682
+## 7: q221 2.6206250 16.42843100
+## 8: q222 4.1955267 75.78894504
+
+
+
+
Clusters Used in Regression
+
The right column of plots shows the corresponding regression (plus
+endpoints and central point) for the order of NNS
+partitioning.
+
for(i in 1 : 3){NNS.part(x, y, order = i, obs.req = 0, Voronoi = TRUE, type = "XONLY") ; NNS.reg(x, y, order = i, ncores = 1)}
+
+
+
+
NNS Regression NNS.reg()
+
NNS.reg can fit any \(f(x)\), for both uni- and multivariate
+cases. NNS.reg returns a self-evident list
+of values provided below.
+
+
Univariate:
+
NNS.reg(x, y, ncores = 1)
+
+
## $R2
+## [1] 0.9999858
+##
+## $SE
+## [1] 0.1822738
+##
+## $Prediction.Accuracy
+## NULL
+##
+## $equation
+## NULL
+##
+## $x.star
+## NULL
+##
+## $derivative
+## Coefficient X.Lower.Range X.Upper.Range
+## <num> <num> <num>
+## 1: 74.25250000 -5.0000000 -4.9750000
+## 2: 72.47650000 -4.9750000 -4.8500000
+## 3: 68.69350000 -4.8500000 -4.7250000
+## 4: 64.88656716 -4.7250000 -4.5854167
+## 5: 61.01480519 -4.5854167 -4.4250000
+## 6: 57.64628788 -4.4250000 -4.2875000
+## 7: 52.29438889 -4.2875000 -4.1000000
+## 8: 55.28971014 -4.1000000 -3.9562500
+## 9: 39.55816092 -3.9562500 -3.7750000
+## 10: 41.08694030 -3.7750000 -3.6354167
+## 11: 38.01863636 -3.6354167 -3.4750000
+## 12: 34.69626866 -3.4750000 -3.3354167
+## 13: 31.88168831 -3.3354167 -3.1750000
+## 14: 29.28265152 -3.1750000 -3.0375000
+## 15: 25.79438889 -3.0375000 -2.8500000
+## 16: 23.18416667 -2.8500000 -2.6875000
+## 17: 20.17544872 -2.6875000 -2.5250000
+## 18: 18.23350000 -2.5250000 -2.4000000
+## 19: 16.36150000 -2.4000000 -2.2750000
+## 20: 15.45634921 -2.2750000 -2.1437500
+## 21: 12.10506173 -2.1437500 -1.9750000
+## 22: 10.84291045 -1.9750000 -1.8354167
+## 23: 9.17837079 -1.8354167 -1.6500000
+## 24: 7.71713333 -1.6500000 -1.4937500
+## 25: 5.67487654 -1.4937500 -1.3250000
+## 26: 4.77015152 -1.3250000 -1.1875000
+## 27: 3.65525641 -1.1875000 -1.0250000
+## 28: 2.71828358 -1.0250000 -0.8854167
+## 29: 1.97577922 -0.8854167 -0.7250000
+## 30: 1.29696970 -0.7250000 -0.5875000
+## 31: 0.71536082 -0.5875000 -0.3854167
+## 32: 0.26031250 -0.3854167 -0.1854167
+## 33: 0.08077922 -0.1854167 -0.1052083
+## 34: 0.01168831 -0.1052083 -0.0250000
+## 35: 0.00625000 -0.0250000 0.0750000
+## 36: 0.05125000 0.0750000 0.1750000
+## 37: 0.17050000 0.1750000 0.3000000
+## 38: 0.40450000 0.3000000 0.4250000
+## 39: 0.68125000 0.4250000 0.5250000
+## 40: 0.99625000 0.5250000 0.6250000
+## 41: 1.30261905 0.6250000 0.7562500
+## 42: 2.23351852 0.7562500 0.9250000
+## 43: 2.85625000 0.9250000 1.0250000
+## 44: 3.47125000 1.0250000 1.1250000
+## 45: 4.21750000 1.1250000 1.2500000
+## 46: 5.19250000 1.2500000 1.3750000
+## 47: 6.18250000 1.3750000 1.5000000
+## 48: 7.35250000 1.5000000 1.6250000
+## 49: 7.76690476 1.6250000 1.7562500
+## 50: 10.84596774 1.7562500 1.9500000
+## 51: 10.93692308 1.9500000 2.1125000
+## 52: 14.30505155 2.1125000 2.3145833
+## 53: 17.95467391 2.3145833 2.5062500
+## 54: 21.46451613 2.5062500 2.7000000
+## 55: 20.50807692 2.7000000 2.8625000
+## 56: 26.01343750 2.8625000 3.0625000
+## 57: 32.71687737 3.0625000 3.2671585
+## 58: 34.19048114 3.2671585 3.5000000
+## 59: 33.57759494 3.5000000 3.6645833
+## 60: 46.95453488 3.6645833 3.8437500
+## 61: 42.67514286 3.8437500 4.0625000
+## 62: 57.09307692 4.0625000 4.2250000
+## 63: 55.24078947 4.2250000 4.3437500
+## 64: 59.68593153 4.3437500 4.5671585
+## 65: 66.33740696 4.5671585 4.8301031
+## 66: 72.01977335 4.8301031 5.0000000
+## Coefficient X.Lower.Range X.Upper.Range
+##
+## $Point.est
+## NULL
+##
+## $pred.int
+## NULL
+##
+## $regression.points
+## x y
+## <num> <num>
+## 1: -5.0000000 -1.250000e+02
+## 2: -4.9750000 -1.231437e+02
+## 3: -4.8500000 -1.140841e+02
+## 4: -4.7250000 -1.054974e+02
+## 5: -4.5854167 -9.644035e+01
+## 6: -4.4250000 -8.665256e+01
+## 7: -4.2875000 -7.872620e+01
+## 8: -4.1000000 -6.892100e+01
+## 9: -3.9562500 -6.097310e+01
+## 10: -3.7750000 -5.380319e+01
+## 11: -3.6354167 -4.806814e+01
+## 12: -3.4750000 -4.196931e+01
+## 13: -3.3354167 -3.712629e+01
+## 14: -3.1750000 -3.201194e+01
+## 15: -3.0375000 -2.798557e+01
+## 16: -2.8500000 -2.314913e+01
+## 17: -2.6875000 -1.938170e+01
+## 18: -2.5250000 -1.610319e+01
+## 19: -2.4000000 -1.382400e+01
+## 20: -2.2750000 -1.177881e+01
+## 21: -2.1437500 -9.750167e+00
+## 22: -1.9750000 -7.707437e+00
+## 23: -1.8354167 -6.193948e+00
+## 24: -1.6500000 -4.492125e+00
+## 25: -1.4937500 -3.286323e+00
+## 26: -1.3250000 -2.328687e+00
+## 27: -1.1875000 -1.672792e+00
+## 28: -1.0250000 -1.078812e+00
+## 29: -0.8854167 -6.993854e-01
+## 30: -0.7250000 -3.824375e-01
+## 31: -0.5875000 -2.041042e-01
+## 32: -0.3854167 -5.954167e-02
+## 33: -0.1854167 -7.479167e-03
+## 34: -0.1052083 -1.000000e-03
+## 35: -0.0250000 -6.250000e-05
+## 36: 0.0750000 5.625000e-04
+## 37: 0.1750000 5.687500e-03
+## 38: 0.3000000 2.700000e-02
+## 39: 0.4250000 7.756250e-02
+## 40: 0.5250000 1.456875e-01
+## 41: 0.6250000 2.453125e-01
+## 42: 0.7562500 4.162813e-01
+## 43: 0.9250000 7.931875e-01
+## 44: 1.0250000 1.078813e+00
+## 45: 1.1250000 1.425938e+00
+## 46: 1.2500000 1.953125e+00
+## 47: 1.3750000 2.602188e+00
+## 48: 1.5000000 3.375000e+00
+## 49: 1.6250000 4.294063e+00
+## 50: 1.7562500 5.313469e+00
+## 51: 1.9500000 7.414875e+00
+## 52: 2.1125000 9.192125e+00
+## 53: 2.3145833 1.208294e+01
+## 54: 2.5062500 1.552425e+01
+## 55: 2.7000000 1.968300e+01
+## 56: 2.8625000 2.301556e+01
+## 57: 3.0625000 2.821825e+01
+## 58: 3.2671585 3.491404e+01
+## 59: 3.5000000 4.287500e+01
+## 60: 3.6645833 4.840131e+01
+## 61: 3.8437500 5.681400e+01
+## 62: 4.0625000 6.614919e+01
+## 63: 4.2250000 7.542681e+01
+## 64: 4.3437500 8.198666e+01
+## 65: 4.5671585 9.532100e+01
+## 66: 4.8301031 1.127641e+02
+## 67: 5.0000000 1.250000e+02
+## x y
+##
+## $Fitted.xy
+## x y y.hat NNS.ID gradient residuals standard.errors
+## <num> <num> <num> <char> <num> <num> <num>
+## 1: -5.00 -125.0000 -125.0000 q1111111 74.25250 0.0000000 0.00000000
+## 2: -4.95 -121.2874 -121.3318 q1111112 72.47650 -0.0444000 0.07380015
+## 3: -4.90 -117.6490 -117.7080 q1111121 72.47650 -0.0589500 0.07380015
+## 4: -4.85 -114.0841 -114.0841 q1111121 68.69350 0.0000000 0.05069967
+## 5: -4.80 -110.5920 -110.6495 q1111122 68.69350 -0.0574500 0.05069967
+## ---
+## 197: 4.80 110.5920 110.7671 q2222221 66.33741 0.1751022 0.27620216
+## 198: 4.85 114.0841 114.1970 q2222222 72.01977 0.1129090 0.12572307
+## 199: 4.90 117.6490 117.7980 q2222222 72.01977 0.1490227 0.12572307
+## 200: 4.95 121.2874 121.3990 q2222222 72.01977 0.1116363 0.12572307
+## 201: 5.00 125.0000 125.0000 q2222222 72.01977 0.0000000 0.12572307
+
+
+
Multivariate:
+
Multivariate regressions return a plot of \(y\) and \(\hat{y}\), as well as the regression points
+($RPM) and partitions ($rhs.partitions) for
+each regressor.
+
f = function(x, y) x ^ 3 + 3 * y - y ^ 3 - 3 * x
+y = x ; z <- expand.grid(x, y)
+g = f(z[ , 1], z[ , 2])
+NNS.reg(z, g, order = "max", plot = FALSE, ncores = 1)
+
## $R2
+## [1] 1
+##
+## $rhs.partitions
+## Var1 Var2
+## <num> <num>
+## 1: -5.00 -5
+## 2: -4.95 -5
+## 3: -4.90 -5
+## 4: -4.85 -5
+## 5: -4.80 -5
+## ---
+## 40397: 4.80 5
+## 40398: 4.85 5
+## 40399: 4.90 5
+## 40400: 4.95 5
+## 40401: 5.00 5
+##
+## $RPM
+## Var1 Var2 y.hat
+## <num> <num> <num>
+## 1: -4.8 -4.80 -7.105427e-15
+## 2: -4.8 -2.55 -8.726063e+01
+## 3: -4.8 -2.50 -8.806700e+01
+## 4: -4.8 -2.45 -8.883587e+01
+## 5: -4.8 -2.40 -8.956800e+01
+## ---
+## 40397: -2.6 -2.80 3.776000e+00
+## 40398: -2.6 -2.75 2.770875e+00
+## 40399: -2.6 -2.70 1.807000e+00
+## 40400: -2.6 -2.65 8.836250e-01
+## 40401: -2.6 -2.60 1.776357e-15
+##
+## $Point.est
+## NULL
+##
+## $pred.int
+## NULL
+##
+## $Fitted.xy
+## Var1 Var2 y y.hat NNS.ID residuals
+## <num> <num> <num> <num> <char> <num>
+## 1: -5.00 -5 0.000000 0.000000 201.201 0
+## 2: -4.95 -5 3.562625 3.562625 402.201 0
+## 3: -4.90 -5 7.051000 7.051000 603.201 0
+## 4: -4.85 -5 10.465875 10.465875 804.201 0
+## 5: -4.80 -5 13.808000 13.808000 1005.201 0
+## ---
+## 40397: 4.80 5 -13.808000 -13.808000 39597.40401 0
+## 40398: 4.85 5 -10.465875 -10.465875 39798.40401 0
+## 40399: 4.90 5 -7.051000 -7.051000 39999.40401 0
+## 40400: 4.95 5 -3.562625 -3.562625 40200.40401 0
+## 40401: 5.00 5 0.000000 0.000000 40401.40401 0
+
+
+
+
NNS Dimension Reduction Regression
+
NNS.reg also provides a dimension
+reduction regression by including a parameter
+NNS.reg(x, y, dim.red.method = "cor", ...).
+Reducing all regressors to a single dimension using the returned
+equation
+NNS.reg(..., dim.red.method = "cor", ...)$equation.
+
NNS.reg(iris[ , 1 : 4], iris[ , 5], dim.red.method = "cor", location = "topleft", ncores = 1)$equation
+
+
## Variable Coefficient
+## <char> <num>
+## 1: Sepal.Length 0.7980781
+## 2: Sepal.Width -0.4402896
+## 3: Petal.Length 0.9354305
+## 4: Petal.Width 0.9381792
+## 5: DENOMINATOR 4.0000000
+
Thus, our model for this regression would be: \[Species = \frac{0.798*Sepal.Length
+-0.44*Sepal.Width +0.935*Petal.Length +0.938*Petal.Width}{4}
+\]
+
+
Threshold
+
NNS.reg(x, y, dim.red.method = "cor", threshold = ...)
+offers a method of reducing regressors further by controlling the
+absolute value of required correlation.
+
NNS.reg(iris[ , 1 : 4], iris[ , 5], dim.red.method = "cor", threshold = .75, location = "topleft", ncores = 1)$equation
+
+
## Variable Coefficient
+## <char> <num>
+## 1: Sepal.Length 0.7980781
+## 2: Sepal.Width 0.0000000
+## 3: Petal.Length 0.9354305
+## 4: Petal.Width 0.9381792
+## 5: DENOMINATOR 3.0000000
+
Thus, our model for this further reduced dimension regression would
+be: \[Species = \frac{\: 0.798*Sepal.Length +
+0*Sepal.Width +0.935*Petal.Length +0.938*Petal.Width}{3} \]
+
and the point.est = (...) operates in the same manner as
+the full regression above, again called with
+NNS.reg(...)$Point.est.
+
NNS.reg(iris[ , 1 : 4], iris[ , 5], dim.red.method = "cor", threshold = .75, point.est = iris[1 : 10, 1 : 4], location = "topleft", ncores = 1)$Point.est
+
+
## [1] 1 1 1 1 1 1 1 1 1 1
+
+
+
+
Classification
+
For a classification problem, we simply set
+NNS.reg(x, y, type = "CLASS", ...).
+
NOTE: Base category of response variable should be 1, not 0
+for classification problems.
+
NNS.reg(iris[ , 1 : 4], iris[ , 5], type = "CLASS", point.est = iris[1 : 10, 1 : 4], location = "topleft", ncores = 1)$Point.est
+
+
## [1] 1 1 1 1 1 1 1 1 1 1
+
+
Cross-Validation NNS.stack()
+
The NNS.stack routine cross-validates
+for a given objective function the n.best parameter in the
+multivariate NNS.reg function as well as
+the threshold parameter in the dimension reduction
+NNS.reg version.
+NNS.stack can be used for
+classification:
+
NNS.stack(..., type = "CLASS", ...)
+
or continuous dependent variables:
+
NNS.stack(..., type = NULL, ...).
+
Any objective function obj.fn can be called using
+expression() with the terms predicted and
+actual, even from external packages such as
+Metrics.
+
NNS.stack(..., obj.fn = expression(Metrics::mape(actual, predicted)), objective = "min").
+
NNS.stack(IVs.train = iris[ , 1 : 4],
+ DV.train = iris[ , 5],
+ IVs.test = iris[1 : 10, 1 : 4],
+ dim.red.method = "cor",
+ obj.fn = expression( mean(round(predicted) == actual) ),
+ objective = "max", type = "CLASS",
+ folds = 1, ncores = 1)
+
Folds Remaining = 0
+Current NNS.reg(... , threshold = 0.9350 ) | eval(obj.fn) = 1.000000 | MAX Iterations Remaining = 2
+Current NNS.reg(... , threshold = 0.7950 ) | eval(obj.fn) = 0.973684 | MAX Iterations Remaining = 1
+Current NNS.reg(... , threshold = 0.4400 ) | eval(obj.fn) = 0.894737 | MAX Iterations Remaining = 0
+Current NNS.reg(. , n.best = 1 ) | eval(obj.fn) = 0.868421 | MAX Iterations Remaining = 12
+Current NNS.reg(. , n.best = 2 ) | eval(obj.fn) = 0.736842 | MAX Iterations Remaining = 11
+Current NNS.reg(. , n.best = 3 ) | eval(obj.fn) = 0.763158 | MAX Iterations Remaining = 10
+Current NNS.reg(. , n.best = 4 ) | eval(obj.fn) = 0.736842 | MAX Iterations Remaining = 9
+$OBJfn.reg
+[1] 0.9733333
+
+$NNS.reg.n.best
+[1] 1
+
+$probability.threshold
+[1] 0.495
+
+$OBJfn.dim.red
+[1] 0.9666667
+
+$NNS.dim.red.threshold
+[1] 0.935
+
+$reg
+ [1] 1 1 1 1 1 1 1 1 1 1
+
+$reg.pred.int
+NULL
+
+$dim.red
+ [1] 1 1 1 1 1 1 1 1 1 1
+
+$dim.red.pred.int
+NULL
+
+$stack
+ [1] 1 1 1 1 1 1 1 1 1 1
+
+$pred.int
+NULL
+
+
Increasing Dimensions
+
Given multicollinearity is not an issue for nonparametric regressions
+as it is for OLS, in the case of an ill-fit univariate model a better
+option may be to increase the dimensionality of regressors with a copy
+of itself and cross-validate the number of clusters n.best
+via:
+
NNS.stack(IVs.train = cbind(x, x), DV.train = y, method = 1, ...).
+
set.seed(123)
+x = rnorm(100); y = rnorm(100)
+
+nns.params = NNS.stack(IVs.train = cbind(x, x),
+ DV.train = y,
+ method = 1, ncores = 1)
+
NNS.reg(cbind(x, x), y,
+ n.best = nns.params$NNS.reg.n.best,
+ point.est = cbind(x, x),
+ residual.plot = TRUE,
+ ncores = 1, confidence.interval = .95)
+
+
+
Smoothing Option
+
Smoothness is not required for curve fitting, but the
+NNS.reg function offers an optional smoothed fit. This
+feature applies a smoothing spline to regression points generated
+internally using the partitioning method described earlier.
+
NNS.reg(x, y, smooth = TRUE)
+
+
+
Imputation
+
Imputation in NNS is a direct application of nearest
+neighbor regression. When values of \(y\) are missing, we use the observed \((X,y)\) pairs as the training set and the
+predictors of the missing rows as point.est.
+
A key insight is that even in univariate regressions,
+NNS.reg benefits from the increasing dimensions trick: by
+duplicating the predictor into a multivariate form,
+e.g. cbind(x, x), the distance function underlying
+NNS.reg operates in a 2-D space. This sharpened distance
+metric allows a more robust donor selection, effectively turning
+univariate imputation into a special case of multivariate nearest
+neighbor regression.
+
For multivariate predictors, the same form applies directly — supply
+the full set of observed predictors in \(x\), the observed responses in \(y\), and the incomplete rows in
+point.est. With order = "max", n.best = 1, the
+imputation is always 1-NN donor-based: each missing \(y\) is filled in by the response of its
+closest donor under the NNS hybrid distance. This ensures
+imputations remain strictly within the support of the observed data.
+
Categorical data is handled analogously, only
+requiring NNS.reg(..., type = "CLASS") in the
+procedure.
+
+
Univariate Imputation
+
set.seed(123)
+
+# Univariate predictor with nonlinear signal
+n <- 400
+x <- sort(runif(n, -3, 3))
+y <- sin(x) + 0.2 * x^2 + rnorm(n, 0, 0.25)
+
+# Induce ~25% MCAR missingness in y
+miss <- rbinom(n, 1, 0.25) == 1
+y_mis <- y
+y_mis[miss] <- NA
+
+# ---- Increasing dimensions trick ----
+# Duplicate x so the distance operates in a 2D space: cbind(x, x).
+# This sharpens nearest-neighbor selection even in a nominally univariate setting.
+x2_train <- cbind(x[!miss], x[!miss])
+x2_miss <- cbind(x[miss], x[miss])
+
+# 1-NN donor imputation with NNS.reg
+y_hat_uni <- NNS::NNS.reg(
+ x = x2_train, # predictors (duplicated x)
+ y = y[!miss], # observed responses
+ point.est = x2_miss, # rows to impute
+ order = "max", # dependence-maximizing order
+ n.best = 1, # 1-NN donor
+ plot = FALSE
+)$Point.est
+
+# Fill back
+y_completed_uni <- y_mis
+y_completed_uni[miss] <- y_hat_uni
+
+# Plot observed vs imputed (NNS 1-NN)
+plot(x, y, pch = 1, col = "steelblue", cex = 1.5, lwd = 2,
+ xlab = "x", ylab = "y", main = "NNS 1-NN Imputation")
+points(x[miss], y_hat_uni, col = "red", pch = 15, cex = 1.3)
+
+legend("topleft",
+ legend = c("Observed", "Imputed (NNS 1-NN)"),
+ col = c("steelblue", "red"),
+ pch = c(1, 15),
+ pt.lwd = c(2, NA),
+ bty = "n")
+
+
+
+
+
Multivariate Imputation
+
set.seed(123)
+
+# Multivariate predictors with nonlinear & interaction structure
+n <- 800
+X <- cbind(
+ x1 = rnorm(n),
+ x2 = runif(n, -2, 2),
+ x3 = rnorm(n, 0, 1)
+)
+
+f <- function(x1, x2, x3) 1.1*x1 - 0.8*x2 + 0.5*x3 + 0.6*x1*x2 - 0.4*x2*x3 + 0.3*sin(1.3*x1)
+y <- f(X[,1], X[,2], X[,3]) + rnorm(n, 0, 0.4)
+
+# Induce ~30% MCAR missingness in y
+miss <- rbinom(n, 1, 0.30) == 1
+y_mis <- y
+y_mis[miss] <- NA
+
+# Training (observed) vs rows to impute
+X_obs <- X[!miss, , drop = FALSE]
+y_obs <- y[!miss]
+X_mis <- X[ miss, , drop = FALSE]
+
+# 1-NN donor imputation with NNS.reg
+y_hat_mv <- NNS::NNS.reg(
+ x = X_obs, # all observed predictors
+ y = y_obs, # observed responses
+ point.est = X_mis, # rows to impute
+ order = "max", # dependence-maximizing order
+ n.best = 1, # 1-NN donor
+ plot = FALSE
+)$Point.est
+
+# Completed vector
+y_completed_mv <- y_mis
+y_completed_mv[miss] <- y_hat_mv
+
+# Plot observed vs imputed (multivariate, NNS 1-NN)
+plot(seq_along(y), y,
+ pch = 1, col = "steelblue", cex = 1.5, lwd = 2,
+ xlab = "Observation index", ylab = "y",
+ main = "NNS 1-NN Multivariate Imputation")
+
+# Overlay imputed values
+points(which(miss), y_hat_mv, pch = 15, col = "red", cex = 1.2)
+
+# Legend
+legend("topleft",
+ legend = c("Observed", "Imputed (NNS 1-NN)"),
+ col = c("steelblue", "red"),
+ pch = c(1, 15),
+ pt.lwd = c(2, NA),
+ bty = "n")
+
+
+
+
+
A Note on Uncertainty Propagation
+
A common concern with local imputation methods is whether imputation
+uncertainty propagates correctly into downstream inference.
+NNS addresses this through bootstrap multiple imputation:
+resampling complete cases across m iterations generates
+between-imputation variance that flows through standard Rubin’s rules
+pooling identically to any classical procedure.
+
Empirically, NNS bootstrap MI outperforms MICE with
+predictive mean matching on nonlinear data — producing a pooled estimate
+closer to the true parameter with a smaller pooled SE. The advantage
+comes not from compressing uncertainty but from a more accurate
+imputation model, which reduces between-imputation variance driven by
+model error rather than genuine data uncertainty.
+
See NNS
+Multiple Imputation vs MICE for the full reproducible
+comparison.
+
+
+
References
+
If the user is so motivated, detailed arguments further examples are
+provided within the following:
+
+
Getting Started with NNS:
+Classification
+Fred Viole
+
+
+
+library(NNS)
+library(data.table)
+require(knitr)
+require(rgl)
+
Classification
+
NNS.reg is a very robust regression
+technique capable of nonlinear regressions of continuous variables and
+classification tasks in machine learning problems.
+
We have extended the NNS.reg
+applications per the use of an ensemble method of classification in
+NNS.boost. In short,
+NNS.reg is the base learner instead of
+trees.
+
One major advantage NNS.boost has over tree
+based methods is the ability to seamlessly extrapolate beyond the
+current range of observations.
+
+
Splits vs. Partitions
+
Popular boosting algorithms take a series of weak learning decision
+tree models, and aggregate their outputs. NNS is also a
+decision tree of sorts, by partitioning each regressor with respect to
+the dependent variable. We can directly control the number of “splits”
+with the NNS.reg(..., order = , ...)
+parameter.
+
+
NNS Partitions
+
We can see how NNS partitions each regressor by calling
+the $rhs.partitions output. You will notice that each
+partition is not an equal interval, nor of equal length, which
+differentiates NNS from other bandwidth or tree-based
+techniques.
+
Higher dependence between a regressor and the dependent variable will
+allow for a larger number of partitions. This is determined internally
+with the NNS.dep measure.
+
NNS.reg(iris[,1:4], iris[,5], residual.plot = FALSE, ncores = 1)$rhs.partitions
+
## V1 V2 V3 V4
+## <num> <num> <num> <num>
+## 1: 4.300000 2.000000 1.000000 0.100000
+## 2: 4.381250 2.645276 1.050000 0.200000
+## 3: 4.577396 2.980556 1.200000 0.300000
+## 4: 4.700000 3.181155 1.300000 0.400000
+## 5: 4.800000 3.552439 1.400000 0.500000
+## 6: 4.900000 4.400000 1.500000 0.600000
+## 7: 5.000000 NA 1.600000 1.000000
+## 8: 5.100000 NA 1.700000 1.100000
+## 9: 5.205000 NA 1.900000 1.200000
+## 10: 5.400000 NA 3.416305 1.300000
+## 11: 5.500000 NA 3.834865 1.400000
+## 12: 5.600000 NA 4.000000 1.500000
+## 13: 5.700000 NA 4.184722 1.600000
+## 14: 5.800000 NA 4.400000 1.700000
+## 15: 5.900000 NA 4.500000 1.800000
+## 16: 6.000000 NA 4.670803 1.900000
+## 17: 6.100000 NA 4.863889 2.000000
+## 18: 6.200000 NA 5.000000 2.117708
+## 19: 6.300000 NA 5.100000 2.300000
+## 20: 6.400000 NA 5.200000 2.435206
+## 21: 6.500000 NA 5.337500 2.500000
+## 22: 6.600000 NA 5.500000 NA
+## 23: 6.700000 NA 5.617708 NA
+## 24: 6.800000 NA 5.849554 NA
+## 25: 6.900000 NA 6.336875 NA
+## 26: 7.050000 NA 6.900000 NA
+## 27: 7.224375 NA NA NA
+## 28: 7.687079 NA NA NA
+## 29: 7.900000 NA NA NA
+## V1 V2 V3 V4
+
+
+
+
NNS.boost()
+
Through resampling of the training set and letting each iterated set
+of data speak for themselves (while paying extra attention to the
+residuals throughout), we can test various regressor combinations in
+these dynamic decision trees…only keeping those combinations that add
+predictive value. From there we simply aggregate the predictions.
+
NNS.boost will automatically search for
+an accuracy threshold from the training set, reporting
+iterations remaining and level obtained in the console. A plot of the
+frequency of the learning accuracy on the training set is also
+provided.
+
Once a threshold is obtained,
+NNS.boost will test various feature
+combinations against different splits of the training set and report
+back the frequency of each regressor used in the final estimate.
+
Let’s have a look and see how it works. We use 140 random
+iris observations as our training set with the 10 holdout
+observations as our test set. For brevity, we set
+epochs = 10, learner.trials = 10, folds = 1.
+
NOTE: Base category of response variable should be 1, not 0
+for classification problems when using
+NNS.boost(..., type = "CLASS").
+
test.set = 141:150
+
+a = NNS.boost(IVs.train = iris[-test.set, 1:4],
+ DV.train = iris[-test.set, 5],
+ IVs.test = iris[test.set, 1:4],
+ epochs = 10, learner.trials = 10,
+ status = FALSE, balance = TRUE,
+ type = "CLASS", folds = 5)
+
+a
+$results
+ [1] 3 3 3 3 3 3 3 3 3 3
+
+$pred.int
+NULL
+
+$feature.weights
+ Petal.Width Petal.Length Sepal.Length
+ 0.4285714 0.4285714 0.1428571
+
+$feature.frequency
+ Petal.Width Petal.Length Sepal.Length
+ 3 3 1
+
+mean( a$results == as.numeric(iris[test.set, 5]) )
+[1] 1
+
A perfect classification, using the features weighted per the output
+above.
+
+
Cross-Validation Classification Using NNS.stack()
+
The NNS.stack() routine cross-validates
+for a given objective function the n.best parameter in the
+multivariate NNS.reg function as well as
+the threshold parameter in the dimension reduction
+NNS.reg version.
+NNS.stack can be used for classification
+via
+NNS.stack(..., type = "CLASS", ...).
+
For brevity, we set folds = 1.
+
NOTE: Base category of response variable should be 1, not 0
+for classification problems when using
+NNS.stack(..., type = "CLASS").
+
b = NNS.stack(IVs.train = iris[-test.set, 1:4],
+ DV.train = iris[-test.set, 5],
+ IVs.test = iris[test.set, 1:4],
+ type = "CLASS", balance = TRUE,
+ ncores = 1, folds = 5)
+
+b
+
$OBJfn.reg
+[1] 0.955787
+
+$NNS.reg.n.best
+[1] 1
+
+$probability.threshold
+[1] 0.6429167
+
+$OBJfn.dim.red
+[1] 0.955787
+
+$NNS.dim.red.threshold
+[1] 0.925
+
+$reg
+ [1] 3 3 3 3 3 3 3 3 3 3
+
+$reg.pred.int
+NULL
+
+$dim.red
+ [1] 3 3 3 3 3 3 3 3 3 3
+
+$dim.red.pred.int
+NULL
+
+$stack
+ [1] 3 3 3 3 3 3 3 3 3 3
+
+$pred.int
+NULL
+
mean( b$stack == as.numeric(iris[test.set, 5]) )
+
+
+
Brief Notes on Other Parameters
+
+depth = "max" will force all observations to be
+their own partition, forcing a perfect fit of the multivariate
+regression. In essence, this is the basis for a kNN nearest
+neighbor type of classification.
n.best = 1 will use the single nearest neighbor.
+When coupled with depth = "max", NNS will
+emulate a kNN = 1 but as the dimensions increase the
+results diverge demonstrating NNS is less sensitive to the
+curse of dimensionality than kNN.
extreme will use the maximum or minimum
+threshold obtained, and may result in errors if that
+threshold cannot be eclipsed by subsequent iterations.
+
+
+
References
+
If the user is so motivated, detailed arguments further examples are
+provided within the following:
+
+
Getting Started with NNS: Forecasting
+Fred Viole
+
+
+
+library(NNS)
+library(data.table)
+require(knitr)
+require(rgl)
+
Forecasting
+
The underlying assumptions of traditional autoregressive models are
+well known. The resulting complexity with these models leads to
+observations such as,
+
``We have found that choosing the wrong model or parameters can
+often yield poor results, and it is unlikely that even experienced
+analysts can choose the correct model and parameters efficiently given
+this array of choices.’’
+
NNS simplifies the forecasting process. Below are some
+examples demonstrating NNS.ARMA and its
+assumption free, minimal parameter forecasting
+method.
+
+
Linear Regression
+
NNS.ARMA has the ability to fit a
+linear regression to the relevant component series, yielding very fast
+results. For our running example we will use the
+AirPassengers dataset loaded in base R.
+
We will forecast 44 periods h = 44 of
+AirPassengers using the first 100 observations
+training.set = 100, returning estimates of the final 44
+observations. We will then test this against our validation set of
+tail(AirPassengers,44).
+
Since this is monthly data, we will try a
+seasonal.factor = 12.
+
Below is the linear fit and associated root mean squared error (RMSE)
+using method = "lin".
+
nns_lin = NNS.ARMA(AirPassengers,
+ h = 44,
+ training.set = 100,
+ method = "lin",
+ plot = TRUE,
+ seasonal.factor = 12,
+ seasonal.plot = FALSE)
+
+
sqrt(mean((nns_lin - tail(AirPassengers, 44)) ^ 2))
+
## [1] 35.39965
+
+
+
Nonlinear Regression
+
Now we can try using a nonlinear regression on the relevant component
+series using method = "nonlin".
+
nns_nonlin = NNS.ARMA(AirPassengers,
+ h = 44,
+ training.set = 100,
+ method = "nonlin",
+ plot = FALSE,
+ seasonal.factor = 12,
+ seasonal.plot = FALSE)
+
+sqrt(mean((nns_nonlin - tail(AirPassengers, 44)) ^ 2))
+
+
+
+
Cross-Validation
+
We can test a series of seasonal.factors and select the
+best one to fit. The largest period to consider would be
+0.5 * length(variable), since we need more than 2 points
+for a regression! Remember, we are testing the first 100 observations of
+AirPassengers, not the full 144 observations.
+
seas = t(sapply(1 : 25, function(i) c(i, sqrt( mean( (NNS.ARMA(AirPassengers, h = 44, training.set = 100, method = "lin", seasonal.factor = i, plot=FALSE) - tail(AirPassengers, 44)) ^ 2) ) ) ) )
+
+colnames(seas) = c("Period", "RMSE")
+seas
+
## Period RMSE
+## [1,] 1 75.67783
+## [2,] 2 75.71250
+## [3,] 3 75.87604
+## [4,] 4 75.16563
+## [5,] 5 76.07418
+## [6,] 6 70.43185
+## [7,] 7 77.98493
+## [8,] 8 75.48997
+## [9,] 9 79.16378
+## [10,] 10 81.47260
+## [11,] 11 106.56886
+## [12,] 12 35.39965
+## [13,] 13 90.98265
+## [14,] 14 95.64979
+## [15,] 15 82.05345
+## [16,] 16 74.63052
+## [17,] 17 87.54036
+## [18,] 18 74.90881
+## [19,] 19 96.96011
+## [20,] 20 88.75015
+## [21,] 21 100.21346
+## [22,] 22 108.68674
+## [23,] 23 85.06430
+## [24,] 24 35.49018
+## [25,] 25 75.16192
+
Now we know seasonal.factor = 12 is our best fit, we can
+see if there’s any benefit from using a nonlinear regression.
+Alternatively, we can define our best fit as the corresponding
+seas$Period entry of the minimum value in our
+seas$RMSE column.
+
a = seas[which.min(seas[ , 2]), 1]
+
Below you will notice the use of seasonal.factor = a
+generates the same output.
+
nns = NNS.ARMA(AirPassengers,
+ h = 44,
+ training.set = 100,
+ method = "nonlin",
+ seasonal.factor = a,
+ plot = TRUE, seasonal.plot = FALSE)
+
+
sqrt(mean((nns - tail(AirPassengers, 44)) ^ 2))
+
## [1] 18.1809
+
Note: You may experience instances with monthly data
+that report seasonal.factor close to multiples of 3, 4, 6
+or 12. For instance, if the reported
+seasonal.factor = {37, 47, 71, 73} use
+(seasonal.factor = c(36, 48, 72)) by setting the
+modulo parameter in
+NNS.seas(..., modulo = 12). The same
+suggestion holds for daily data and multiples of 7, or any other time
+series with logically inferred cyclical patterns. The nearest periods to
+that modulo will be in the expanded output.
+
NNS.seas(AirPassengers, modulo = 12, plot = FALSE)
+
## $all.periods
+## Period Coefficient.of.Variation Variable.Coefficient.of.Variation
+## 1 48 0.4002249 0.4279947
+## 2 12 0.4059923 0.4279947
+## 3 24 0.4279947 0.4279947
+## 4 36 0.4279947 0.4279947
+## 5 60 0.4279947 0.4279947
+##
+## $best.period
+## [1] 48
+##
+## $periods
+## [1] 48 12 24 36 60
+
+
+
Cross-Validating All Combinations of
+seasonal.factor
+
NNS also offers a wrapper function
+NNS.ARMA.optim() to test a given vector of
+seasonal.factor and returns the optimized objective
+function (in this case RMSE written as
+obj.fn = expression( sqrt(mean((predicted - actual)^2)) ))
+and the corresponding periods, as well as the
+NNS.ARMA regression method used.
+Alternatively, using external package objective functions work as well
+such as
+obj.fn = expression(Metrics::rmse(actual, predicted)).
+
NNS.ARMA.optim() will also test whether
+to regress the underlying data first, shrink the estimates
+to their subset mean values, include a bias.shift based on
+its internal validation errors, and compare different
+weights of both linear and nonlinear estimates.
+
Given our monthly dataset, we will try multiple years by setting
+seasonal.factor = seq(12, 60, 6) every 6 months based on
+our NNS.seas() insights above.
+
nns.optimal = NNS.ARMA.optim(AirPassengers,
+ training.set = 100,
+ seasonal.factor = seq(12, 60, 6),
+ obj.fn = expression( sqrt(mean((predicted - actual)^2)) ),
+ objective = "min",
+ pred.int = .95, plot = TRUE)
+
+nns.optimal
+
[1] "CURRNET METHOD: lin"
+[1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
+[1] "NNS.ARMA(... method = 'lin' , seasonal.factor = c( 12 ) ...)"
+[1] "CURRENT lin OBJECTIVE FUNCTION = 35.3996540135277"
+[1] "BEST method = 'lin', seasonal.factor = c( 12 )"
+[1] "BEST lin OBJECTIVE FUNCTION = 35.3996540135277"
+[1] "CURRNET METHOD: nonlin"
+[1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
+[1] "NNS.ARMA(... method = 'nonlin' , seasonal.factor = c( 12 ) ...)"
+[1] "CURRENT nonlin OBJECTIVE FUNCTION = 18.1809033101955"
+[1] "BEST method = 'nonlin' PATH MEMBER = c( 12 )"
+[1] "BEST nonlin OBJECTIVE FUNCTION = 18.1809033101955"
+[1] "CURRNET METHOD: both"
+[1] "COPY LATEST PARAMETERS DIRECTLY FOR NNS.ARMA() IF ERROR:"
+[1] "NNS.ARMA(... method = 'both' , seasonal.factor = c( 12 ) ...)"
+[1] "CURRENT both OBJECTIVE FUNCTION = 22.7363330823967"
+[1] "BEST method = 'both' PATH MEMBER = c( 12 )"
+[1] "BEST both OBJECTIVE FUNCTION = 22.7363330823967"
+>
+> nns.optimal
+$periods
+[1] 12
+
+$weights
+NULL
+
+$obj.fn
+[1] 18.1809
+
+$method
+[1] "nonlin"
+
+$shrink
+[1] FALSE
+
+$nns.regress
+[1] FALSE
+
+$bias.shift
+[1] 0
+
+$errors
+ [1] -6.0626221 -10.8434613 -10.7646998 -22.7134790 -15.3519569 -12.9673866 -9.1626428 3.9393939 7.4882812 12.3750000 29.1132812 34.3281250 19.7002739
+[14] 20.0656989 11.8833952 -15.1389735 24.1108241 7.4289721 15.2385271 38.3826941 19.2903993 17.4644272 19.3331767 19.8155057 -4.0856291 26.3260739
+[27] 2.6153110 -24.3491085 3.9057436 -8.8271346 -7.9236143 5.9867956 -3.9068174 -0.7986170 42.1995863 -10.1324609 -20.0852820 8.6573328 -21.3067790
+[40] -24.3403514 -0.6332912 -29.8418247 -5.8572216 14.8998761
+
+$results
+ [1] 348.9374 411.1565 454.2353 444.2865 388.6480 334.0326 295.8374 339.9394 347.4883 330.3750 391.1133 382.3281 382.7003 455.0657 502.8834 489.8610 428.1108
+[18] 366.4290 325.2385 375.3827 379.2904 359.4644 425.3332 415.8155 415.9144 498.3261 550.6153 534.6509 466.9057 398.1729 354.0764 410.9868 413.0932 390.2014
+[35] 461.1996 450.8675 451.9147 543.6573 600.6932 581.6596 507.3667 431.1582 384.1428 446.8999
+
+$lower.pred.int
+ [1] 310.8588 373.0779 416.1567 406.2079 350.5694 295.9540 257.7588 301.8608 309.4097 292.2964 353.0347 344.2495 344.6217 416.9871 464.8048 451.7824 390.0322
+[18] 328.3504 287.1599 337.3041 341.2118 321.3858 387.2546 377.7369 377.8358 460.2475 512.5367 496.5723 428.8271 360.0943 315.9978 372.9082 375.0146 352.1228
+[35] 423.1210 412.7889 413.8361 505.5787 562.6146 543.5810 469.2881 393.0796 346.0642 408.8213
+
+$upper.pred.int
+ [1] 387.0160 449.2351 492.3139 482.3651 426.7266 372.1112 333.9160 378.0180 385.5669 368.4536 429.1919 420.4067 420.7789 493.1443 540.9620 527.9396 466.1894
+[18] 404.5076 363.3171 413.4613 417.3690 397.5430 463.4118 453.8941 453.9930 536.4047 588.6939 572.7295 504.9843 436.2515 392.1550 449.0654 451.1718 428.2800
+[35] 499.2782 488.9461 489.9933 581.7359 638.7718 619.7382 545.4453 469.2368 422.2214 484.9785
+
+
+
+
+
+
Extension of Estimates
+
We can forecast another 50 periods out-of-sample
+(h = 50), by dropping the training.set
+parameter while generating the 95% prediction intervals.
+
NNS.ARMA.optim(AirPassengers,
+ seasonal.factor = seq(12, 60, 6),
+ obj.fn = expression( sqrt(mean((predicted - actual)^2)) ),
+ objective = "min",
+ pred.int = .95, h = 50, plot = TRUE)
+
+
+
+
+
+
Brief Notes on Other Parameters
+
+seasonal.factor = c(1, 2, ...)
+
We included the ability to use any number of specified seasonal
+periods simultaneously, weighted by their strength of seasonality.
+Computationally expensive when used with nonlinear regressions and large
+numbers of relevant periods.
+
+
Instead of weighting by the seasonal.factor strength of
+seasonality, we offer the ability to weight each per any defined
+compatible vector summing to 1.
+Equal weighting would be weights = "equal".
+
+
Provides the values for the specified prediction intervals within
+[0,1] for each forecasted point and plots the bootstrapped replicates
+for the forecasted points.
+
+seasonal.factor = FALSE
+
We also included the ability to use all detected seasonal periods
+simultaneously, weighted by their strength of seasonality.
+Computationally expensive when used with nonlinear regressions and large
+numbers of relevant periods.
+
+
This parameter restricts the number of detected seasonal periods to
+use, again, weighted by their strength. To be used in conjunction with
+seasonal.factor = FALSE.
+
+
To be used in conjunction with seasonal.factor = FALSE.
+This parameter will ensure logical seasonal patterns (i.e.,
+modulo = 7 for daily data) are included along with the
+results.
+
+
To be used in conjunction with
+seasonal.factor = FALSE & modulo != NULL. This
+parameter will ensure empirical patterns are kept along with the logical
+seasonal patterns.
+
+
This setting generates a new seasonal period(s) using the estimated
+values as continuations of the variable, either with or without a
+training.set. Also computationally expensive due to the
+recalculation of seasonal periods for each estimated value.
+
+
These are the plotting arguments, easily enabled or disabled with
+TRUE or FALSE.
+seasonal.plot = TRUE will not plot without
+plot = TRUE. If a seasonal analysis is all that is desired,
+NNS.seas is the function specifically suited for that
+task.
+
+
+
Multivariate Time Series Forecasting
+
The extension to a generalized multivariate instance is provided in
+the following documentation of the
+NNS.VAR() function:
+
+
+
References
+
If the user is so motivated, detailed arguments and proofs are
+provided within the following:
+
+