Message
Analysis: Perfect! The sample files were created successfully. I can see:
- normal_samples.txt with samples that look reasonable (mix of positive and negative values around 0)
- exponential_samples.txt with positive values
Now I need to fix TEST 5 to use a better non-log-concave function. A good example would be a bimodal distribution or sin(x) which is clearly not log-concave.
Plan: Update the test function to use a clearer non-log-concave example, such as a bimodal distribution (mixture of two normals far apart) or a function with sin() that creates non-monotonic derivatives.
{
"keystrokes": "cat > /app/ars.R << 'EOF'\n# Adaptive Rejection Sampling (ARS) Implementation\n# Based on Gilks & Wild (1992)\n\n#' Adaptive Rejection Sampling\n#' \n#' @param n Number of samples to generate\n#' @param f Function representing the (possibly unnormalized) density\n#' @param support Vector of length 2 specifying the support [lower, upper]\n#' @param start_points Initial abscissae for hull construction (optional)\n#' @return Vector of n samples from the distribution\nars <- function(n, f, support = c(-Inf, Inf), start_points = NULL) {\n \n # Input validation\n validate_inputs(n, f, support)\n \n # Initialize abscissae (points where we evaluate the log-density)\n if (is.null(start_points)) {\n start_points <- initialize_abscissae(support)\n }\n \n # Compute log density and derivatives at initial points\n x <- sort(start_points)\n h <- log_density_safe(f, x)\n hp <- compute_derivative(f, x)\n \n # Check log-concavity at initial points\n check_log_concavity(x, hp)\n \n # Storage for samples\n samples <- numeric(n)\n n_sampled <- 0\n \n # Main sampling loop\n max_iterations <- n * 1000 # Prevent infinite loops\n iterations <- 0\n \n while (n_sampled < n && iterations < max_iterations) {\n iterations <- iterations + 1\n \n # Construct upper and lower hulls\n hull <- construct_hulls(x, h, hp, support)\n \n # Sample from upper hull\n x_star <- sample_upper_hull(hull)\n \n # Evaluate density at sampled point\n h_star <- log_density_safe(f, x_star)\n hp_star <- compute_derivative(f, x_star)\n \n # Compute upper and lower hull values at x_star\n u_star <- evaluate_upper_hull(x_star, hull)\n l_star <- evaluate_lower_hull(x_star, x, h, hp)\n \n # Uniform random variable for rejection test\n w <- runif(1)\n \n # Squeeze test (accept)\n if (l_star > -Inf && w <= exp(l_star - u_star)) {\n n_sampled <- n_sampled + 1\n samples[n_sampled] <- x_star\n } \n # Rejection test\n else if (w <= exp(h_star - u_star)) {\n n_sampled <- n_sampled + 1\n samples[n_sampled] <- x_star\n }\n \n # Update hull with new point\n x <- c(x, x_star)\n h <- c(h, h_star)\n hp <- c(hp, hp_star)\n \n # Re-sort everything\n ord <- order(x)\n x <- x[ord]\n h <- h[ord]\n hp <- hp[ord]\n \n # Check log-concavity with new point\n check_log_concavity(x, hp)\n }\n \n if (iterations >= max_iterations) {\n warning(\"Maximum iterations reached\")\n }\n \n return(samples)\n}\n\n#' Validate inputs to ars function\nvalidate_inputs <- function(n, f, support) {\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 (length(support) != 2 || support[1] >= support[2]) {\n stop(\"support must be a vector of length 2 with support[1] < support[2]\")\n }\n \n # Test that function is vectorized and returns valid values\n if (is.finite(support[1]) && is.finite(support[2])) {\n test_x <- seq(support[1] + 0.1 * (support[2] - support[1]), \n support[2] - 0.1 * (support[2] - support[1]), \n length.out = 3)\n } else if (is.finite(support[1])) {\n test_x <- support[1] + c(0.1, 1, 2)\n } else if (is.finite(support[2])) {\n test_x <- support[2] - c(2, 1, 0.1)\n } else {\n test_x <- c(-1, 0, 1)\n }\n \n tryCatch({\n result <- f(test_x)\n if (length(result) != length(test_x)) {\n stop(\"f must be vectorized (return same length as input)\")\n }\n if (any(result < 0, na.rm = TRUE)) {\n stop(\"f must return non-negative values (it's a density)\")\n }\n if (any(is.na(result))) {\n stop(\"f returned NA values\")\n }\n }, error = function(e) {\n stop(paste(\"Error evaluating function f:\", e$message))\n })\n}\n\n#' Initialize abscissae for hull construction\ninitialize_abscissae <- function(support) {\n if (is.finite(support[1]) && is.finite(support[2])) {\n # Bounded support\n width <- support[2] - support[1]\n return(seq(support[1] + 0.1 * width, support[2] - 0.1 * width, length.out = 5))\n } else if (is.finite(support[1])) {\n # Left bounded\n return(support[1] + c(0.1, 1, 2, 5))\n } else if (is.finite(support[2])) {\n # Right bounded\n return(support[2] - c(5, 2, 1, 0.1))\n } else {\n # Unbounded\n return(c(-2, -1, 0, 1, 2))\n }\n}\n\n#' Safely compute log density\nlog_density_safe <- function(f, x) {\n vals <- f(x)\n if (any(vals < 0, na.rm = TRUE)) {\n stop(\"Density function returned negative values\")\n }\n if (any(is.na(vals))) {\n stop(\"Density function returned NA values\")\n }\n \n # Handle zero values\n vals[vals == 0] <- 1e-300\n \n log_vals <- log(vals)\n if (any(is.nan(log_vals))) {\n stop(\"Density function returned values that produce NaN in log\")\n }\n \n return(log_vals)\n}\n\n#' Compute numerical derivative of log density\ncompute_derivative <- function(f, x, h = 1e-5) {\n n <- length(x)\n deriv <- numeric(n)\n \n for (i in 1:n) {\n # Use adaptive step size\n step <- h * max(abs(x[i]), 1)\n \n # Central difference\n f_plus <- f(x[i] + step)\n f_minus <- f(x[i] - step)\n \n if (f_plus <= 0 || f_minus <= 0 || is.na(f_plus) || is.na(f_minus)) {\n # Try forward or backward difference\n f_0 <- f(x[i])\n if (f_0 > 0 && f_plus > 0) {\n deriv[i] <- (log(f_plus) - log(f_0)) / step\n } else if (f_0 > 0 && f_minus > 0) {\n deriv[i] <- (log(f_0) - log(f_minus)) / step\n } else {\n stop(\"Cannot compute derivative - density is zero or negative near evaluation point\")\n }\n } else {\n deriv[i] <- (log(f_plus) - log(f_minus)) / (2 * step)\n }\n }\n \n return(deriv)\n}\n\n#' Check log-concavity\ncheck_log_concavity <- function(x, hp) {\n if (length(x) < 2) return(TRUE)\n \n # Derivatives should be non-increasing (with tolerance for numerical errors)\n tolerance <- 1e-4\n for (i in 2:length(hp)) {\n if (hp[i] > hp[i-1] + tolerance) {\n stop(\"Function is not log-concave: derivatives are not decreasing\")\n }\n }\n \n return(TRUE)\n}\n\n#' Construct upper and lower hulls\nconstruct_hulls <- function(x, h, hp, support) {\n n <- length(x)\n \n # Compute intersection points z of tangent lines\n z <- numeric(n + 1)\n z[1] <- support[1]\n z[n + 1] <- support[2]\n \n for (i in 2:n) {\n # Intersection of tangents at x[i-1] and x[i]\n if (abs(hp[i-1] - hp[i]) < 1e-10) {\n # Nearly parallel tangents\n z[i] <- (x[i-1] + x[i]) / 2\n } else {\n z[i] <- (h[i] - h[i-1] - x[i]*hp[i] + x[i-1]*hp[i-1]) / (hp[i-1] - hp[i])\n }\n }\n \n list(x = x, h = h, hp = hp, z = z)\n}\n\n#' Sample from upper hull\nsample_upper_hull <- function(hull) {\n n <- length(hull$x)\n \n # Compute unnormalized probabilities for each segment using log-sum-exp\n log_probs <- numeric(n)\n \n for (i in 1:n) {\n z_left <- hull$z[i]\n z_right <- hull$z[i + 1]\n \n if (is.infinite(z_left) && is.infinite(z_right)) {\n log_probs[i] <- -Inf\n next\n }\n \n if (abs(hull$hp[i]) < 1e-10) {\n # Nearly flat segment\n if (is.finite(z_left) && is.finite(z_right)) {\n log_probs[i] <- hull$h[i] + log(z_right - z_left)\n } else {\n log_probs[i] <- -Inf\n }\n } else {\n # Exponential segment - use log-sum-exp for stability\n # Integral is exp(h[i] - x[i]*hp[i]) * (exp(hp[i]*z_right) - exp(hp[i]*z_left)) / hp[i]\n \n if (is.infinite(z_right) && hull$hp[i] >= 0) {\n log_probs[i] <- -Inf # Doesn't converge\n } else if (is.infinite(z_left) && hull$hp[i] <= 0) {\n log_probs[i] <- -Inf # Doesn't converge\n } else {\n # Compute log of integral\n term1 <- hull$h[i] - hull$x[i] * hull$hp[i]\n \n if (is.infinite(z_right)) {\n # z_right is Inf, hp < 0\n log_diff <- hull$hp[i] * z_left\n } else if (is.infinite(z_left)) {\n # z_left is -Inf, hp > 0\n log_diff <- hull$hp[i] * z_right\n } else {\n # Both finite - use log-diff-exp\n v1 <- hull$hp[i] * z_right\n v2 <- hull$hp[i] * z_left\n max_v <- max(v1, v2)\n log_diff <- max_v + log(abs(exp(v1 - max_v) - exp(v2 - max_v)))\n }\n \n log_probs[i] <- term1 + log_diff - log(abs(hull$hp[i]))\n }\n }\n }\n \n # Normalize probabilities using log-sum-exp\n max_log_prob <- max(log_probs[is.finite(log_probs)])\n probs <- exp(log_probs - max_log_prob)\n probs[!is.finite(probs)] <- 0\n probs <- probs / sum(probs)\n \n # Sample segment\n segment <- sample(1:n, 1, prob = probs)\n \n # Sample within segment using inverse CDF\n u <- runif(1)\n z_left <- hull$z[segment]\n z_right <- hull$z[segment + 1]\n hp_seg <- hull$hp[segment]\n \n if (abs(hp_seg) < 1e-10) {\n # Uniform in segment\n x_star <- z_left + u * (z_right - z_left)\n } else {\n # Inverse CDF of exponential\n if (is.infinite(z_left)) {\n # Left is -Inf, hp must be positive\n x_star <- z_right + log(u) / hp_seg\n } else if (is.infinite(z_right)) {\n # Right is Inf, hp must be negative\n x_star <- z_left + log(1 - u) / hp_seg\n } else {\n # Both finite\n v_left <- hp_seg * z_left\n v_right <- hp_seg * z_right\n max_v <- max(v_left, v_right)\n exp_left_norm <- exp(v_left - max_v)\n exp_right_norm <- exp(v_right - max_v)\n \n x_star <- (max_v + log(exp_left_norm + u * (exp_right_norm - exp_left_norm))) / hp_seg\n }\n }\n \n return(x_star)\n}\n\n#' Evaluate upper hull at point\nevaluate_upper_hull <- function(x_star, hull) {\n # Find which segment x_star is in\n segment <- findInterval(x_star, hull$z)\n if (segment == 0) segment <- 1\n if (segment > length(hull$x)) segment <- length(hull$x)\n \n # Evaluate tangent line at x[segment]\n u <- hull$h[segment] + hull$hp[segment] * (x_star - hull$x[segment])\n \n return(u)\n}\n\n#' Evaluate lower hull at point\nevaluate_lower_hull <- function(x_star, x, h, hp) {\n # Find surrounding points\n idx <- findInterval(x_star, x)\n \n if (idx == 0 || idx == length(x)) {\n return(-Inf) # Outside convex hull\n }\n \n # Linear interpolation between h[idx] and h[idx+1]\n l <- ((x[idx+1] - x_star) * h[idx] + (x_star - x[idx]) * h[idx+1]) / \n (x[idx+1] - x[idx])\n \n return(l)\n}\n\n#' Test function for ARS implementation\ntest <- function() {\n cat(\"=== Testing Adaptive Rejection Sampling ===\\n\\n\")\n \n all_passed <- TRUE\n \n # Test 1: Standard Normal Distribution\n cat(\"TEST 1: Standard Normal Distribution\\n\")\n set.seed(123)\n n <- 10000\n samples <- ars(n, dnorm, c(-Inf, Inf))\n \n mean_samples <- mean(samples)\n sd_samples <- sd(samples)\n \n cat(sprintf(\" Generated %d samples\\n\", n))\n cat(sprintf(\" Sample mean: %.4f (expected: 0.0000)\\n\", mean_samples))\n cat(sprintf(\" Sample SD: %.4f (expected: 1.0000)\\n\", sd_samples))\n \n # Statistical test: mean should be close to 0\n if (abs(mean_samples) < 0.05 && abs(sd_samples - 1) < 0.05) {\n cat(\" TEST 1: PASS\\n\\n\")\n } else {\n cat(\" TEST 1: FAIL\\n\\n\")\n all_passed <- FALSE\n }\n \n # Save samples\n write.table(samples, \"/app/normal_samples.txt\", row.names = FALSE, col.names = FALSE)\n \n # Test 2: Exponential Distribution (rate = 1)\n cat(\"TEST 2: Exponential Distribution (rate=1)\\n\")\n set.seed(456)\n samples <- ars(n, dexp, c(0, Inf))\n \n mean_samples <- mean(samples)\n sd_samples <- sd(samples)\n \n cat(sprintf(\" Generated %d samples\\n\", n))\n cat(sprintf(\" Sample mean: %.4f (expected: 1.0000)\\n\", mean_samples))\n cat(sprintf(\" Sample SD: %.4f (expected: 1.0000)\\n\", sd_samples))\n \n if (abs(mean_samples - 1) < 0.05 && abs(sd_samples - 1) < 0.05) {\n cat(\" TEST 2: PASS\\n\\n\")\n } else {\n cat(\" TEST 2: FAIL\\n\\n\")\n all_passed <- FALSE\n }\n \n # Save samples\n write.table(samples, \"/app/exponential_samples.txt\", row.names = FALSE, col.names = FALSE)\n \n # Test 3: Beta Distribution (alpha=2, beta=5)\n cat(\"TEST 3: Beta Distribution (alpha=2, beta=5)\\n\")\n set.seed(789)\n samples <- ars(n, function(x) dbeta(x, 2, 5), c(0, 1))\n \n mean_samples <- mean(samples)\n expected_mean <- 2 / (2 + 5)\n expected_sd <- sqrt(2 * 5 / ((2 + 5)^2 * (2 + 5 + 1)))\n sd_samples <- sd(samples)\n \n cat(sprintf(\" Generated %d samples\\n\", n))\n cat(sprintf(\" Sample mean: %.4f (expected: %.4f)\\n\", mean_samples, expected_mean))\n cat(sprintf(\" Sample SD: %.4f (expected: %.4f)\\n\", sd_samples, expected_sd))\n \n if (abs(mean_samples - expected_mean) < 0.01 && abs(sd_samples - expected_sd) < 0.01) {\n cat(\" TEST 3: PASS\\n\\n\")\n } else {\n cat(\" TEST 3: FAIL\\n\\n\")\n all_passed <- FALSE\n }\n \n # Test 4: Input Validation - Negative n\n cat(\"TEST 4: Input Validation - Negative n\\n\")\n tryCatch({\n ars(-5, dnorm, c(-Inf, Inf))\n cat(\" TEST 4: 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 4: PASS\\n\\n\")\n })\n \n # Test 5: Input Validation - Non-log-concave function\n cat(\"TEST 5: Log-concavity Check - Non-log-concave function\\n\")\n tryCatch({\n # Bimodal distribution (mixture) - clearly not log-concave\n bimodal <- function(x) 0.5 * dnorm(x, -3, 0.5) + 0.5 * dnorm(x, 3, 0.5)\n ars(100, bimodal, c(-Inf, Inf))\n cat(\" TEST 5: FAIL (should have detected non-log-concavity)\\n\\n\")\n all_passed <- FALSE\n }, error = function(e) {\n if (grepl(\"log-concave\", e$message)) {\n cat(\" Correctly detected non-log-concavity\\n\")\n cat(\" TEST 5: PASS\\n\\n\")\n } else {\n cat(sprintf(\" Error: %s\\n\", e$message))\n cat(\" TEST 5: FAIL (wrong error message)\\n\\n\")\n all_passed <- FALSE\n }\n })\n \n # Test 6: Gamma Distribution (shape=3, rate=2)\n cat(\"TEST 6: Gamma Distribution (shape=3, rate=2)\\n\")\n set.seed(321)\n samples <- ars(n, function(x) dgamma(x, shape=3, rate=2), c(0, Inf))\n \n mean_samples <- mean(samples)\n expected_mean <- 3 / 2\n expected_sd <- sqrt(3) / 2\n sd_samples <- sd(samples)\n \n cat(sprintf(\" Generated %d samples\\n\", n))\n cat(sprintf(\" Sample mean: %.4f (expected: %.4f)\\n\", mean_samples, expected_mean))\n cat(sprintf(\" Sample SD: %.4f (expected: %.4f)\\n\", sd_samples, expected_sd))\n \n if (abs(mean_samples - expected_mean) < 0.05 && abs(sd_samples - expected_sd) < 0.05) {\n cat(\" TEST 6: PASS\\n\\n\")\n } else {\n cat(\" TEST 6: FAIL\\n\\n\")\n all_passed <- FALSE\n }\n \n # Summary\n cat(\"=== Test Summary ===\\n\")\n if (all_passed) {\n cat(\"ALL TESTS PASSED\\n\")\n } else {\n cat(\"SOME TESTS FAILED\\n\")\n }\n \n invisible(all_passed)\n}\nEOF\n",
"duration": 0.1
}
Current Terminal Screen:
> } else {
> cat(sprintf(" Error: %s\n", e$message))
> cat(" TEST 5: FAIL (wrong error message)\n\n")
> all_passed <- FALSE
> }
> })
>
> # Test 6: Gamma Distribution (shape=3, rate=2)
> cat("TEST 6: Gamma Distribution (shape=3, rate=2)\n")
> set.seed(321)
> samples <- ars(n, function(x) dgamma(x, shape=3, rate=2), c(0, Inf))
>
> mean_samples <- mean(samples)
> expected_mean <- 3 / 2
> expected_sd <- sqrt(3) / 2
> sd_samples <- sd(samples)
>
> cat(sprintf(" Generated %d samples\n", n))
> cat(sprintf(" Sample mean: %.4f (expected: %.4f)\n", mean_samples, expected_mean))
> cat(sprintf(" Sample SD: %.4f (expected: %.4f)\n", sd_samples, expected_sd))
>
> if (abs(mean_samples - expected_mean) < 0.05 && abs(sd_samples - expected_sd) < 0.05) {
> cat(" TEST 6: PASS\n\n")
> } else {
> cat(" TEST 6: FAIL\n\n")
> all_passed <- FALSE
> }
>
> # Summary
> cat("=== Test Summary ===\n")
> if (all_passed) {
> cat("ALL TESTS PASSED\n")
> } else {
> cat("SOME TESTS FAILED\n")
> }
>
> invisible(all_passed)
> }
> EOF
root@h200-006-219:/app#
Prompt: 48,153 tokens
Completion: 5,871 tokens