Add scaled option for matrix generation

R/singular_value_plot.R

@@ -0,0 +1,234 @@
+# Load function ----------------------------------------------------------------
+source(here::here("R", "singular_values.R"))
+source(here::here("R", "graphon_distribution.R"))
+
+
+# Convert a call or character to a nicely formatted character string.
+#   * If the user supplied a character, we keep it unchanged.
+#   * If the user supplied a call (e.g. quote(20 / sqrt(x))) we deparse it
+#     and collapse the result to a single line.
+expr_to_label <- function(expr) {
+  if (is.character(expr)) {
+    expr
+  } else {
+    # deparse returns a character vector; collapse to one line
+    paste(deparse(expr), collapse = "")
+  }
+}
+
+
+#' Compute the smallest singular value of a sequence of matrices Q(K)
+#'
+#' @title Smallest singular values for a family of matrices Q(K)
+#' @description
+#'   For a given vector of coefficients `a`, sample size `n` and a maximum
+#'   rank `maxK`, this function repeatedly calls `compute_matrix()` to build a
+#'   matrix `Q` for each rank `K = 1, …, maxK`.  The smallest singular value of
+#'   each `Q` is extracted with `compute_minmax_sv()`.  The result can be
+#'   returned as a numeric vector and, if desired, plotted together with a
+#'   user‑supplied reference curve.
+#'
+#' @param a Numeric vector of coefficients that are passed to `compute_matrix()`.
+#' @param n Positive integer, the sample size used inside the sampling function.
+#' @param maxK Positive integer, the largest rank `K` for which a matrix `Q`
+#'   will be built.  The function will evaluate `K = 1:maxK`.
+#' @param seed Integer (default = 1) used as the random‑seed argument for
+#'   `compute_matrix()`.  Supplying a seed makes the whole procedure
+#'   reproducible.
+#' @param sampler_fn Function that draws a sample of size `n`.  It must accept a
+#'   single argument `n` and return a **numeric matrix** with `n` rows and
+#'   one column (the shape used in the original script).  The default is a
+#'   thin wrapper around `rnorm()`.
+#' @param fv Function giving the density of the underlying distribution
+#'   (default = `dnorm`).  It is passed unchanged to `compute_matrix()`.
+#' @param Fv Function giving the cumulative distribution function
+#'   (default = `pnorm`).  It is passed unchanged to `compute_matrix()`.
+#' @param guard Positive numeric guard used inside `compute_matrix()` to avoid
+#'   division‑by‑zero or log‑of‑zero problems (default = `1e-12`).
+#' @param plot Logical, whether to produce a quick base‑R plot of the
+#'   smallest singular values (`TRUE` by default).  If `FALSE` the function
+#'   only returns the numeric vector.
+#' @param add_curve Logical, whether to overlay a reference curve on the plot.
+#'   Ignored when `plot = FALSE`.  Default = `TRUE`.
+#' @param curve_expr Expression (as a *character* or *call*) that defines the
+#'   reference curve.  The default reproduces the line you used
+#'   `20 / sqrt(x)`.  The expression must be a valid R expression in which the
+#'   variable `x` stands for the horizontal axis.
+#' @param curve_from,curve_to Numeric limits for the reference curve.  By
+#'   default they are set to the range of `K` (`1:maxK`).
+#' @param curve_col Colour of the reference curve (default = `"red"`).
+#' @param curve_lwd Line width of the reference curve (default = 2).
+#' @param log_plot If True, then the y-axis is on a log scale.
+#' @param main_title Main title for the plot
+#' @return A list with the following components  
+#'   \item{K}{Integer vector `1:maxK`.}  
+#'   \item{sv}{Numeric vector of the smallest singular values for each `K`.}  
+#'   \item{plot}{If `plot = TRUE`, the value returned by `graphics::plot()`
+#'   (normally `NULL`).  If `plot = FALSE` this element is omitted.}
+#' @examples
+#' ## Run the whole routine with the defaults
+#' res <- smallest_sv_sequence(
+#'   a = c(0.4), n = 400, maxK = 20,
+#'   seed = 3, guard = 1e-12
+#' )
+#' head(res$sv)
+#'
+#' ## Supply a custom sampler (e.g., uniform)
+#' my_sampler <- function(m) matrix(runif(m, -1, 1), ncol = 1)
+#' res2 <- smallest_sv_sequence(
+#'   a = c(0.4), n = 400, maxK = 10,
+#'   sampler_fn = my_sampler,
+#'   plot = FALSE
+#' )
+#' plot(res2$K, res2$sv, type = "b")
+#' @importFrom graphics plot lines curve
+#' @export
+smallest_sv_sequence <- function(
+    a,
+    n,
+    maxK,
+    seed = 1L,
+    sampler_fn = function(m) matrix(rnorm(m), ncol = 1L),
+    fv = dnorm,
+    Fv = pnorm,
+    guard = 1e-12,
+    plot = TRUE,
+    add_curve = TRUE,
+    curve_expr = quote(20 / sqrt(x)),
+    curve_from = NULL,
+    curve_to = NULL,
+    curve_col = "red",
+    curve_lwd = 2, 
+    log_plot = FALSE,
+    main_title = "Smallest singular value vs. K"
+) {
+  ## 1. Input validation =======================================================
+  if (!is.numeric(a) || length(a) == 0) {
+    stop("`a` must be a non‑empty numeric vector.")
+  }
+  if (!is.numeric(n) || length(n) != 1L || n <= 0 || n != as.integer(n)) {
+    stop("`n` must be a positive integer.")
+  }
+  if (!is.numeric(maxK) || length(maxK) != 1L || maxK <= 0 ||
+      maxK != as.integer(maxK)) {
+    stop("`maxK` must be a positive integer.")
+  }
+  if (!is.function(sampler_fn)) {
+    stop("`sampler_fn` must be a function that takes a single integer argument `n`.")
+  }
+  if (!is.function(fv) || !is.function(Fv)) {
+    stop("`fv` and `Fv` must be corresponding density functions  and cdf (e.g., dnorm/ pnorm).")
+  }
+  if (!is.numeric(guard) || length(guard) != 1L || guard <= 0) {
+    stop("`guard` must be a single positive numeric value.")
+  }
+  if (!is.logical(plot) || length(plot) != 1L) {
+    stop("`plot` must be a single logical value.")
+  }
+  if (!is.logical(add_curve) || length(add_curve) != 1L) {
+    stop("`add_curve` must be a single logical value.")
+  }
+  if (!inherits(curve_expr, "call") && !is.character(curve_expr)) {
+    stop("`curve_expr` must be a call (e.g., quote(20/sqrt(x))) or a character string.")
+  }
+  if (!is.character(main_title)){
+    stop("`main_title` must be a character vector.")
+  }
+  
+  ## 2. Prepare storage ========================================================
+  K_vec <- seq_len(maxK)
+  smallest_sv <- numeric(maxK)
+  
+  ## 3. Main loop – build Q(K) and extract the smallest singular value =========
+  for (K in K_vec) {
+    Q <- compute_matrix(
+      seed = seed,
+      a = a,
+      n = n,
+      K = K,
+      sample_X_fn = sampler_fn,
+      fv = fv,
+      Fv = Fv,
+      guard = guard
+    )
+    
+    Q <- 1 /sqrt(n) * Q
+    
+    sv_res <- compute_minmax_sv(Q)
+    if (!is.list(sv_res) || is.null(sv_res$smallest_singular_value)) {
+      stop("`compute_minmax_sv()` must return a list containing `$smallest_singular_value`.")
+    }
+    smallest_sv[K] <- sv_res$smallest_singular_value
+  }
+  
+  ## 4. Plotting (optional) ====================================================
+  if (plot) {
+    ## Basic scatter/line plot of the singular values
+    par(mar = c(5, 4, 4, 8))  # extra space on the right for the legend
+    plot_args <- list(
+      x = K_vec,
+      y = smallest_sv,
+      type = "b",
+      pch = 19,
+      col = "steelblue",
+      xlab = "K subdivisions",
+      ylab = "Smallest singular value of Q",
+      main = main_title
+    )
+    if (log_plot) plot_args$log <- "y"
+    
+    do.call(graphics::plot, plot_args)
+    # add legend. The par(xpd = ...) allows drawing outside of the plot region.
+    par(xpd = TRUE)
+    legend("topright",
+           inset=c(-0.2,0), 
+           legend=c("SV of Q"), 
+           col="steelblue", 
+           title="Legend",
+           pch = 16,
+           bty = "n")
+    par(xpd = FALSE)
+    ## Add the reference curve if requested
+    if (add_curve) {
+      ## Determine sensible defaults for the curve limits
+      if (is.null(curve_from)) curve_from <- min(K_vec)
+      if (is.null(curve_to))   curve_to   <- max(K_vec)
+      
+      ## `curve()` expects an *expression* where `x` is the variable.
+      ## If the user supplied a character string we turn it into a call.
+      if (is.character(curve_expr)) curve_expr <- parse(text = curve_expr)[[1L]]
+      curve_fun <- function(x) eval(curve_expr)
+      graphics::curve(
+        expr = curve_fun,
+        from = curve_from,
+        to   = curve_to,
+        add  = TRUE,
+        col  = curve_col,
+        lwd  = curve_lwd
+      )
+      
+      # add label with the curve expression
+      label_txt <- expr_to_label(curve_expr)
+      x_pos <- curve_from + 0.8 * (curve_to - curve_from)
+      y_pos <- 0.85 * max(smallest_sv)
+      graphics::text(
+        x = x_pos, y = y_pos,
+        labels = label_txt,
+        col    = curve_col,
+        pos    = 4,          # 4 = right‑justified (so the text sits left of the point)
+        offset = 0.5,        # a small gap between the point and the text
+        cex    = 0.9        # slightly smaller than default text size
+      )
+      
+    }
+  }
+  ## 5. Return a tidy list =====================================================
+  out <- list(
+    K = K_vec,
+    sv = smallest_sv
+  )
+  # if (plot) out$plot <- NULL   ## the plot is already drawn; we return NULL for consistency
+  out
+}
+
+

R/singular_values.R

@@ -48,6 +48,8 @@ source("R/graphon_distribution.R")
 #'   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.
+#' @param scaled Rescales the matrix according to the dimension with the factor
+#' 1 / sqrt(n). Default value is FALSE.
 #'
 #' @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
@@ -105,7 +107,8 @@ compute_matrix <- function(
     sample_X_fn,
     fv,
     Fv,
-    guard = sqrt(.Machine$double.eps)
+    guard = sqrt(.Machine$double.eps),
+    scaled = FALSE
 ) {
   ## 1.1 Check inputs ==========================================================
   if (!is.numeric(seed) || length(seed) != 1) stop("'seed' must be a single number")
@@ -147,6 +150,7 @@ compute_matrix <- function(
   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
 }
 
@@ -193,6 +197,9 @@ compute_matrix <- function(
 compute_minmax_sv <- function(M) {
   s <- svd(M, nu=0, nv=0)$d
   
+  # just a check if we compute the right thing
+  # s <- sqrt(eigen(M %*% t(M), symmetric = TRUE, only.value=TRUE)$values)
+
   list(
     largest_singular_value  = max(s),
     smallest_singular_value = min(s) # smallest non zero singular value

scripts/plot_densities.R

@@ -50,14 +50,14 @@ title(
 )
 
 ## 3.3 Plot one dimensional X_i ================================================
-p3 <- create_cond_density(a1, dnorm, pnorm, X1)
-givenX <- c(-0.5)
+p3 <- create_cond_density(2, dnorm, pnorm, X1)
+givenX <- c(-3.0)
 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 = ", "))) ~ "," ~
+    a == (.(paste(c(2), collapse = ", "))) ~ "," ~
       N == .(N) ~ "," ~
       x == (.(paste(givenX, collapse = ", ")))
   ),

scripts/plot_sv.R

@@ -0,0 +1,364 @@
+# Load function ----------------------------------------------------------------
+source(here::here("R", "singular_values.R"))
+source(here::here("R", "graphon_distribution.R"))
+source(here::here("R", "singular_value_plot.R"))
+
+# https://stackoverflow.com/a/5790430
+resetPar <- function() {
+  dev.new()
+  op <- par(no.readonly = TRUE)
+  dev.off()
+  op
+}
+
+calc_conv_rate <- function(x,y) {
+  if (!is.numeric(x) || length(x) == 0) {
+    stop("`x` must be a non‑empty numeric vector.")
+  }
+  if (!is.numeric(y) || length(y) == 0) {
+    stop("`y` must be a non‑empty numeric vector.")
+  }
+  if (length(x) != length(y)) {
+    stop("`x` and `y` must have the same length.")
+  }
+  
+  df <- data.frame("x" = x, "y" = y)
+  lm_model <- lm(log(y) ~ log(x), data=df)
+  C <- exp(coefficients(lm_model)[[1]])
+  alpha <- coefficients(lm_model)[[2]]
+  
+  df[, "y_pred"] <- C * df[, "x"]^alpha
+  df[, "residual"] <- df[, "y"] - df[, "y_pred"]
+  
+  out <- list(
+    "C" = C,
+    "alpha" = alpha,
+    "obs" = df
+  )
+  out
+}
+
+calc_exp_conv_rate <- function(x,y) {
+  if (!is.numeric(x) || length(x) == 0) {
+    stop("`x` must be a non‑empty numeric vector.")
+  }
+  if (!is.numeric(y) || length(y) == 0) {
+    stop("`y` must be a non‑empty numeric vector.")
+  }
+  if (length(x) != length(y)) {
+    stop("`x` and `y` must have the same length.")
+  }
+  
+  df <- data.frame("x" = x, "y" = y)
+  fit <- nls(y ~ C * exp(r * x^m),
+             start = list(C = min(y), r= -0.5, m = 1))
+}
+
+# Nearly match with sample function --------------------------------------------
+# v ~ N(0,1) and X ~ discrete Uniform on [1:n]
+out <- smallest_sv_sequence(
+  a = c(0.5),
+  n = 9,
+  maxK= 3,
+  sampler_fn = sample,
+  guard=1e-12,
+  plot=TRUE,
+  log_plot = TRUE,
+  curve_expr = quote(1 / x^0.545)
+)
+conv_rate <- calc_conv_rate(out$K, out$sv)
+
+# Normally distributed X ~ N(0,1) and v ~ N(0,1) -------------------------------
+out <- smallest_sv_sequence(
+  a = c(0.5),
+  n = 1200,
+  maxK = 20,
+  sampler_fn =function(n) matrix(rnorm(n), ncol = 1L),
+  guard=1e-12,
+  plot=TRUE,
+  log_plot = TRUE,
+  curve_expr = quote(1.5 * exp(-0.95 * x^1.34))
+  #curve_expr = quote( 1/exp(x^1.32))
+)
+# convergence rate does not work here, probably because the underlying model
+# does not work well
+conv_rate <- calc_conv_rate(out$K[1:20], out$sv[1:20])
+
+# Uniform distributed X ~ U[0,1] and v ~ N(0,1) --------------------------------
+out <- smallest_sv_sequence(
+  a = c(0.5),
+  n = 400,
+  maxK = 20,
+  sampler_fn =function(n) matrix(runif(n), ncol = 1L),
+  guard=1e-12,
+  plot=TRUE,
+  log_plot = TRUE,
+  curve_expr = quote(1* exp(-1.1 * x^1.5))
+)
+# here the optimal fit does not work too, probably other model
+calc_conv_rate(out$K[1:9], out$sv[1:9])
+
+# Compare of parameters of Normal distribution ----------------------------------
+#
+out_sd0_5 <- smallest_sv_sequence(
+  a = c(0.5),
+  n = 400,
+  maxK = 20,
+  sampler_fn =function(n) matrix(rnorm(n), ncol = 1L),
+  guard=1e-12,
+  plot=TRUE,
+  log_plot = TRUE,
+  fv = function(x) {dnorm(x, mean=0, sd=0.5)},
+  Fv = function(x) {pnorm(x, mean=0, sd=0.5)},
+  main_title="Smallest SV of v~ N(0,0.5^2) distribution"
+)
+
+out_sd1 <- smallest_sv_sequence(
+  a = c(0.5),
+  n = 400,
+  maxK = 20,
+  sampler_fn =function(n) matrix(rnorm(n), ncol = 1L),
+  guard=1e-12,
+  plot=TRUE,
+  log_plot = TRUE,
+  fv = function(x) {dnorm(x, mean=0, sd=1)},
+  Fv = function(x) {pnorm(x, mean=0, sd=1)},
+  main_title="Smallest SV of v~ N(0,1) distribution"
+)
+
+out_sd2 <- smallest_sv_sequence(
+  a = c(0.5),
+  n = 400,
+  maxK = 20,
+  sampler_fn =function(n) matrix(rnorm(n), ncol = 1L),
+  guard=1e-12,
+  plot=TRUE,
+  log_plot = TRUE,
+  fv = function(x) {dnorm(x, mean=0, sd=2)},
+  Fv = function(x) {pnorm(x, mean=0, sd=2)},
+  main_title="Smallest SV of v~ N(0,2^2) distribution"
+)
+
+out_sd4 <- smallest_sv_sequence(
+  a = c(0.5),
+  n = 400,
+  maxK = 20,
+  sampler_fn =function(n) matrix(rnorm(n), ncol = 1L),
+  guard=1e-12,
+  plot=TRUE,
+  log_plot = TRUE,
+  fv = function(x) {dnorm(x, mean=0, sd=4)},
+  Fv = function(x) {pnorm(x, mean=0, sd=4)},
+  main_title="Smallest SV of v~ N(0,4^2) distribution"
+)
+
+par(mar = c(5, 4, 4, 8)) 
+plot(out_sd0_5$K, out_sd0_5$sv,
+     type = "b",
+     pch = 19,
+     col = "#D3BA68FF",
+     xlab = "K subdivisions",
+     ylab = "Smallest singular value of Q",
+     main="smallest SV for different variances of a normal distribution",
+     sub = "n = 400, a = 0.5",
+     log="y")
+lines(out_sd1$K, out_sd1$sv,
+     type="b", pch=19, col="#D5695DFF", add=TRUE)
+lines(out_sd2$K, out_sd2$sv,
+      type = "b", pch=19, col="#5D8CA8FF", add=TRUE)
+lines(out_sd4$K, out_sd4$sv,
+      type = "b", pch=19, col="#65A479FF", add=TRUE)
+par(xpd = TRUE)
+legend("topright",
+       inset=c(-0.2,0), 
+       legend=c("sd=0.5", "sd=1", "sd=2", "sd=4"), 
+       col=c("#D3BA68FF", "#D5695DFF","#5D8CA8FF", "#65A479FF" ), 
+       title="Legend",
+       pch = 16,
+       bty = "n")
+
+# Break phenomena of Exp-distribution ------------------------------------------
+out_exp <- smallest_sv_sequence(
+  a = c(0.5),
+  n = 400,
+  maxK = 20,
+  sampler_fn =function(n) matrix(rnorm(n), ncol = 1L),
+  guard=1e-12,
+  plot=TRUE,
+  log_plot = TRUE,
+  fv = function(x) {dexp(x, rate=1)},
+  Fv = function(x) {pexp(x, rate=1)},
+  main_title="Smallest SV of v~ Exp(1) distribution"
+)
+
+par(resetPar)
+plot(out_exp$K, out_exp$sv,
+     log="y",
+     xlab="K subdivsions",
+     ylab="Smallest singular value of Q",
+     col="steelblue",
+     type="b",
+     main="Smallest singular value for v ~ Exp(1)",
+     sub="a = 0.5, n = 400")
+arrows(8, 1e-5, 6.5, 1e-7, angle=20, lty = 1, lwd=2)
+text(8.5, 1e-5, "Break only seen for exp-distribution", pos=4)
+
+## Observations of the break point depending the rate --------------------------
+#
+# Note: this also depends on the the parameter n of samples
+out_exp_1 <- smallest_sv_sequence(
+  a = c(0.5),
+  n = 400,
+  maxK = 80,
+  sampler_fn =function(n) matrix(rnorm(n), ncol = 1L),
+  guard=1e-12,
+  plot=FALSE,
+  fv = function(x) {dexp(x, rate=1)},
+  Fv = function(x) {pexp(x, rate=1)}
+)
+
+out_exp_2 <- smallest_sv_sequence(
+  a = c(0.5),
+  n = 400,
+  maxK = 80,
+  sampler_fn =function(n) matrix(rnorm(n), ncol = 1L),
+  guard=1e-12,
+  plot=FALSE,
+  fv = function(x) {dexp(x, rate=2)},
+  Fv = function(x) {pexp(x, rate=2)}
+)
+
+out_exp_3 <- smallest_sv_sequence(
+  a = c(0.5),
+  n = 400,
+  maxK = 80,
+  sampler_fn =function(n) matrix(rnorm(n), ncol = 1L),
+  guard=1e-12,
+  plot=FALSE,
+  fv = function(x) {dexp(x, rate=3)},
+  Fv = function(x) {pexp(x, rate=3)}
+)
+
+out_exp_4 <- smallest_sv_sequence(
+  a = c(0.5),
+  n = 400,
+  maxK = 80,
+  sampler_fn =function(n) matrix(rnorm(n), ncol = 1L),
+  guard=1e-12,
+  plot=FALSE,
+  fv = function(x) {dexp(x, rate=4)},
+  Fv = function(x) {pexp(x, rate=4)}
+)
+
+par(mar = c(5, 4, 4, 8)) 
+plot(out_exp_1$K, out_exp_1$sv,
+     type="b", col="#D3BA68FF", log="y",
+     main="Smallest SV of Q for different rates of Exp-distribution",
+     ylab="Smallest singular value of Q",
+     xlab="K subdivisions",
+     sub="a = 0.5, n = 400, depending also on n")
+lines(out_exp_2$K, out_exp_2$sv, type="b", col="#D5695DFF")
+lines(out_exp_3$K, out_exp_3$sv, type="b", col="#5D8CA8FF")
+lines(out_exp_4$K, out_exp_4$sv, type="b", col="#65A479FF")
+par(xpd=TRUE)
+legend("topright",
+       inset=c(-0.2,0), 
+       legend=c(expression(lambda == 1), expression(lambda == 2), expression(lambda == 3), expression(lambda == 4)),
+       col=c("#D3BA68FF", "#D5695DFF","#5D8CA8FF", "#65A479FF" ), 
+       title="Rate",
+       pch = 16,
+       bty = "n")
+
+
+# Use of the gamma distribution ------------------------------------------------
+# Fix the parameters, such that the mean stays the same and the variance is
+# changing.
+# From the documentation
+# Note that for smallish values of shape (and moderate scale) a large parts of 
+# the mass of the Gamma distribution is on values of x
+# so near zero that they will be represented as zero in computer arithmetic.
+# So rgamma may well return values which will be represented as zero. 
+# (This will also happen for very large values of scale since the actual generation is done for scale = 1.)
+#
+# Take E(X) = 5, so sigma = 5 / alpha, and with this we have
+# Var(X) = sigma^2 * alpha = 25 / alpha. 
+# -> Increasing alpha yields lower variance
+
+# alpha = 1
+alpha <- 1.0
+out_gamma_1 <- smallest_sv_sequence(
+  a = c(0.5),
+  n = 400,
+  maxK = 20,
+  sampler_fn =function(n) matrix(rnorm(n), ncol = 1L),
+  guard=1e-12,
+  plot=TRUE,
+  log_plot = TRUE,
+  fv = function(x) {dgamma(x, shape = alpha, rate = 5 / alpha)},
+  Fv = function(x) {pexp(x, rate=1)},
+  main_title="Smallest SV of v~ Gamma(1, 5) distribution"
+)
+
+# alpha = 2
+alpha <- 2.0
+out_gamma_2 <- smallest_sv_sequence(
+  a = c(0.5),
+  n = 400,
+  maxK = 20,
+  sampler_fn =function(n) matrix(rnorm(n), ncol = 1L),
+  guard=1e-12,
+  plot=TRUE,
+  log_plot = TRUE,
+  fv = function(x) {dgamma(x, shape = alpha, rate = 5 / alpha)},
+  Fv = function(x) {pexp(x, rate=1)},
+  main_title="Smallest SV of v~ Gamma(2, 2.5) distribution"
+)
+
+# alpha = 3
+alpha <- 3.0
+out_gamma_3 <- smallest_sv_sequence(
+  a = c(0.5),
+  n = 400,
+  maxK = 20,
+  sampler_fn =function(n) matrix(rnorm(n), ncol = 1L),
+  guard=1e-12,
+  plot=TRUE,
+  log_plot = TRUE,
+  fv = function(x) {dgamma(x, shape = alpha, rate = 5 / alpha)},
+  Fv = function(x) {pexp(x, rate=1)},
+  main_title="Smallest SV of v~ Gamma(3, 5/3) distribution"
+)
+
+# alpha = 4
+alpha <- 4.0
+out_gamma_4 <- smallest_sv_sequence(
+  a = c(0.5),
+  n = 400,
+  maxK = 20,
+  sampler_fn =function(n) matrix(rnorm(n), ncol = 1L),
+  guard=1e-12,
+  plot=TRUE,
+  log_plot = TRUE,
+  fv = function(x) {dgamma(x, shape = alpha, rate = 5 / alpha)},
+  Fv = function(x) {pexp(x, rate=1)},
+  main_title="Smallest SV of v~ Gamma(3, 5/3) distribution"
+)
+
+par(mar = c(5, 4, 4, 8)) 
+plot(out_gamma_1$K, out_gamma_1$sv,
+     type="b", col="#D3BA68FF", log="y",
+     main="Smallest SV of Q for variance of the Gamma distribution",
+     ylab="Smallest singular value of Q",
+     xlab="K subdivisions",
+     sub="a = 0.5, n = 400")
+lines(out_gamma_2$K, out_gamma_2$sv, type="b", col="#D5695DFF")
+lines(out_gamma_3$K, out_gamma_3$sv, type="b", col="#5D8CA8FF")
+lines(out_gamma_4$K, out_gamma_4$sv, type="b", col="#65A479FF")
+par(xpd=TRUE)
+legend("topright",
+       inset=c(-0.2,0), 
+       legend=c(expression(alpha == 1), expression(alpha == 2), expression(alpha == 3), expression(alpha == 4)),
+       col=c("#D3BA68FF", "#D5695DFF","#5D8CA8FF", "#65A479FF" ), 
+       title="Rate",
+       pch = 16,
+       bty = "n")