added floating point guard

R/graphon_distribution.R

@@ -67,8 +67,7 @@ pgraphon <- function(
   dev_mat <- outer(y, inner_products, "-")
   
   # Apply CDF Fv to each deviation
-  # With sapply we do not have to assume that Fv is vectorized
-  cdf_values <- sapply(dev_mat, Fv)
+  cdf_vals <- Fv(dev_mat)
   
   # Row means give the empirical expectation for each y_j
   out <- rowMeans(cdf_vals)
@@ -140,8 +139,7 @@ dgraphon <- function(
   dev_mat <- outer(y, inner_products, "-")
   
   # Apply the density function to each y[j]
-  # With sapply we do not have to assume that fv is vectorized
-  dens_vals <- sapply(dev_mat, fv)
+  dens_vals <- fv(dev_mat)
   
   # Row means give the empirical expectation for each y_j
   out <- rowMeans(dens_vals)

R/singular_values.R

@@ -44,6 +44,10 @@ source("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 guard Positive numeric guard value.  Default is `sqrt(.Machine$double.eps)`,
+#'   which is about `1.5e‑8` on most platforms – small enough to be negligible
+#'   for most computations. If it is null, then it is not used. 
+#'   The guard is used for the value k = 0, which can cause arithmetic errors.
 #'
 #' @return A numeric matrix **Q** of dimension `K × n`.  The \eqn{j}-th row
 #'         (for `j = 1,…,K`) contains the increments of the CDF evaluated at
@@ -100,7 +104,8 @@ compute_matrix <- function(
     K,
     sample_X_fn,
     fv,
-    Fv
+    Fv,
+    guard = sqrt(.Machine$double.eps)
 ) {
   ## 1.1 Check inputs ==========================================================
   if (!is.numeric(seed) || length(seed) != 1) stop("'seed' must be a single number")
@@ -127,6 +132,9 @@ compute_matrix <- function(
   
   ## 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)
   
   ## 1.5 Build the matrix Q ====================================================

scripts/plot_densities.R

@@ -0,0 +1,177 @@
+# 1. Load functions and Libraries ----------------------------------------------
+source("./R/graphon_distribution.R")
+
+# 2. Set constants -------------------------------------------------------------
+withr::with_seed(42) # set seed in a local environment
+
+a0 <- c(0, 0)
+a1 <- c(1.0)
+a2 <- c(0.5, -0.3)
+
+N <- 400 # number of samples for X
+
+X1 = matrix(rnorm(N))
+X2 = matrix(rnorm(2 * N), ncol = 2)
+
+
+# 3. Normally distributed v's --------------------------------------------------
+## 3.1 Plot a = 0 ==============================================================
+# In this section we plot the conditional density p(u | x) for v ~ N(0,1) and
+# a = (0, 0)
+p1 <- create_cond_density(a0, dnorm, pnorm, X2)
+givenX <- c(0, 0)
+p1_plot <- \(x) p1(x, givenX)
+plot.function(p1_plot, xlab="u", ylab="p(u |x) ")
+title(main=expression("Conditional density for" ~ v %~% N(0, 1)), line=2)
+title(
+  main = bquote(
+    a == (.(paste(a0, collapse = ", "))) ~ "," ~
+      N == .(N) ~ "," ~
+      x == (.(paste(givenX, collapse = ", ")))
+  ),
+  line = 1,
+  cex = 0.65
+)
+
+## 3.2 Plot a != 0, X != 0 =====================================================
+p2 <- create_cond_density(a2, dnorm, pnorm, X2)
+givenX <- c(-0.5, 0.5)
+p2_plot <- \(x) p2(x, givenX)
+plot.function(p2_plot, xlab="u", ylab="p(u |x) ")
+title(main=expression("Conditional density for" ~ v %~% N(0, 1)), line=2)
+title(
+  main = bquote(
+    a == (.(paste(a2, collapse = ", "))) ~ "," ~
+      N == .(N) ~ "," ~
+      x == (.(paste(givenX, collapse = ", ")))
+  ),
+  line = 1,
+  cex = 0.65
+)
+
+## 3.3 Plot one dimensional X_i ================================================
+p3 <- create_cond_density(a1, dnorm, pnorm, X1)
+givenX <- c(-0.5)
+p3_plot <- \(x) p3(x, givenX)
+plot.function(p3_plot, xlab="u", ylab="p(u |x) ")
+title(main=expression("Conditional density for" ~ v %~% N(0, 1)), line=2)
+title(
+  main = bquote(
+    a == (.(paste(a1, collapse = ", "))) ~ "," ~
+      N == .(N) ~ "," ~
+      x == (.(paste(givenX, collapse = ", ")))
+  ),
+  line = 1,
+  cex = 0.65
+)
+
+## 3.4 Plot v ~ N(0, 2) ========================================================
+dens_func <-\(x) dnorm(x, sd=sqrt(2))
+cdf <- \(x) pnorm(x, sd=sqrt(2))
+p4 <- create_cond_density(a1, dens_func, cdf, X1)
+givenX <- c(-0.5)
+p4_plot <- \(x) p4(x, givenX)
+plot.function(p4_plot, xlab="u", ylab="p(u |x) ")
+title(main=expression("Conditional density for" ~ v %~% N(0, 2)), line=2)
+title(
+  main = bquote(
+    a == (.(paste(a1, collapse = ", "))) ~ "," ~
+      N == .(N) ~ "," ~
+      x == (.(paste(givenX, collapse = ", ")))
+  ),
+  line = 1,
+  cex = 0.65
+)
+
+# 4. Expontial distribution for v's --------------------------------------------
+## 4.1 One dimensional X_i =====================================================
+# This example has a jump in [0.14, 0.15].
+# The zero can not be explained and it's not a probability as:
+# integrate(p5_plot, 0.15, 1, subdivisions = 500)
+# => 0.952492 with absolute error < 0.00011
+lambda <- 1.0
+p5 <- create_cond_density(a1, dexp, pexp, X1)
+givenX <- c(-0.5)
+p5_plot <- \(x) p5(x, givenX)
+plot.function(p5_plot, xlab="u", ylab="p(u |x) ")
+title(main=expression("Conditional density for" ~ v %~% Exp(1)), line=2)
+title(
+  main = bquote(
+    a == (.(paste(a1, collapse = ", "))) ~ "," ~
+      N == .(N) ~ "," ~
+      x == (.(paste(givenX, collapse = ", ")))
+  ),
+  line = 1,
+  cex = 0.65
+)
+
+## 4.2 Two dimensional X_i =====================================================
+lambda <- 1.0
+p6 <- create_cond_density(a2, dexp, pexp, X2)
+givenX <- c(-0.5)
+p6_plot <- \(x) p6(x, givenX)
+plot.function(p6_plot, xlab="u", ylab="p(u |x) ")
+title(main=expression("Conditional density for" ~ v %~% Exp(1)), line=2)
+title(
+  main = bquote(
+    a == (.(paste(a2, collapse = ", "))) ~ "," ~
+      N == .(N) ~ "," ~
+      x == (.(paste(givenX, collapse = ", ")))
+  ),
+  line = 1,
+  cex = 0.65
+)
+
+## 4.3 Higher Rate of exponential distribtion ==================================
+lambda <- 2.0
+dens_func <- \(x) dexp(x, rate=lambda)
+cdf <- \(x) pexp(x, rate=lambda)
+p7 <- create_cond_density(a2, dens_func, cdf, X2)
+givenX <- c(-0.5, 0.5)
+p7_plot <- \(x) p7(x, givenX)
+plot.function(p7_plot, xlab="u", ylab="p(u |x) ")
+title(main=expression("Conditional density for" ~ v %~% Exp(2.0)), line=2)
+title(
+  main = bquote(
+    a == (.(paste(a2, collapse = ", "))) ~ "," ~
+      N == .(N) ~ "," ~
+      x == (.(paste(givenX, collapse = ", ")))
+  ),
+  line = 1,
+  cex = 0.65
+)
+
+# 5. Uniform Distribution ------------------------------------------------------
+
+## Two dimensional X_i =========================================================
+p8 <- create_cond_density(a2, dunif, punif, X2)
+givenX <- c(-0.5, 0.5)
+p8_plot <- \(x) p8(x, givenX)
+plot.function(p8_plot, xlab="u", ylab="p(u |x) ")
+title(main=expression("Conditional density for" ~ v %~% U(0,1)), line=2)
+title(
+  main = bquote(
+    a == (.(paste(a2, collapse = ", "))) ~ "," ~
+      N == .(N) ~ "," ~
+      x == (.(paste(givenX, collapse = ", ")))
+  ),
+  line = 1,
+  cex = 0.65
+)
+
+## 5.2 One dimensional X_i =====================================================
+p9 <- create_cond_density(a1, dunif, punif, X1)
+givenX <- c(-0.5)
+p9_plot <- \(x) p9(x, givenX)
+plot.function(p9_plot, xlab="u", ylab="p(u |x) ")
+title(main=expression("Conditional density for" ~ v %~% U(0,1)), line=2)
+title(
+  main = bquote(
+    a == (.(paste(a1, collapse = ", "))) ~ "," ~
+      N == .(N) ~ "," ~
+      x == (.(paste(givenX, collapse = ", ")))
+  ),
+  line = 1,
+  cex = 0.65
+)
+