neuster Stand

R/estimators.R

@@ -7,9 +7,9 @@ source(here::here("R", "build_network.R"))
 # Helper functions -------------------------------------------------------------
 # helper function for wrapping the parameters of the Q_a creation function
 # TODO rename this function
-make_matrix_creation <- function(seed, n, K, matrix_X, fv, Fv, guard) {
+make_matrix_creation <- function(seed, n, K, matrix_X, fv, Fv, guard, fX=NULL) {
   function(a) {
-    compute_matrix(seed=seed, a, n=n, K=K, matrix_X = matrix_X, fv=fv, Fv=Fv, guard=guard)
+    compute_matrix(seed=seed, a, n=n, K=K, matrix_X = matrix_X, fv=fv, Fv=Fv, guard=guard, fX=fX)
   }
 }
 
@@ -200,7 +200,7 @@ calculate_edge_density <- function(adj_matrix) {
   return(rho)
 }
 # test the estimator routines
-seed <- 17L #  121L this seed works exceptionally well
+seed <- 121L #  121L this seed works exceptionally well
 set.seed(seed)
 #X  <- matrix(seq(-1, 1, length.out = 5), ncol = 1)
 a  <- 20
@@ -226,8 +226,8 @@ adj
 
 # Q_a matrix
 Qa <- compute_matrix(seed, a=a, n=n, K=K, fv=fv, Fv=Fv, guard=guard, matrix_X=X)
-
+Qa2 <- compute_matrix(seed, a=a, n=n, K=K, fv=fv, Fv=Fv, guard=guard, matrix_X =X, fX= dnorm)
-calc_Q_a <- make_matrix_creation(seed, n=n, K=K, matrix_X = X, fv=fv, Fv=Fv, guard=guard)
+calc_Q_a <- make_matrix_creation(seed, n=n, K=K, matrix_X = X, fv=fv, Fv=Fv, guard=guard, fX=dnorm)
 
 loss_func <- function(a) {
   Q_a <- calc_Q_a(a)
@@ -236,7 +236,6 @@ loss_func <- function(a) {
   norm(pinv_Qa %*% Q_a %*% adj %*% pinv_Qa %*% Q_a - adj, type="F")^2
 }
 
-
 plot_as <- seq(-10, 100, length.out=500)
 loss_vals <- sapply(plot_as, loss_func)
 plot(plot_as , loss_vals, type="b")

R/graphon_distribution.R

@@ -2,7 +2,67 @@
 # using the here package.
 source(here::here("R", "qinf.R"))
 
-
+# 0. Analytical test case for d = 1 -----------------------------
+#' --------------------------------------------------------------
+#'  Fa(y) = ∫ Fv(y - z) fZ(z) dz   with  Z = a * X
+#' --------------------------------------------------------------
+#' Compute the graphon CDF for a scalar X and scalar coefficient a.
+#'
+#' @param y   numeric vector of points where Fa(y) is required.
+#' @param a   numeric scalar (coefficient).  a != 0.
+#' @param fX  function(z) returning the pdf of X (vectorised).
+#' @param Fv  function(x) returning the CDF of the noise v (vectorised).
+#' @param lower,upper numeric bounds for the integration over z.
+#'        By default they are set to the (practical) support of X
+#'        transformed by a.  You can tighten them for speed.
+#' @param rel.tol relative tolerance passed to `integrate`.
+#' @return numeric vector of the same length as y.
+#' @examples
+#' ## Example 1: X ~ Gamma(2,1), a = 0.8, v ~ N(0,1)
+#' fX  <- function(z) dgamma(z, shape = 2, rate = 1)   # pdf of X
+#' Fv  <- function(x) pnorm(x, mean = 0, sd = 1)      # cdf of v
+#' a   <- 0.8
+#' y   <- seq(-3, 6, length.out = 200)
+#' Fa_vals <- Fa_one_dim(y, a, fX, Fv)
+#' plot(y, Fa_vals, type = "l", col = "steelblue",
+#'      main = "Fa(y) for Gamma X, Normal noise",
+#'      xlab = "y", ylab = expression(F[a](y)))
+#'
+#' ## Example 2: both X and v are normal – analytic check
+#' fX2 <- function(z) dnorm(z, mean = 1, sd = 2)
+#' Fv2 <- function(x) pnorm(x, mean = 0, sd = 1)
+#' a2  <- -1.5
+#' y2  <- seq(-5, 8, length.out = 100)
+#' Fa2_num <- Fa_one_dim(y2, a2, fX2, Fv2)
+#' ## analytic result: Y = a*X + v ~ N(a*mu_X, a^2*sd_X^2 + sd_v^2)
+#' muY  <- a2 * 1
+#' sdY  <- sqrt(a2^2 * 2^2 + 1^2)
+#' Fa2_ana <- pnorm(y2, mean = muY, sd = sdY)
+#' max(abs(Fa2_num - Fa2_ana))   # should be ~ 1e-6
+#' --------------------------------------------------------------
+pgraphon_analytical_1d <- function(y, a, fX, Fv,
+                       lower = -Inf, upper = Inf,
+                       rel.tol = .Machine$double.eps^0.5) {
+  stopifnot(is.numeric(a), length(a) == 1, a != 0)
+  stopifnot(is.function(fX), is.function(Fv))
+  stopifnot(is.numeric(y))
+  
+  # pdf of Z = a*X
+  fZ <- function(z) {
+    # Jacobian |a|^{-1} and argument scaling
+    (1 / abs(a)) * fX(z / a)
+  }
+  
+  # Helper that integrates for a single y
+  one_y <- function(yy) {
+    integrand <- function(z) Fv(yy - z) * fZ(z)
+    integrate(integrand, lower = lower, upper = upper,
+              rel.tol = rel.tol)$value
+  }
+  
+  # Vectorised over the whole y vector
+  vapply(y, one_y, numeric(1))
+}
 # 1. Distribution Function -----------------------------------------------------
 
 #' Empirical Distribution function for the graphon model.
@@ -148,108 +208,110 @@ dgraphon <- function(
 
 # 3. Quantile Function ---------------------------------------------------------
 
-#' Quantile of the empirical graphon distribution
+#' Quantile of the (empirical or analytic) graphon distribution
 #'
-#' This is a thin wrapper around the generic infimum‑type quantile routine
+#' If the covariate matrix has a single column and the coefficient vector
-#' \code{qinf()}.  It builds the CDF of the graphon,
+#' has length 1, the function can use an *analytic* CDF (provided the user
-#' \eqn{p_{\text{graphon}}(x) = \frac{1}{n}\sum_{i=1}^{n}F_v\bigl(x-a^\top X_i\bigr)},
+#' supplies the pdf of the covariate via `fX`).  Otherwise it falls back to the
-#' and then finds the quantile(s) for the supplied probability(ies) \code{p}.
+#' original empirical estimator `pgraphon`.
-#'
-#' @param p Numeric vector of probabilities in \eqn{[0,1]}.  Values outside this
-#'   interval are rejected.
-#' @param a Numeric coefficient vector.  Its length must equal the number of
-#'   columns of \code{X_matrix}.
-#' @param Fv Function.  The CDF of the latent variable $v$.  It must be
-#'   vectorised (i.e. accept a numeric vector and return a numeric vector of
-#'   the same length).  Typical examples are \code{pnorm}, \code{pexp}, etc.
-#' @param X_matrix Numeric matrix of dimension \eqn{n\times p}.  Each row
-#'   corresponds to an observation $X_i$.
-#' @param lower,upper Numeric scalars giving the search interval for the root
-#'   finder.  By default they are set to \code{-Inf} and \code{Inf}.
-#' @param tol Numeric tolerance for the bisection algorithm used inside
-#'   \code{qinf()}.  The default is the square‑root of machine epsilon.
-#' @param max.iter Maximum number of bisection iterations (safety guard).
-#' @param ... Additional arguments that are passed **directly** to the
-#'   internal CDF \code{pgraphon}.  This makes the wrapper flexible – you do not
-#'   have to list every argument (e.g. \code{a}, \code{X_matrix}, \code{Fv})
-#'   explicitly.
-#'
-#' @return A numeric vector of the same length as \code{p} containing the
-#'   quantiles.  The vector is named by the probabilities.
 #'
+#' @param p Numeric vector of probabilities in \eqn{[0,1]}.
+#' @param a Numeric coefficient vector (length = ncol(X_matrix)).
+#' @param Fv Function. CDF of the latent variable $v$ (vectorised).
+#' @param X_matrix Numeric matrix of size $n\times p$; each row is an observation.
+#' @param fX Optional function(z) returning the *known* pdf of the *scalar*
+#'            covariate when `ncol(X_matrix)==1 && length(a)==1`.  If supplied,
+#'            the analytic CDF `pgraphon_analytical_1d` is used.
+#' @param lower,upper Numeric bounds for the root‑finding algorithm.
+#' @param tol Numeric tolerance for the bisection inside `qinf`.
+#' @param max.iter Maximum number of bisection iterations.
+#' @param ... Additional arguments passed to the internal CDF (`pgraphon` or the
+#'            analytic version).  For the analytic case you may want to pass
+#'            `lower`, `upper`, or `rel.tol`.
+#' @return Numeric vector of quantiles, named by the probabilities.
 #' @examples
-#' ## ---- simple normal example ------------------------------------------------
+#' ## 1. Empirical case
 #' set.seed(123)
-#' X <- matrix(rnorm(200), ncol = 2)          # n = 100, p = 2
+#' X <- matrix(rnorm(200), ncol = 2)
 #' a <- c(0.7, -0.3)
-#' qgraphon(p = c(0.25, 0.5, 0.75),
+#' qgraphon(p = c(0.25, 0.5, 0.75), a = a,
-#'          a = a,
+#'          Fv = pnorm, X_matrix = X)
-#'          Fv = pnorm,
-#'          X_matrix = X)
-#'
-#' ## ---- mixture example ------------------------------------------------------
-#' mix_cdf <- function(x, w = 0.4, mu = 0, sigma = 1) {
-#'   w * (x >= 0) + (1 - w) * pnorm(x, mean = mu, sd = sigma)
-#' }
-#' qgraphon(p = seq(0, 1, 0.2),
-#'          a = a,
-#'          Fv = mix_cdf,
-#'          X_matrix = X)
 #'
+#' ## 2. Analytic 1‑D case
+#' ## X ~ Gamma(2,1), a = 0.8, v ~ N(0,1)
+#' fX  <- function(z) dgamma(z, shape = 2, rate = 1)
+#' a1  <- 0.8
+#' X1  <- matrix(rnorm(100), ncol = 1)   # dummy – not used by the analytic branch
+#' qgraphon(p = c(0.1, 0.5, 0.9), a = a1,
+#'          Fv = pnorm, X_matrix = X1,
+#'          fX = fX, lower = -10, upper = 20)
 #' @export
-qgraphon <- function(
+qgraphon <- function(p, a, Fv, X_matrix,
-    p,                               # probabilities
+                     fX = NULL,               # <- new optional argument
-    a,                               # coefficient vector
+                     lower = -Inf, upper = Inf,
-    Fv,                              # CDF of the v's
+                     tol = .Machine$double.eps^0.5,
-    X_matrix,                         # X_i samples, matrix
+                     max.iter = 100, ...) {
-    lower = -Inf,                     # lower bound roots
-    upper =  Inf,                     # upper bound roots
-    tol   = .Machine$double.eps^0.5,  # parameter root finding
-    max.iter = 100) {
   
-  ## 3.1 Check inputs ----------------------------------------------------------
+  ## ---- input checks -------------------------------------------------
   if (!is.numeric(p) || any(is.na(p))) {
     stop("'p' must be a numeric vector without NA values")
   }
   if (any(p < 0 | p > 1)) {
     stop("All probabilities in 'p' must lie in [0, 1]")
   }
-  
   if (!is.numeric(a) || !is.vector(a)) {
     stop("'a' must be a numeric vector")
   }
-  
   if (!is.matrix(X_matrix) || !is.numeric(X_matrix)) {
     stop("'X_matrix' must be a numeric matrix")
   }
-  
   if (ncol(X_matrix) != length(a)) {
     stop("Number of columns of 'X_matrix' (", ncol(X_matrix),
          ") must equal length of 'a' (", length(a), ")")
   }
-  
   if (!is.function(Fv)) {
     stop("'Fv' must be a function (the CDF of the latent variable)")
   }
   
-  ## 3.2 Call the generic quantile function ------------------------------------
+  ## ---- decide which CDF to use --------------------------------------
-  out <- qinf(F = pgraphon,
+  use_analytic <- (ncol(X_matrix) == 1L) && (length(a) == 1L) && (!is.null(fX))
-            p = p,
+  
-            lower = lower,
+  if (use_analytic) {
-            upper = upper,
+    # Build a *wrapper* that has the same signature as pgraphon()
-            tol = tol,
+    # (y, a, Fv, X_matrix, ...) so that qinf() can call it unchanged.
-            max.iter = max.iter,
+    analytic_cdf <- function(y, a, Fv, X_matrix, ...) {
-            a = a,
+      # X_matrix is ignored – the distribution of X is supplied via fX
-            X_matrix = X_matrix,
+      pgraphon_analytical_1d(y = y,
-            Fv = Fv)
+                             a = a,
+                             fX = fX,
+                             Fv = Fv,
+                             lower = lower,
+                             upper = upper,
+                             rel.tol = tol,
+                             ...)
+    }
+    CDF_to_use <- analytic_cdf
+  } else {
+    # Fall back to the empirical estimator 
+    CDF_to_use <- pgraphon
+  }
+  
+  ## ---- call the generic quantile routine ----------------------------
+  out <- qinf(F = CDF_to_use,
+              p = p,
+              lower = lower,
+              upper = upper,
+              tol = tol,
+              max.iter = max.iter,
+              a = a,
+              X_matrix = X_matrix,
+              Fv = Fv,
+              ...)   # forward any extra arguments (e.g., rel.tol)
   
-  # Name the result by the probabilities – this mirrors the behaviour of
-  # base‑R quantile functions (e.g. qnorm, qbeta).
   names(out) <- as.character(p)
   out
 }
 
-
 # 4. Create conditional density ------------------------------------------------
 
 #' Create a Conditional Density Function for the Graphon Model.

R/singular_values.R

@@ -40,6 +40,9 @@ source(here::here("R", "graphon_distribution.R"))
 #' @param Fv     Cumulative distribution function of the latent variable
 #'               \eqn{v}.  Also has to be vectorised.  Typical examples are
 #'               `pnorm`, `pexp`, ….
+#' @param fX Optional function(z) returning the *known* pdf of the *scalar*
+#'            covariate when `ncol(X_matrix)==1 && length(a)==1`.  If supplied,
+#'            the analytic CDF `pgraphon_analytical_1d` is used.
 #' @param sample_X_fn
 #'               Function with a single argument `n`.  It must return an
 #'               \eqn{n \times p} matrix (or an object coercible to a matrix) of
@@ -111,6 +114,7 @@ compute_matrix <- function(
     K,
     fv,
     Fv,
+    fX=NULL,
     sample_X_fn=NULL,
     matrix_X = NULL,
     guard = sqrt(.Machine$double.eps),
@@ -127,6 +131,7 @@ compute_matrix <- function(
   if (!is.null(matrix_X) && !is.matrix(matrix_X)) stop("matrix_X must be either null or a matrix")
   if (is.null(matrix_X) && is.null(sample_X_fn)) stop("Either 'matrix_X' or 'sample_X_fn' must be supplied!")
   if (!is.null(matrix_X) && !is.null(sample_X_fn)) warning("Both arguments 'matrix_X' and `sample_X_fn` is given. Priority is given by to the first!")
+  if (!is.null(fX) && !is.function(fX)) stop("'fX' must be a density function")
   
   ## 1.2 Generate the Matrix X of covariates ===================================
   # If the argument matrix_X is present, use this matrix, otherwise generate one
@@ -146,14 +151,18 @@ compute_matrix <- function(
   }
   
   ## 1.3 Create conditional density ============================================
-  empir_cond_density <- create_cond_density(a, fv, Fv, X)
+  # this is not used in the computation
+  # empir_cond_density <- create_cond_density(a, fv, Fv, X)
   
   ## 1.4 Compute the graphon quantiles =========================================
   k <- seq(0, K) / K
   if (!is.null(guard)) {
     k[1] <- guard
   }
-  graphon_quantiles <- qgraphon(k, a = a, Fv = Fv, X_matrix = X)
+  # here there is an automatic switch included, if fX is not null and we have a
+  # scalar case, then qpgrahon automatically switches to the analytical 
+  # expression. The intended use is for small values of n
+  graphon_quantiles <- qgraphon(k, a = a, Fv = Fv, X_matrix = X, fX= fX)
   
   ## 1.5 Build the matrix Q ====================================================
   inner_products = as.vector(X %*% a)