🧠 Injected Memory Context
The debug output confirms exactly what was diagnosed in memory entry [PshfUZWW]: After the first iteration, a new point x=0.1167288 is inserted, and the z values become [0, -0.67743, 0.740992, 1.518016, Inf] — z[2] = -0.67743 is LESS than z[1] = 0, making the z array unsorted.
The root cause identified in memory: when all dh values are nearly -1 (exponential distribution), the intersection formula produces z values outside the support and not monotonically increasing. The midpoint fix with threshold 1e-8 is NOT being triggered because the dh differences are small but larger than 1e-8 (floating point differences in numerical derivatives).
The memory entry [PshfUZWW] recommends this specific fix for construct_envelope:
1. **Clamp all interior z[i] to [support[1], support[2]]** — this would clamp -0.67743 to 0
2. **Ensure z[i] >= z[i-1] for all i** (enforce monotonicity)
3. OR: when dh values are nearly equal, ensure z stays within [x[i-1], x[i]]
The construct_envelope function needs the `support` parameter to be used for clamping, and after computing all z values, enforce monotonicity. This would fix both TEST 2 (exponential) and TEST 7 (bimodal), since the bimodal test also hits the same unsorted z issue before log-concavity can be detected.
Additionally, the `findInterval(x, env$z)` call in `eval_upper_envelope` crashes when z is unsorted — this is the actual error path now visible in the debug output.
Message
Analysis: Excellent! The debug output reveals the actual problem: the z values become unsorted! After adding a new point at x=0.1167288, the z array becomes [0, -0.67743, 0.740992, 1.518016, Inf], where z[2] = -0.67743 < z[1] = 0. This violates the requirement that z must be sorted, causing `findInterval` to fail.
The issue is that when all derivatives are nearly -1 (exponential), the intersection formula produces z values that can be outside the valid range and not monotonically increasing. I need to:
1. Clamp z values to the support range
2. Ensure z values are monotonically increasing
Plan: Fix the `construct_envelope` function to:
1. Clamp all interior z values to be within [support[1], support[2]]
2. Ensure monotonicity: z[i] >= z[i-1] for all i
3. This will fix both TEST 2 and TEST 7
{
"keystrokes": "cat > /app/ars.R << 'EOF'\n################################################################################\n# Adaptive Rejection Sampler (ARS)\n# Based on Gilks & Wild (1992)\n# \n# Implementation of adaptive rejection sampling for log-concave densities\n################################################################################\n\n#' Adaptive Rejection Sampler\n#'\n#' @param n Number of samples to generate\n#' @param f Function that computes the (possibly unnormalized) density (vectorized)\n#' @param support Vector of length 2 giving the support of the distribution\n#' @param start_points Initial abscissae for constructing envelopes (optional)\n#' @param ... Additional arguments passed to f\n#'\n#' @return Vector of n samples from the target distribution\n#' @export\nars <- function(n, f, support = c(-Inf, Inf), start_points = NULL, ...) {\n \n # Input validation\n if (!is.numeric(n) || length(n) != 1 || n <= 0 || n != floor(n)) {\n stop(\"n must be a positive integer\")\n }\n \n if (!is.function(f)) {\n stop(\"f must be a function\")\n }\n \n if (!is.numeric(support) || length(support) != 2 || support[1] >= support[2]) {\n stop(\"support must be a numeric vector of length 2 with support[1] < support[2]\")\n }\n \n # Initialize abscissae (starting points)\n if (is.null(start_points)) {\n start_points <- initialize_abscissae(support)\n } else {\n if (!is.numeric(start_points) || length(start_points) < 2) {\n stop(\"start_points must be a numeric vector with at least 2 elements\")\n }\n if (any(start_points <= support[1]) || any(start_points >= support[2])) {\n stop(\"start_points must be within the support\")\n }\n start_points <- sort(start_points)\n }\n \n # Evaluate log-density and derivatives at starting points\n h_vals <- log(f(start_points, ...))\n if (any(!is.finite(h_vals))) {\n stop(\"Density function returns non-positive or non-finite values at starting points\")\n }\n \n dh_vals <- compute_derivatives(start_points, f, ...)\n \n # Check log-concavity at starting points\n if (!check_log_concavity(start_points, h_vals, dh_vals)) {\n stop(\"Function does not appear to be log-concave at starting points\")\n }\n \n # Store abscissae and their properties\n abscissae <- list(\n x = start_points,\n h = h_vals,\n dh = dh_vals\n )\n \n # Generate samples\n samples <- numeric(n)\n n_generated <- 0\n max_iter <- n * 1000 # Safety limit\n iter <- 0\n \n while (n_generated < n && iter < max_iter) {\n iter <- iter + 1\n \n # Update envelope\n env <- construct_envelope(abscissae, support)\n \n # Sample from upper envelope\n x_star <- sample_upper_envelope(env, support)\n \n # Evaluate densities\n h_star <- log(f(x_star, ...))\n if (!is.finite(h_star)) {\n warning(\"Sampled point has non-finite log-density, skipping\")\n next\n }\n \n u_star <- eval_upper_envelope(x_star, env)\n l_star <- eval_lower_envelope(x_star, abscissae)\n \n # Acceptance test\n w <- runif(1)\n \n # Squeeze test (accept)\n if (log(w) <= l_star - u_star) {\n n_generated <- n_generated + 1\n samples[n_generated] <- x_star\n } else {\n # Rejection test\n if (log(w) <= h_star - u_star) {\n n_generated <- n_generated + 1\n samples[n_generated] <- x_star\n }\n \n # Update abscissae (add rejected point to improve envelope)\n dh_star <- compute_derivatives(x_star, f, ...)\n \n # Check log-concavity\n if (!check_log_concavity_update(abscissae, x_star, h_star, dh_star)) {\n stop(\"Function is not log-concave: detected during sampling at x = \", \n round(x_star, 4))\n }\n \n # Insert new point\n abscissae <- insert_abscissa(abscissae, x_star, h_star, dh_star)\n }\n }\n \n if (iter >= max_iter) {\n warning(\"Maximum iterations reached, returning \", n_generated, \" samples\")\n samples <- samples[1:n_generated]\n }\n \n return(samples)\n}\n\n#' Initialize starting abscissae\ninitialize_abscissae <- function(support) {\n if (is.finite(support[1]) && is.finite(support[2])) {\n # Bounded support\n range <- support[2] - support[1]\n return(c(support[1] + 0.25 * range, \n support[1] + 0.5 * range, \n support[1] + 0.75 * range))\n } else if (is.finite(support[1])) {\n # Left-bounded\n return(c(support[1] + 0.5, support[1] + 1.5, support[1] + 3))\n } else if (is.finite(support[2])) {\n # Right-bounded\n return(c(support[2] - 3, support[2] - 1.5, support[2] - 0.5))\n } else {\n # Unbounded\n return(c(-1, 0, 1))\n }\n}\n\n#' Compute numerical derivatives\ncompute_derivatives <- function(x, f, h = 1e-6, ...) {\n n <- length(x)\n dh <- numeric(n)\n \n for (i in 1:n) {\n f_plus <- f(x[i] + h, ...)\n f_minus <- f(x[i] - h, ...)\n \n if (f_plus <= 0 || f_minus <= 0) {\n stop(\"Density function returns non-positive values near x = \", x[i])\n }\n \n dh[i] <- (log(f_plus) - log(f_minus)) / (2 * h)\n }\n \n return(dh)\n}\n\n#' Check log-concavity\ncheck_log_concavity <- function(x, h, dh) {\n n <- length(x)\n if (n < 2) return(TRUE)\n \n for (i in 1:(n-1)) {\n # Derivative should be decreasing for log-concave functions\n if (dh[i+1] > dh[i] + 1e-6) {\n return(FALSE)\n }\n }\n \n return(TRUE)\n}\n\n#' Check log-concavity when adding a new point\ncheck_log_concavity_update <- function(abscissae, x_new, h_new, dh_new) {\n # Find where to insert\n idx <- findInterval(x_new, abscissae$x)\n \n # Check with left neighbor\n if (idx > 0) {\n if (dh_new > abscissae$dh[idx] + 1e-6) {\n return(FALSE)\n }\n }\n \n # Check with right neighbor\n if (idx < length(abscissae$x)) {\n if (abscissae$dh[idx + 1] > dh_new + 1e-6) {\n return(FALSE)\n }\n }\n \n return(TRUE)\n}\n\n#' Insert new abscissa\ninsert_abscissa <- function(abscissae, x_new, h_new, dh_new) {\n idx <- findInterval(x_new, abscissae$x)\n \n if (idx == 0) {\n abscissae$x <- c(x_new, abscissae$x)\n abscissae$h <- c(h_new, abscissae$h)\n abscissae$dh <- c(dh_new, abscissae$dh)\n } else if (idx == length(abscissae$x)) {\n abscissae$x <- c(abscissae$x, x_new)\n abscissae$h <- c(abscissae$h, h_new)\n abscissae$dh <- c(abscissae$dh, dh_new)\n } else {\n abscissae$x <- c(abscissae$x[1:idx], x_new, abscissae$x[(idx+1):length(abscissae$x)])\n abscissae$h <- c(abscissae$h[1:idx], h_new, abscissae$h[(idx+1):length(abscissae$h)])\n abscissae$dh <- c(abscissae$dh[1:idx], dh_new, abscissae$dh[(idx+1):length(abscissae$dh)])\n }\n \n return(abscissae)\n}\n\n#' Construct piecewise linear envelope\nconstruct_envelope <- function(abscissae, support) {\n k <- length(abscissae$x)\n z <- numeric(k + 1)\n \n # Set boundary points\n z[1] <- support[1]\n z[k + 1] <- support[2]\n \n # Compute interior intersection points\n for (i in 2:k) {\n # Intersection of tangent lines at x[i-1] and x[i]\n dh_diff <- abscissae$dh[i-1] - abscissae$dh[i]\n \n # Handle case where derivatives are nearly equal (parallel tangent lines)\n if (abs(dh_diff) < 1e-8) {\n # Use midpoint when tangent lines are nearly parallel\n z[i] <- (abscissae$x[i-1] + abscissae$x[i]) / 2\n } else {\n z[i] <- (abscissae$h[i] - abscissae$h[i-1] - \n abscissae$x[i] * abscissae$dh[i] + \n abscissae$x[i-1] * abscissae$dh[i-1]) / dh_diff\n }\n \n # Clamp z[i] to support boundaries\n if (is.finite(support[1])) {\n z[i] <- max(z[i], support[1])\n }\n if (is.finite(support[2])) {\n z[i] <- min(z[i], support[2])\n }\n \n # Ensure monotonicity: z[i] must be >= z[i-1]\n if (z[i] <= z[i-1]) {\n # Place z[i] between x[i-1] and x[i] to maintain monotonicity\n z[i] <- (abscissae$x[i-1] + abscissae$x[i]) / 2\n # Ensure it's still greater than z[i-1]\n if (z[i] <= z[i-1]) {\n z[i] <- z[i-1] + (abscissae$x[i] - abscissae$x[i-1]) * 0.1\n }\n }\n }\n \n return(list(\n x = abscissae$x,\n h = abscissae$h,\n dh = abscissae$dh,\n z = z\n ))\n}\n\n#' Evaluate upper envelope at a point\neval_upper_envelope <- function(x, env) {\n k <- length(env$x)\n idx <- findInterval(x, env$z)\n \n if (idx == 0) idx <- 1\n if (idx > k) idx <- k\n \n return(env$h[idx] + env$dh[idx] * (x - env$x[idx]))\n}\n\n#' Evaluate lower envelope (piecewise linear through points)\neval_lower_envelope <- function(x, abscissae) {\n k <- length(abscissae$x)\n \n if (x <= abscissae$x[1] || x >= abscissae$x[k]) {\n return(-Inf)\n }\n \n idx <- findInterval(x, abscissae$x)\n \n if (idx == 0 || idx >= k) {\n return(-Inf)\n }\n \n # Linear interpolation\n slope <- (abscissae$h[idx + 1] - abscissae$h[idx]) / \n (abscissae$x[idx + 1] - abscissae$x[idx])\n \n return(abscissae$h[idx] + slope * (x - abscissae$x[idx]))\n}\n\n#' Sample from upper envelope\nsample_upper_envelope <- function(env, support) {\n k <- length(env$x)\n \n # Compute normalization constants for each piece\n integrals <- numeric(k)\n \n for (i in 1:k) {\n z_left <- env$z[i]\n z_right <- env$z[i + 1]\n \n # Handle infinite boundaries\n if (!is.finite(z_left)) {\n if (env$dh[i] <= 0) {\n stop(\"Cannot have non-positive derivative at -Inf boundary\")\n }\n z_left <- env$x[i] - 10 / env$dh[i]\n }\n \n if (!is.finite(z_right)) {\n if (env$dh[i] >= 0) {\n stop(\"Cannot have non-negative derivative at +Inf boundary\")\n }\n z_right <- env$x[i] - 10 / env$dh[i]\n }\n \n # Ensure z_left < z_right\n if (z_left >= z_right) {\n z_right <- z_left + 0.01\n }\n \n # Compute integral of exp(h + dh*(x - x_i)) from z_left to z_right\n if (abs(env$dh[i]) < 1e-10) {\n # Nearly flat - use rectangular approximation\n integrals[i] <- exp(env$h[i]) * (z_right - z_left)\n } else {\n # Exponential piece: integral of exp(h[i] + dh[i]*(x - x[i])) from z_left to z_right\n h_left <- env$h[i] + env$dh[i] * (z_left - env$x[i])\n h_right <- env$h[i] + env$dh[i] * (z_right - env$x[i])\n \n # Avoid overflow by factoring out max\n h_max <- max(h_left, h_right)\n \n # The integral should always be positive\n integrals[i] <- abs(exp(h_max) * (exp(h_right - h_max) - exp(h_left - h_max)) / env$dh[i])\n }\n \n if (integrals[i] <= 0 || !is.finite(integrals[i])) {\n # Fallback to simple approximation\n integrals[i] <- exp(env$h[i]) * (z_right - z_left)\n }\n }\n \n # Sample piece\n probs <- integrals / sum(integrals)\n piece <- sample(1:k, 1, prob = probs)\n \n # Sample from selected piece\n z_left <- env$z[piece]\n z_right <- env$z[piece + 1]\n \n # Handle infinite boundaries (same as above)\n if (!is.finite(z_left)) {\n z_left <- env$x[piece] - 10 / env$dh[piece]\n }\n if (!is.finite(z_right)) {\n z_right <- env$x[piece] - 10 / env$dh[piece]\n }\n \n # Ensure z_left < z_right\n if (z_left >= z_right) {\n z_right <- z_left + 0.01\n }\n \n if (abs(env$dh[piece]) < 1e-10) {\n # Uniform sampling\n x <- runif(1, z_left, z_right)\n } else {\n # Inverse CDF sampling from exponential\n u <- runif(1)\n h_left <- env$h[piece] + env$dh[piece] * (z_left - env$x[piece])\n h_right <- env$h[piece] + env$dh[piece] * (z_right - env$x[piece])\n \n # Avoid overflow\n h_max <- max(h_left, h_right)\n exp_left_norm <- exp(h_left - h_max)\n exp_right_norm <- exp(h_right - h_max)\n \n h_sample <- h_max + log(exp_left_norm + u * (exp_right_norm - exp_left_norm))\n x <- env$x[piece] + (h_sample - env$h[piece]) / env$dh[piece]\n }\n \n return(x)\n}\n\n################################################################################\n# Testing Function\n################################################################################\n\n#' Test the ARS implementation\n#'\n#' @export\ntest <- function() {\n cat(\"\\n========================================\\n\")\n cat(\"ADAPTIVE REJECTION SAMPLER TEST SUITE\\n\")\n cat(\"========================================\\n\\n\")\n \n all_passed <- TRUE\n \n # Test 1: Standard Normal Distribution\n cat(\"TEST 1: Standard Normal Distribution\\n\")\n cat(\"--------------------------------------\\n\")\n set.seed(123)\n tryCatch({\n samples <- ars(n = 1000, f = dnorm, support = c(-Inf, Inf))\n \n # Save samples\n write.table(samples, \"/app/normal_samples.txt\", \n row.names = FALSE, col.names = FALSE)\n \n mean_val <- mean(samples)\n sd_val <- sd(samples)\n \n cat(sprintf(\" Generated %d samples\\n\", length(samples)))\n cat(sprintf(\" Sample mean: %.4f (expected: 0.0000)\\n\", mean_val))\n cat(sprintf(\" Sample SD: %.4f (expected: 1.0000)\\n\", sd_val))\n \n # Statistical test (mean should be close to 0)\n if (abs(mean_val) < 0.1 && abs(sd_val - 1) < 0.1) {\n cat(\" TEST 1: PASS\\n\\n\")\n } else {\n cat(\" TEST 1: FAIL (statistics out of expected range)\\n\\n\")\n all_passed <- FALSE\n }\n }, error = function(e) {\n cat(\" TEST 1: FAIL -\", e$message, \"\\n\\n\")\n all_passed <<- FALSE\n })\n \n # Test 2: Exponential Distribution\n cat(\"TEST 2: Exponential Distribution (rate=1)\\n\")\n cat(\"------------------------------------------\\n\")\n set.seed(456)\n tryCatch({\n samples <- ars(n = 1000, f = dexp, support = c(0, Inf), rate = 1)\n \n # Save samples\n write.table(samples, \"/app/exponential_samples.txt\", \n row.names = FALSE, col.names = FALSE)\n \n mean_val <- mean(samples)\n sd_val <- sd(samples)\n \n cat(sprintf(\" Generated %d samples\\n\", length(samples)))\n cat(sprintf(\" Sample mean: %.4f (expected: 1.0000)\\n\", mean_val))\n cat(sprintf(\" Sample SD: %.4f (expected: 1.0000)\\n\", sd_val))\n \n if (abs(mean_val - 1) < 0.15 && abs(sd_val - 1) < 0.15) {\n cat(\" TEST 2: PASS\\n\\n\")\n } else {\n cat(\" TEST 2: FAIL (statistics out of expected range)\\n\\n\")\n all_passed <- FALSE\n }\n }, error = function(e) {\n cat(\" TEST 2: FAIL -\", e$message, \"\\n\\n\")\n all_passed <<- FALSE\n })\n \n # Test 3: Gamma Distribution\n cat(\"TEST 3: Gamma Distribution (shape=2, rate=1)\\n\")\n cat(\"--------------------------------------------\\n\")\n set.seed(789)\n tryCatch({\n samples <- ars(n = 1000, f = dgamma, support = c(0, Inf), \n shape = 2, rate = 1)\n \n mean_val <- mean(samples)\n sd_val <- sd(samples)\n \n cat(sprintf(\" Generated %d samples\\n\", length(samples)))\n cat(sprintf(\" Sample mean: %.4f (expected: 2.0000)\\n\", mean_val))\n cat(sprintf(\" Sample SD: %.4f (expected: 1.4142)\\n\", sd_val))\n \n if (abs(mean_val - 2) < 0.2 && abs(sd_val - sqrt(2)) < 0.2) {\n cat(\" TEST 3: PASS\\n\\n\")\n } else {\n cat(\" TEST 3: FAIL (statistics out of expected range)\\n\\n\")\n all_passed <- FALSE\n }\n }, error = function(e) {\n cat(\" TEST 3: FAIL -\", e$message, \"\\n\\n\")\n all_passed <<- FALSE\n })\n \n # Test 4: Beta Distribution (log-concave case)\n cat(\"TEST 4: Beta Distribution (a=2, b=2)\\n\")\n cat(\"-------------------------------------\\n\")\n set.seed(101)\n tryCatch({\n samples <- ars(n = 1000, f = dbeta, support = c(0, 1), \n shape1 = 2, shape2 = 2)\n \n mean_val <- mean(samples)\n sd_val <- sd(samples)\n \n cat(sprintf(\" Generated %d samples\\n\", length(samples)))\n cat(sprintf(\" Sample mean: %.4f (expected: 0.5000)\\n\", mean_val))\n cat(sprintf(\" Sample SD: %.4f (expected: 0.2236)\\n\", sd_val))\n \n if (abs(mean_val - 0.5) < 0.1 && abs(sd_val - sqrt(1/20)) < 0.05) {\n cat(\" TEST 4: PASS\\n\\n\")\n } else {\n cat(\" TEST 4: FAIL (statistics out of expected range)\\n\\n\")\n all_passed <- FALSE\n }\n }, error = function(e) {\n cat(\" TEST 4: FAIL -\", e$message, \"\\n\\n\")\n all_passed <<- FALSE\n })\n \n # Test 5: Input Validation - Negative n\n cat(\"TEST 5: Input Validation - Negative n\\n\")\n cat(\"--------------------------------------\\n\")\n tryCatch({\n samples <- ars(n = -5, f = dnorm, support = c(-Inf, Inf))\n cat(\" TEST 5: FAIL (should have thrown error)\\n\\n\")\n all_passed <- FALSE\n }, error = function(e) {\n cat(\" Correctly rejected negative n\\n\")\n cat(\" TEST 5: PASS\\n\\n\")\n })\n \n # Test 6: Input Validation - Invalid support\n cat(\"TEST 6: Input Validation - Invalid support\\n\")\n cat(\"-------------------------------------------\\n\")\n tryCatch({\n samples <- ars(n = 100, f = dnorm, support = c(5, 2))\n cat(\" TEST 6: FAIL (should have thrown error)\\n\\n\")\n all_passed <- FALSE\n }, error = function(e) {\n cat(\" Correctly rejected invalid support\\n\")\n cat(\" TEST 6: PASS\\n\\n\")\n })\n \n # Test 7: Non-log-concave detection (bimodal distribution)\n cat(\"TEST 7: Non-log-concave Detection\\n\")\n cat(\"----------------------------------\\n\")\n tryCatch({\n # Bimodal distribution - clearly non-log-concave\n bimodal <- function(x) {\n dnorm(x, mean = -2, sd = 0.5) + dnorm(x, mean = 2, sd = 0.5)\n }\n samples <- ars(n = 100, f = bimodal, support = c(-5, 5))\n cat(\" TEST 7: FAIL (should have detected non-log-concavity)\\n\\n\")\n all_passed <- FALSE\n }, error = function(e) {\n if (grepl(\"log-concave\", e$message, ignore.case = TRUE)) {\n cat(\" Correctly detected non-log-concave function\\n\")\n cat(\" TEST 7: PASS\\n\\n\")\n } else {\n cat(\" TEST 7: FAIL (wrong error):\", e$message, \"\\n\\n\")\n all_passed <- FALSE\n }\n })\n \n # Test 8: Non-function input\n cat(\"TEST 8: Input Validation - Non-function\\n\")\n cat(\"----------------------------------------\\n\")\n tryCatch({\n samples <- ars(n = 100, f = \"not a function\", support = c(-Inf, Inf))\n cat(\" TEST 8: FAIL (should have thrown error)\\n\\n\")\n all_passed <- FALSE\n }, error = function(e) {\n cat(\" Correctly rejected non-function input\\n\")\n cat(\" TEST 8: PASS\\n\\n\")\n })\n \n # Test 9: Kolmogorov-Smirnov test for Normal\n cat(\"TEST 9: K-S Test for Normal Distribution\\n\")\n cat(\"-----------------------------------------\\n\")\n set.seed(999)\n tryCatch({\n samples <- ars(n = 500, f = dnorm, support = c(-Inf, Inf))\n ks_result <- ks.test(samples, \"pnorm\")\n \n cat(sprintf(\" K-S statistic: %.4f\\n\", ks_result$statistic))\n cat(sprintf(\" K-S p-value: %.4f\\n\", ks_result$p.value))\n \n if (ks_result$p.value > 0.05) {\n cat(\" TEST 9: PASS\\n\\n\")\n } else {\n cat(\" TEST 9: FAIL (samples don't match normal distribution)\\n\\n\")\n all_passed <- FALSE\n }\n }, error = function(e) {\n cat(\" TEST 9: FAIL -\", e$message, \"\\n\\n\")\n all_passed <<- FALSE\n })\n \n # Summary\n cat(\"========================================\\n\")\n if (all_passed) {\n cat(\"ALL TESTS PASSED\\n\")\n } else {\n cat(\"SOME TESTS FAILED\\n\")\n }\n cat(\"========================================\\n\\n\")\n \n return(invisible(all_passed))\n}\n\n# Run tests if this script is sourced\nif (!interactive()) {\n test()\n}\nEOF\n",
"duration": 0.1
}
New Terminal Output:
root@h200-011-033:/app# Rscript /app/ars.R
========================================
ADAPTIVE REJECTION SAMPLER TEST SUITE
========================================
TEST 1: Standard Normal Distribution
--------------------------------------
Generated 1000 samples
Sample mean: 0.0191 (expected: 0.0000)
Sample SD: 0.9702 (expected: 1.0000)
TEST 1: PASS
TEST 2: Exponential Distribution (rate=1)
------------------------------------------
TEST 2: FAIL - negative probability
TEST 3: Gamma Distribution (shape=2, rate=1)
--------------------------------------------
Generated 1000 samples
Sample mean: 2.0183 (expected: 2.0000)
Sample SD: 1.4379 (expected: 1.4142)
TEST 3: PASS
TEST 4: Beta Distribution (a=2, b=2)
-------------------------------------
Generated 1000 samples
Sample mean: 0.4984 (expected: 0.5000)
Sample SD: 0.2203 (expected: 0.2236)
TEST 4: PASS
TEST 5: Input Validation - Negative n
--------------------------------------
Correctly rejected negative n
TEST 5: PASS
TEST 6: Input Validation - Invalid support
-------------------------------------------
Correctly rejected invalid support
TEST 6: PASS
TEST 7: Non-log-concave Detection
----------------------------------
TEST 7: FAIL (wrong error): negative probability
TEST 8: Input Validation - Non-function
----------------------------------------
Correctly rejected non-function input
TEST 8: PASS
========================================
SOME TESTS FAILED
========================================
root@h200-011-033:/app#
Prompt: 50,853 tokens
Completion: 7,465 tokens