145 lines
5.1 KiB
R
145 lines
5.1 KiB
R
# Load required helper (the file must define `pgraphon`)
|
||
if (requireNamespace("here", quietly = TRUE)) {
|
||
source(here::here("R", "graphon_distribution.R"))
|
||
} else {
|
||
stop("Package `here` is required to locate the helper file.")
|
||
}
|
||
|
||
# TODO: Check the naming consistency across the documentation
|
||
|
||
#' Compute an adjacency matrix from a graphon model
|
||
#'
|
||
#' @title Generate a the adjacency matrix using a graphon specification
|
||
#' @description
|
||
#' Given a covariate matrix `X_matrix`, a vector of intercepts `v`, a
|
||
#' coefficient vector `a`, a symmetric graphon kernel `phi`, and a scaling
|
||
#' factor `rho_n`, this function draws a random undirected adjacency matrix.
|
||
#' The construction follows the steps:
|
||
#' 1. Compute latent variables \eqn{\xi_i = F_v(a^\top X_i + v_i)}.
|
||
#' 2. Generate a symmetric matrix of i.i.d. Uniform$[0,1]$ random numbers
|
||
#' with ones on the diagonal.
|
||
#' 3. Evaluate the graphon probabilities \eqn{p_{ij}= \rho_n \, \phi(\xi_i,\xi_j)}.
|
||
#' 4. Set \eqn{A_{ij}=1} if the uniform draw is smaller than the probability,
|
||
#' otherwise \eqn{0}. The resulting matrix is symmetric with zeros on the
|
||
#' diagonal.
|
||
#'
|
||
#' @param X_matrix A numeric matrix of dimension *n \times d* (covariates).
|
||
#' Must not contain `NA` or `NaN` values.
|
||
#' @param v A numeric vector of length *n* (intercepts).
|
||
#' Must be the same length as the number of rows of `X_matrix`.
|
||
#' @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]$.
|
||
#' @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`.
|
||
#' Must accept a numeric vector and return a numeric vector of the same length.
|
||
#'
|
||
#' @return An *n* × *n* integer matrix (`0L`/`1L`) representing the adjacency
|
||
#' matrix of an undirected graph.
|
||
#'
|
||
#' @examples
|
||
#' ## Simple example with a multiplicative graphon
|
||
#' set.seed(1)
|
||
#' X <- matrix(seq(-1, 1, length.out = 5), ncol = 1)
|
||
#' a <- 2
|
||
#' v <- seq(0, 0.8, length.out = 5)
|
||
#' phi_fun <- function(x, y) x * y # multiplicative kernel
|
||
#' adj <- compute_adj_matrix(
|
||
#' X_matrix = X,
|
||
#' v = v,
|
||
#' a = a,
|
||
#' phi = phi_fun,
|
||
#' rho_n = 0.5,
|
||
#' Fv = punif
|
||
#' )
|
||
#' adj
|
||
#'
|
||
#' @export
|
||
compute_adj_matrix <- function(
|
||
X_matrix,
|
||
v,
|
||
a,
|
||
phi,
|
||
rho_n,
|
||
Fv = pnorm
|
||
) {
|
||
## -------------------------- Guard rails -------------------------- ##
|
||
# Missing / NULL checks
|
||
if (missing(X_matrix) || is.null(X_matrix)) {
|
||
stop("`X_matrix` is missing or NULL.")
|
||
}
|
||
if (missing(v) || is.null(v)) {
|
||
stop("`v` is missing or NULL.")
|
||
}
|
||
if (missing(a) || is.null(a)) {
|
||
stop("`a` is missing or NULL.")
|
||
}
|
||
if (missing(phi) || is.null(phi)) {
|
||
stop("`phi` is missing or NULL.")
|
||
}
|
||
if (missing(rho_n) || is.null(rho_n)) {
|
||
stop("`rho_n` is missing or NULL.")
|
||
}
|
||
if (missing(Fv) || is.null(Fv)) {
|
||
stop("`Fv` is missing or NULL.")
|
||
}
|
||
|
||
# Type / dimension checks
|
||
if (!is.matrix(X_matrix) || !is.numeric(X_matrix)) {
|
||
stop("`X_matrix` must be a numeric matrix.")
|
||
}
|
||
if (any(is.na(X_matrix) | is.nan(X_matrix))) {
|
||
stop("`X_matrix` must not contain NA or NaN values.")
|
||
}
|
||
n <- nrow(X_matrix)
|
||
d <- ncol(X_matrix)
|
||
|
||
if (!is.numeric(v) || length(v) != n) {
|
||
stop("`v` must be a numeric vector of length equal to nrow(`X_matrix`).")
|
||
}
|
||
if (!is.numeric(a) || length(a) != d) {
|
||
stop("`a` must be a numeric vector of length equal to ncol(`X_matrix`).")
|
||
}
|
||
if (!is.function(phi)) {
|
||
stop("`phi` must be a function of two numeric arguments.")
|
||
}
|
||
if (!is.numeric(rho_n) || length(rho_n) != 1 || rho_n < 0 || rho_n > 1) {
|
||
stop("`rho_n` must be a numeric scalar in [0, 1].")
|
||
}
|
||
if (!is.function(Fv)) {
|
||
stop("`Fv` must be a function (CDF) that accepts a numeric vector.")
|
||
}
|
||
|
||
## -------------------------- Core computation -------------------------- ##
|
||
# 1. Latent variables ξ_i = Fv(aᵀ X_i + v_i)
|
||
lin_pred <- as.vector(X_matrix %*% a) + v
|
||
xi <- pgraphon(y = lin_pred, a = a, Fv = Fv, X_matrix = X_matrix)
|
||
|
||
# 2. Symmetric matrix of Uniform[0,1] draws with 1 on the diagonal
|
||
# (only the upper triangle is sampled, then mirrored)
|
||
rand_mat <- matrix(0, n, n)
|
||
upper_idx <- which(upper.tri(rand_mat, diag = FALSE), arr.ind = TRUE)
|
||
rand_mat[upper_idx] <- runif(n * (n - 1) / 2)
|
||
rand_mat <- rand_mat + t(rand_mat) + diag(1, n)
|
||
|
||
# 3. Graphon probabilities p_ij = rho_n * phi(ξ_i, ξ_j)
|
||
graphon_fun <- function(x, y) rho_n * phi(x, y)
|
||
probabilities <- outer(xi, xi, FUN = graphon_fun)
|
||
|
||
# 4. Adjacency matrix: A_ij = 1 if U_ij < p_ij, else 0
|
||
adj_mat <- (rand_mat < probabilities) * 1L
|
||
|
||
## Return the adjacency matrix
|
||
adj_mat
|
||
}
|
||
|
||
set.seed(1)
|
||
X <- matrix(c(-1, -0.5, 0, 0.5, 1), nrow = 5, ncol=1)
|
||
a <- 2.0
|
||
v <- c(0, 0.2, 0.4, 0.6, 0.8)
|
||
adj <- compute_adj_matrix(X, v, a, phi = function(x,y) {x * y}, 0.5, Fv = punif)
|
||
adj
|
||
|