From 37d9682e27bbd25079ac3cb1af5c97bf27f36dc3 Mon Sep 17 00:00:00 2001 From: Niclas Date: Mon, 29 Jun 2026 14:46:22 +0200 Subject: [PATCH] neuster Stand --- R/estimators.R | 11 +-- R/graphon_distribution.R | 202 +++++++++++++++++++++++++-------------- R/singular_values.R | 13 ++- 3 files changed, 148 insertions(+), 78 deletions(-) diff --git a/R/estimators.R b/R/estimators.R index e977c94..3e65ae1 100644 --- a/R/estimators.R +++ b/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) - -calc_Q_a <- make_matrix_creation(seed, n=n, K=K, matrix_X = X, fv=fv, Fv=Fv, guard=guard) +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, 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") diff --git a/R/graphon_distribution.R b/R/graphon_distribution.R index 72fedea..12fe35e 100644 --- a/R/graphon_distribution.R +++ 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 -#' \code{qinf()}. It builds the CDF of the graphon, -#' \eqn{p_{\text{graphon}}(x) = \frac{1}{n}\sum_{i=1}^{n}F_v\bigl(x-a^\top X_i\bigr)}, -#' and then finds the quantile(s) for the supplied probability(ies) \code{p}. -#' -#' @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. +#' If the covariate matrix has a single column and the coefficient vector +#' has length 1, the function can use an *analytic* CDF (provided the user +#' supplies the pdf of the covariate via `fX`). Otherwise it falls back to the +#' original empirical estimator `pgraphon`. #' +#' @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), -#' a = a, -#' 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) +#' qgraphon(p = c(0.25, 0.5, 0.75), a = a, +#' Fv = pnorm, 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( - p, # probabilities - a, # coefficient vector - Fv, # CDF of the v's - X_matrix, # X_i samples, matrix - lower = -Inf, # lower bound roots - upper = Inf, # upper bound roots - tol = .Machine$double.eps^0.5, # parameter root finding - max.iter = 100) { +qgraphon <- function(p, a, Fv, X_matrix, + fX = NULL, # <- new optional argument + lower = -Inf, upper = Inf, + tol = .Machine$double.eps^0.5, + 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 ------------------------------------ - out <- qinf(F = pgraphon, - p = p, - lower = lower, - upper = upper, - tol = tol, - max.iter = max.iter, - a = a, - X_matrix = X_matrix, - Fv = Fv) + ## ---- decide which CDF to use -------------------------------------- + use_analytic <- (ncol(X_matrix) == 1L) && (length(a) == 1L) && (!is.null(fX)) + + if (use_analytic) { + # Build a *wrapper* that has the same signature as pgraphon() + # (y, a, Fv, X_matrix, ...) so that qinf() can call it unchanged. + analytic_cdf <- function(y, a, Fv, X_matrix, ...) { + # X_matrix is ignored – the distribution of X is supplied via fX + pgraphon_analytical_1d(y = y, + 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. diff --git a/R/singular_values.R b/R/singular_values.R index 339f141..3503228 100644 --- a/R/singular_values.R +++ b/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)