diff --git a/R/build_network.R b/R/build_network.R
index 567e230..2a80bc8 100644
--- a/R/build_network.R
+++ b/R/build_network.R
@@ -30,7 +30,8 @@ if (requireNamespace("here", quietly = TRUE)) {
 #' @param a A numeric vector of length *d* (coefficients).  
 #'   Must be the same length as the number of columns of `X_matrix`.
 #' @param phi A binary function `phi(x, y)` returning a numeric value.  
-#'   It should be symmetric (`phi(x, y) == phi(y, x)`) and bounded in $[0,1]$.
+#'   It should be symmetric (`phi(x, y) == phi(y, x)`) and integrate to $1$ over 
+#'   $[0,1]^2$.
 #' @param rho_n A numeric scalar in $[0,1]$ that scales the graphon.
 #' @param Fv A cumulative distribution function (CDF) used to transform the
 #'   linear predictor.  Defaults to the standard normal CDF `pnorm`.  
diff --git a/R/create_Q.R b/R/create_Q.R
new file mode 100644
index 0000000..7cad462
--- /dev/null
+++ b/R/create_Q.R
@@ -0,0 +1,52 @@
+create_matrix_Q <- function(
+    qgraphon,
+    a,
+    K,
+    Fv,
+    matrix_X = NULL,
+    scaled = FALSE
+) {
+  ## 1.1 Check inputs ==========================================================
+  if (!is.numeric(a) || !is.vector(a)) stop("'a' must be a numeric vector")
+  if (!is.numeric(K) || length(K) != 1 || K <= 0) stop("'K' must be a positive integer")
+
+  if (!is.function(Fv)) stop("'F_v' must be a function")
+  if (!is.matrix(matrix_X)) stop("matrix_X must be a matrix")
+  if (!is.logical(scaled)) stop("`scaled` must be a logical!")
+  if (ncol(matrix_X) != length(a)) {
+    stop("Number of columns of `matrix_X` (", ncol(matrix_X), ") must equal length(a) (", length(a), ")")
+  }
+  
+  ## 1.3 Compute the graphon quantiles =========================================
+  k <- seq(0, K) / K
+  n <- nrow(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)
+  
+  ## 1.4 Build the matrix Q ====================================================
+  inner_products = as.vector(matrix_X %*% a)
+  
+  # outer(y, x, "-")  gives a matrix with entry (j,i) = y[j] - x[i]
+  #   then we apply the CDF `F_v` to the whole matrix at once.
+  #   finally we take the difference of successive rows (j) to obtain the
+  #   increments required by equation (3.1).
+  cdf_mat <- Fv(outer(graphon_quantiles, inner_products, "-")) # (K +1) x n matrix
+  Q <- diff(cdf_mat, lag=1)                           # operates along rows
+  
+  if (scaled) { Q <- 1 / sqrt(n) * Q }
+  Q
+}
+
+source(here::here("R", "distributionfunctions.R"))
+
+n <- 100
+K <- 3
+a <- c(2.0, -0.5)
+X <- matrix(rnorm(2 * n), nrow = n, ncol = 2)
+Fv <- function(x) {dnorm(x, mean=0, sd=1)}
+
+qgraphon <- make_distribution_func(a=a, Fv=Fv, X_matrix=X)
+Q <- create_matrix_Q(qgraphon, a, K, Fv, X)
+Q
diff --git a/R/distributionfunctions.R b/R/distributionfunctions.R
new file mode 100644
index 0000000..0686e8e
--- /dev/null
+++ b/R/distributionfunctions.R
@@ -0,0 +1,338 @@
+# Function for creating functions
+
+source(here::here("R", "qinf.R"))
+source(here::here("R", "graphon_distribution.R"))
+
+#' Factory for a Graphon CDF (empirical or analytical)
+#'
+#' `make_distribution_func()` builds a **callable CDF** that can be used in
+#' downstream code (e.g., for simulation, inference, or plotting).  
+#' The function works in two mutually exclusive modes:
+#'
+#' * **Empirical mode** – the CDF is estimated from observed covariates
+#'   `X_matrix` and a noise CDF `Fv`.  For each evaluation point `y` the
+#'   expectation \eqn{E[F_v(y - a^\top X_i)]} is approximated by the empirical
+#'   average over the rows of `X_matrix`.
+#' * **Analytical mode** – the user supplies a pre‑computed analytical
+#'   function `Fa(y)` (e.g., a closed‑form expression for \eqn{F_a(y) = P( X\cdot a + v \leq y)}.
+#'   The returned function simply forwards its argument to `Fa`.
+#'
+#' Exactly **one** of the two mode specifications must be provided; otherwise
+#' an informative error is thrown.
+#'
+#' @param a Numeric vector (or scalar) of coefficients. Required only in the
+#'   empirical mode. Its length must match the number of columns of
+#'   `X_matrix`.
+#' @param X_matrix Numeric matrix (`n × p`) of observed covariates. Required
+#'   only in the empirical mode. Each row corresponds to an observation
+#'   \eqn{X_i}.
+#' @param Fv Function. Cumulative distribution function of the noise variable
+#'   `v` (e.g., `pnorm`, `pexp`). Must be vectorised and is required only in
+#'   the empirical mode.
+#' @param Fa Function. Analytical expression for the graphon CDF
+#'   \eqn{F_a(y)}. Required only in the analytical mode. Must be vectorised.
+#'
+#' @return A **function** `cdf(y)` that evaluates the graphon CDF at the
+#'   supplied numeric vector or scalar `y`. The returned function is fully
+#'   vectorised and performs input validation on each call.
+#'
+#' @details
+#' * **Empirical mode**  
+#'   The inner products \eqn{a^\top X_i} are pre‑computed once for efficiency.
+#'   For a vector of evaluation points `y`, the function builds an
+#'   \eqn{|y| \times n} matrix of deviations `y - a^\top X_i`, applies the noise
+#'   CDF `Fv` element‑wise, and returns the row means (i.e., the empirical
+#'   expectation).  
+#' * **Analytical mode**  
+#'   The supplied `Fa` is assumed to already implement the desired CDF; the
+#'   wrapper simply validates its input and forwards the call via `vapply`
+#'   to guarantee a numeric return of the same length as `y`.
+#'
+#' @examples
+#' ## -----------------------------------------------------------------
+#' ## Empirical mode ----------------------------------------------------
+#' set.seed(42)
+#' X <- matrix(rnorm(200), ncol = 2)   # 100 observations, p = 2
+#' a <- c(1, -0.5)
+#' emp_cdf <- make_distribution_func(a = a, X_matrix = X, Fv = pnorm)
+#' emp_cdf(seq(-3, 3, length.out = 50))
+#'
+#' ## Analytical mode --------------------------------------------------
+#' ## Example: closed‑form CDF for Y = a*X + V where X ~ Gamma(2,1),
+#' ## V ~ N(0,1) and a = 0.8
+#' fX  <- function(z) dgamma(z, shape = 2, rate = 1)   # pdf of X
+#' Fv  <- function(x) pnorm(x, mean = 0, sd = 1)      # CDF of V
+#' a2  <- 0.8
+#' ## Analytic CDF (normal convolution)
+#' Fa  <- function(y) {
+#'   muY <- a2 * 2          # E[X] = shape/rate = 2
+#'   sdY <- sqrt(a2^2 * 2 + 1)  # Var(Y) = a^2*Var(X) + Var(V)
+#'   pnorm(y, mean = muY, sd = sdY)
+#' }
+#' ana_cdf <- make_distribution_func(Fa = Fa)
+#' ana_cdf(seq(-3, 6, length.out = 100))
+#'
+#' @export
+make_distribution_func <- function(a=NULL, X_matrix=NULL, Fv=NULL, Fa=NULL){
+  # ----------------------------------------------------------------
+  # 1. Switch between empirical and analytical mode
+  # ---------------------------------------------------------------
+  analytical_mode <- !is.null(Fa)
+  empirical_mode <- all(!is.null(X_matrix), !is.null(Fv), !is.null(a))
+  
+  if(xor(analytical_mode, empirical_mode) == FALSE){
+    stop("Provide either the arguments for the empiral mode or analytical mode.
+         See documentation for details.")
+  }
+  
+  # ------------------------------------------------------
+  # Empirical mode
+  # ------------------------------------------------------
+  if(empirical_mode) {
+    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 in `X_matrix` must equal length of `a`.")
+    }
+    if (!is.function(Fv)) {
+      stop("'Fv' must be a function (a CDF).")
+    }
+    
+    # Pre‑compute the inner products aᵀX_i once (vectorised)
+    inner_prod <- as.vector(X_matrix %*% a)   # length n
+    
+    # Return a function that only needs y
+    cdf <- function(y) {
+      if (!is.numeric(y)) {
+        stop("'y' must be numeric (scalar or vector).")
+      }
+      # Compute matrix of y - aᵀX_i
+      dev_mat <- outer(y, inner_prod, "-")
+      # Apply the noise CDF and average over rows
+      rowMeans(Fv(dev_mat))
+    }
+    
+    return(cdf)
+  } else {
+    # -------------------------------------------------
+    # Analytical mode
+    # ------------------------------------------------
+    # input validation
+
+    if (!is.function(Fa)){
+      stop("`Fa` must be a function!")
+    }
+    
+    cdf <- function(y) {
+      if (!is.numeric(y)) {
+        stop("'y' must be numeric (scalar or vector).")
+      }
+      vapply(y, Fa, numeric(1))
+    }
+    
+    return(cdf)
+  }
+}
+
+
+#' Factory for a Graphon **density** function (empirical or analytical)
+#'
+#' `make_density_func()` builds a **callable density** that can be used in
+#' simulation, likelihood evaluation, or plotting.  The function operates in
+#' two mutually exclusive modes:
+#'
+#' * **Empirical mode** – the density is estimated from observed covariates
+#'   `X_matrix` and a noise density `fv`.  For each evaluation point `y` the
+#'   expectation \eqn{E[f_v(y - a^\top X_i)]} is approximated by the empirical
+#'   average over the rows of `X_matrix`.
+#' * **Analytical mode** – the user supplies a pre‑computed analytical density
+#'   function `fa(y)`.  The returned function simply forwards its argument to
+#'   `fa`.
+#'
+#' Exactly **one** of the two mode specifications must be provided; otherwise
+#' an informative error is raised.
+#'
+#' @param a Numeric vector (or scalar) of coefficients. Required only in the
+#'   empirical mode. Its length must match the number of columns of
+#'   `X_matrix`.
+#' @param X_matrix Numeric matrix (`n × p`) of observed covariates. Required
+#'   only in the empirical mode. Each row corresponds to an observation
+#'   \eqn{X_i}.
+#' @param fv Function. Density of the noise variable `v` (e.g., `dnorm`,
+#'   `dgamma`). Must be vectorised and is required only in the empirical mode.
+#' @param fa Function. Analytical expression for the graphon density
+#'   \eqn{f_a(y)}. Required only in the analytical mode. Must be vectorised.
+#'
+#' @return A **function** `dens(y)` that evaluates the graphon density at the
+#'   supplied numeric vector or scalar `y`.  The returned function is fully
+#'   vectorised and performs input validation on each call.
+#'
+#' @details
+#' * **Empirical mode**  
+#'   The inner products \eqn{a^\top X_i} are pre‑computed once for efficiency.
+#'   For a vector of evaluation points `y`, the function builds an
+#'   \eqn{|y| \times n} matrix of deviations `y - a^\top X_i`, applies the noise
+#'   density `fv` element‑wise, and returns the row means (i.e., the empirical
+#'   expectation).  
+#' * **Analytical mode**  
+#'   The supplied `fa` is assumed to already implement the desired density; the
+#'   wrapper validates its input and forwards the call via `vapply` to guarantee
+#'   a numeric return of the same length as `y`.
+#'
+#' @examples
+#' ## -----------------------------------------------------------------
+#' ## Empirical mode ----------------------------------------------------
+#' set.seed(123)
+#' X <- matrix(rnorm(200), ncol = 2)   # 100 observations, p = 2
+#' a <- c(1, -0.5)
+#' ## Noise density: standard normal
+#' emp_dens <- make_density_func(a = a, X_matrix = X, fv = dnorm)
+#' emp_dens(seq(-3, 3, length.out = 50))
+#'
+#' ## Analytical mode --------------------------------------------------
+#' ## Example: closed‑form density for Y = a*X + V where
+#' ## X ~ Gamma(2,1), V ~ N(0,1) and a = 0.8
+#' a2  <- 0.8
+#' fa  <- function(y) {
+#'   ## Convolution of Gamma and Normal does not have a simple closed form,
+#'   ## but for illustration we use the normal approximation:
+#'   muY <- a2 * 2                # E[X] = shape/rate = 2
+#'   sdY <- sqrt(a2^2 * 2 + 1)    # Var(Y) = a^2*Var(X) + Var(V)
+#'   dnorm(y, mean = muY, sd = sdY)
+#' }
+#' ana_dens <- make_density_func(fa = fa)
+#' ana_dens(seq(-3, 6, length.out = 100))
+#'
+#' @export
+make_density_func <- function(a=NULL, X_matrix=NULL, fv=NULL, fa=NULL){
+  # ----------------------------------------------------------------
+  # 1. Switch between empirical and analytical mode
+  # ---------------------------------------------------------------
+  analytical_mode <- !is.null(fa)
+  empirical_mode <- all(!is.null(X_matrix), !is.null(fv), !is.null(a))
+  
+  if(xor(analytical_mode, empirical_mode) == FALSE){
+    stop("Provide either the arguments for the empiral mode or analytical mode.
+         See documentation for details.")
+  }
+  
+  # ------------------------------------------------------
+  # Empirical mode
+  # ------------------------------------------------------
+  if(empirical_mode) {
+    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 in `X_matrix` must equal length of `a`.")
+    }
+    if (!is.function(Fv)) {
+      stop("'Fv' must be a function (a density).")
+    }
+    
+    # Pre‑compute the inner products aᵀX_i once (vectorised)
+    inner_prod <- as.vector(X_matrix %*% a)   # length n
+    
+    # Return a function that only needs y
+    cdf <- function(y) {
+      if (!is.numeric(y)) {
+        stop("'y' must be numeric (scalar or vector).")
+      }
+      # Compute matrix of y - aᵀX_i
+      dev_mat <- outer(y, inner_prod, "-")
+      # Apply the noise CDF and average over rows
+      rowMeans(fv(dev_mat))
+    }
+    
+    return(cdf)
+  } else {
+    # -------------------------------------------------
+    # Analytical mode
+    # ------------------------------------------------
+    # input validation
+    
+    if (!is.function(fa)){
+      stop("`fa` must be a function!")
+    }
+    
+    cdf <- function(y) {
+      if (!is.numeric(y)) {
+        stop("'y' must be numeric (scalar or vector).")
+      }
+      vapply(y, fa, numeric(1))
+    }
+    
+    return(cdf)
+  }
+}
+
+#' Factory for a Quantile Function from a Cumulative Distribution Function
+#'
+#' `make_quantile_function()` builds a **quantile (inverse CDF) function** by
+#' wrapping the root‑finding routine `qinf()`.  The returned function accepts a
+#' probability (or a vector of probabilities) `p` and returns the corresponding
+#' quantile(s) of the distribution defined by the supplied CDF `cdf`.
+#'
+#' @param cdf Function.  A cumulative distribution function that takes a numeric
+#'   argument `x` and returns \eqn{F(x) \in [0,1]}.  The function must be
+#'   **monotone non‑decreasing** and vectorised.
+#' @param lower Numeric.  Lower bound of the search interval for the root‑finding
+#'   algorithm.  Defaults to $-\infty$ (practically the left‑most support of
+#'   the distribution).
+#' @param upper Numeric.  Upper bound of the search interval.  Defaults to
+#'   $\infty$ (right‑most support).
+#' @param tol Numeric.  Relative tolerance passed to `qinf()`.  The default is
+#'   `sqrt(.Machine$double.eps)`, which provides double‑precision accuracy.
+#' @param max.iter Integer.  Maximum number of iterations allowed in `qinf()`.
+#'   The default (`100`) is generous for most smooth CDFs.
+#'
+#' @return A **function** `q(p)` that takes a numeric probability (or vector of
+#'   probabilities) `p` with values in $[0,1]$ and returns the corresponding
+#'   quantile(s).  The output is named with the character representation of the
+#'   input probabilities.
+#'
+#' @details
+#' The heavy lifting is performed by `qinf()`, which solves for the root of
+#' \eqn{F(x) - p = 0}` using a bracketing method (e.g., bisection or Brent’s
+#' algorithm).  `make_quantile_function()` simply captures the user‑specified
+#' arguments and returns a thin wrapper that forwards them to `qinf()`.  This
+#' design makes it easy to generate quantile functions for custom CDFs without
+#' having to rewrite the root‑finding code each time.
+#'
+#' @examples
+#' ## -----------------------------------------------------------------
+#' ## Normal distribution (using the built‑in CDF)
+#' qnorm_custom <- make_quantile_function(cdf = pnorm,
+#'                                         lower = -10, upper = 10)
+#' qnorm_custom(c(0.025, 0.5, 0.975))
+#'
+#' ## Gamma distribution – custom CDF defined on the fly
+#' my_gamma_cdf <- function(x) pgamma(x, shape = 2, rate = 1)
+#' qgamma_custom <- make_quantile_function(cdf = my_gamma_cdf,
+#'                                          lower = 0, upper = 20)
+#' qgamma_custom(c(0.1, 0.5, 0.9))
+#'
+#' ## Verify against the base quantile function
+#' all.equal(qnorm_custom(0.8), qnorm(0.8))
+#'
+#' @importFrom stats pnorm pgamma
+#' @export
+make_quantile_function <- function(cdf, lower = -Inf, upper = Inf,
+                                   tol = .Machine$double.eps^0.5,
+                                   max.iter = 100) {
+  
+  quantile_function <- function(p) {
+  out <- qinf(F = cdf,
+       p = p,
+       lower = lower,
+       upper = upper,
+       tol = tol,
+       max.iter = max.iter)
+  names(out) <- as.character(p)
+  out
+  }
+  
+  return(quantile_function)
+}
+
diff --git a/R/estimators.R b/R/estimators.R
index 3e65ae1..5ab7222 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, fX=NULL) {
+make_matrix_creation <- function(seed, n, K, matrix_X, 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, fX=fX)
+    compute_matrix(seed=seed, a, n=n, K=K, matrix_X = matrix_X, Fv=Fv, guard=guard, fX=fX)
   }
 }
 
@@ -200,7 +200,7 @@ calculate_edge_density <- function(adj_matrix) {
   return(rho)
 }
 # test the estimator routines
-seed <- 121L #  121L this seed works exceptionally well
+seed <- 141 #  121L this seed works exceptionally well
 set.seed(seed)
 #X  <- matrix(seq(-1, 1, length.out = 5), ncol = 1)
 a  <- 20
@@ -225,9 +225,9 @@ adj <- compute_adj_matrix(
 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, fX=dnorm)
+Qa <- compute_matrix(seed, a=a, n=n, K=K, Fv=Fv, guard=guard, matrix_X=X)
+Qa2 <- compute_matrix(seed, a=a, n=n, K=K, 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, guard=guard, fX=dnorm)
 
 loss_func <- function(a) {
   Q_a <- calc_Q_a(a)
diff --git a/R/graphon_distribution.R b/R/graphon_distribution.R
index 12fe35e..849f627 100644
--- a/R/graphon_distribution.R
+++ b/R/graphon_distribution.R
@@ -95,7 +95,7 @@ pgraphon_analytical_1d <- function(y, a, fX, Fv,
 #'
 #' @export
 pgraphon <- function(
-    y,           # scalar
+    y,           # scalar or vector
     a,           # vector (coefficient vector, can be one-dimensional)
     Fv,          # CDF function of the v's (e.g., pnorm, pexp, etc.)
     X_matrix     # n x p matrix: each row is X_i
@@ -117,23 +117,112 @@ pgraphon <- function(
     stop("'Fv' must be a function (a CDF)")
   }
   ## 1.2 Compute the expected value ============================================
-  
+
   # Compute a^T X_i as in the paper. Here we switched to the R convention that each
   # observation is a row. The paper assumes that each observation is a vector
   inner_products <- as.vector(X_matrix %*% a)  # vector of length n
-  
+
   # outer(y, inner_prod, "-") creates an |y| × n matrix where element (j,i)
   # equals y[j] - inner_prod[i].
   dev_mat <- outer(y, inner_products, "-")
-  
+
   # Apply CDF Fv to each deviation
   cdf_vals <- Fv(dev_mat)
-  
+
   # Row means give the empirical expectation for each y_j
   out <- rowMeans(cdf_vals)
   out
 }
 
+#' Empirical CDF of a linear predictor plus a constant
+#'
+#' This helper computes the empirical cumulative distribution function (ECDF)
+#' of the quantity \eqn{z = a^\top X + v} and evaluates it at one or more
+#' points `y`.  It adds a few safety checks, supports vectorised `y`,
+#' and optionally returns the underlying ECDF function for later use.
+#'
+#' @param y Numeric vector of values at which the ECDF should be evaluated.
+#' @param a Numeric vector of coefficients (length must equal \code{ncol(X_matrix)}).
+#' @param v Either a numeric scalar or a numeric vector of length \code{nrow(X_matrix)}.
+#'          It is added to the linear predictor.
+#' @param X_matrix Numeric matrix (or something coercible to a matrix) that
+#'          contains the covariates. Rows correspond to observations.
+#' @param return_ecdf Logical (default \code{FALSE}).  If \code{TRUE} the
+#'          function returns a list with two components:
+#'          \itemize{
+#'            \item \code{prob} – the ECDF values at `y`.
+#'            \item \code$ecdf – the step‑function object returned by \code{ecdf()}.
+#'          }
+#'          If \code{FALSE} only the vector of probabilities is returned.
+#' @return Either a numeric vector of ECDF values (default) or a list
+#'         containing that vector and the ECDF function (if \code{return_ecdf=TRUE}).
+#' @examples
+#' ## Simple illustration with the built‑in mtcars data
+#' X <- as.matrix(mtcars[, c("wt", "hp")])
+#' a <- c(-0.5, 0.02)          # coefficients
+#' v <- 0.1                    # constant shift
+#' y_vals <- seq(-5, 5, length.out = 11)
+#' pgraphon_ecdf(y = y_vals, a = a, v = v, X_matrix = X)
+#'
+#' ## Get the ECDF function for later use
+#' out <- pgraphon_ecdf(y = y_vals, a = a, v = v, X_matrix = X,
+#'                      return_ecdf = TRUE)
+#' out$ecdf(0)   # probability that z <= 0
+#' @export
+pgraphon_ecdf <- function(y,
+                          a,
+                          v,
+                          X_matrix,
+                          return_ecdf = FALSE) {
+  ## ---- Input validation -------------------------------------------------
+  # y -------------------------------------------------
+  if (!is.numeric(y)) {
+    stop("`y` must be a numeric vector.")
+  }
+  
+  # X_matrix -------------------------------------------------
+  X_mat <- as.matrix(X_matrix)
+  if (ncol(X_mat) == 0L || nrow(X_mat) == 0L) {
+    stop("`X_matrix` must have at least one row and one column.")
+  }
+  
+  # a -------------------------------------------------
+  if (!is.numeric(a) || length(a) != ncol(X_mat)) {
+    stop("`a` must be a numeric vector with length equal to ncol(X_matrix).")
+  }
+  
+  # v -------------------------------------------------
+  if (!is.numeric(v)) {
+    stop("`v` must be numeric (scalar or vector).")
+  }
+  if (length(v) == 1L) {
+    v_vec <- rep(v, nrow(X_mat))
+  } else if (length(v) == nrow(X_mat)) {
+    v_vec <- v
+  } else {
+    stop("`v` must be either a scalar or a vector of length nrow(X_matrix).")
+  }
+  
+  ## ---- Compute linear predictor -----------------------------------------
+  # z = X %*% a + v
+  # Using as.numeric to coerce the 1‑column matrix result to a plain vector
+  z <- as.numeric(X_mat %*% a) + v_vec
+  
+  ## ---- Build ECDF -------------------------------------------------------
+  # ecdf() returns a step function that can be evaluated at any numeric vector
+  emp_cdf <- stats::ecdf(z)
+  
+  ## ---- Evaluate at y ----------------------------------------------------
+  prob <- emp_cdf(y)   # vectorised automatically
+  
+  ## ---- Return -----------------------------------------------------------
+  if (return_ecdf) {
+    list(prob = prob, ecdf = emp_cdf)
+  } else {
+    prob
+  }
+}
+
 # 2. Density Function ----------------------------------------------------------
 
 #' Empirical Graphon Density Estimate
diff --git a/R/singular_values.R b/R/singular_values.R
index 3503228..728dbbd 100644
--- a/R/singular_values.R
+++ b/R/singular_values.R
@@ -33,10 +33,6 @@ source(here::here("R", "graphon_distribution.R"))
 #'               generated.
 #' @param K      Positive integer.  Number of divisions of the unit interval;
 #'               the resulting grid has length `K+1`.
-#' @param fv     Density function of the latent variable \eqn{v}.  Must be
-#'               vectorised (i.e. accept a numeric vector and return a numeric
-#'               vector of the same length).  Typical examples are
-#'               `dnorm`, `dexp`, ….
 #' @param Fv     Cumulative distribution function of the latent variable
 #'               \eqn{v}.  Also has to be vectorised.  Typical examples are
 #'               `pnorm`, `pexp`, ….
@@ -112,7 +108,6 @@ compute_matrix <- function(
     a,
     n,
     K,
-    fv,
     Fv,
     fX=NULL,
     sample_X_fn=NULL,
@@ -132,6 +127,7 @@ compute_matrix <- function(
   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")
+  if (!is.logical(scaled)) stop("`scaled` must be a logical!")
   
   ## 1.2 Generate the Matrix X of covariates ===================================
   # If the argument matrix_X is present, use this matrix, otherwise generate one
@@ -150,11 +146,7 @@ compute_matrix <- function(
     stop("Number of columns of X (", ncol(X), ") must equal length(a) (", length(a), ")")
   }
   
-  ## 1.3 Create conditional density ============================================
-  # this is not used in the computation
-  # empir_cond_density <- create_cond_density(a, fv, Fv, X)
-  
-  ## 1.4 Compute the graphon quantiles =========================================
+  ## 1.3 Compute the graphon quantiles =========================================
   k <- seq(0, K) / K
   if (!is.null(guard)) {
     k[1] <- guard
@@ -164,7 +156,7 @@ compute_matrix <- function(
   # 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 ====================================================
+  ## 1.4 Build the matrix Q ====================================================
   inner_products = as.vector(X %*% a)
   
   # outer(y, x, "-")  gives a matrix with entry (j,i) = y[j] - x[i]