# Clean environment and load libraries rm(list = ls()) # Clear environment to avoid conflicts # Source files library(MixStable) set.seed(123) cat("=== Testing McCulloch Functions ===\n\n") # Test 1: Generate test data with known parameters cat("Test 1: Generating test data\n") n_samples <- 5000 true_alpha <- 1.5 true_beta <- 0.3 true_gamma <- 1.2 true_delta <- 0.5 test_data <- rstable(n_samples, true_alpha, true_beta, true_gamma, true_delta) cat("Generated", length(test_data), "samples\n") cat("True parameters: α =", true_alpha, ", β =", true_beta, ", γ =", true_gamma, ", δ =", true_delta, "\n") cat("Data summary:", summary(test_data), "\n\n") # Test 2: Test compute_quantile_ratios cat("Test 2: Computing quantile ratios\n") ratios <- compute_quantile_ratios(test_data) cat("v_alpha =", round(ratios$v_alpha, 4), "\n") cat("v_beta =", round(ratios$v_beta, 4), "\n") cat("v_gamma =", round(ratios$v_gamma, 4), "\n") cat("v_delta =", round(ratios$v_delta, 4), "\n\n") # Test 3: Test McCulloch quantile initialization cat("Test 3: McCulloch quantile initialization\n") init_params <- mcculloch_quantile_init(test_data) cat("Initialized parameters:\n") cat("α =", round(init_params$alpha, 4), "(true:", true_alpha, ")\n") cat("β =", round(init_params$beta, 4), "(true:", true_beta, ")\n") cat("γ =", round(init_params$gamma, 4), "(true:", true_gamma, ")\n") cat("δ =", round(init_params$delta, 4), "(true:", true_delta, ")\n\n") # Test 4: Generate McCulloch lookup table (small version for testing) cat("Test 4: Generating McCulloch lookup table\n") alpha_grid <- seq(0.8, 2.0, by = 0.2) beta_grid <- seq(-0.5, 0.5, by = 0.25) cat("Grid sizes: α =", length(alpha_grid), "points, β =", length(beta_grid), "points\n") cat("Total combinations:", length(alpha_grid) * length(beta_grid), "\n") # Generate table (with smaller sample size for speed) lookup_table <- generate_mcculloch_table(alpha_grid, beta_grid, size = 1000, seed = 42) cat("Generated lookup table with", length(lookup_table), "entries\n") # Show a few entries cat("Sample entries:\n") for (i in 1:min(3, length(lookup_table))) { key <- names(lookup_table)[i] entry <- lookup_table[[key]] cat(sprintf("Key %s: v_α=%.4f, v_β=%.4f\n", key, entry$v_alpha, entry$v_beta)) } cat("\n") # Test 5: Build interpolators cat("Test 5: Building McCulloch interpolators\n") interpolators <- build_mcculloch_interpolators(lookup_table) cat("Built interpolation functions\n") # Test interpolation with some values test_va <- 2.5 test_vb <- 0.1 interp_alpha_val <- interpolators$interp_alpha(test_va, test_vb) interp_beta_val <- interpolators$interp_beta(test_va, test_vb) cat(sprintf("Test interpolation: v_α=%.2f, v_β=%.2f → α=%.4f, β=%.4f\n", test_va, test_vb, interp_alpha_val, interp_beta_val)) cat("\n") # Test 6: Full McCulloch lookup estimation cat("Test 6: McCulloch lookup estimation\n") lookup_params <- mcculloch_lookup_estimate(test_data,interp_alpha = interpolators$interp_alpha,interp_beta = interpolators$interp_beta) cat("McCulloch lookup results:\n") cat("α =", round(lookup_params$alpha, 4), "(true:", true_alpha, ", error:", round(abs(lookup_params$alpha - true_alpha), 4), ")\n") cat("β =", round(lookup_params$beta, 4), "(true:", true_beta, ", error:", round(abs(lookup_params$beta - true_beta), 4), ")\n") cat("γ =", round(lookup_params$gamma, 4), "(true:", true_gamma, ", error:", round(abs(lookup_params$gamma - true_gamma), 4), ")\n") cat("δ =", round(lookup_params$delta, 4), "(true:", true_delta, ", error:", round(abs(lookup_params$delta - true_delta), 4), ")\n\n") # Test 7: Test with different parameter sets cat("Test 7: Testing with different parameter combinations\n") test_cases <- list( list(alpha = 0.9, beta = -0.8, gamma = 0.5, delta = 1.0, name = "Heavy-tailed, left-skewed"), list(alpha = 1.8, beta = 0.6, gamma = 2.0, delta = -1.0, name = "Light-tailed, right-skewed"), list(alpha = 2.0, beta = 0.0, gamma = 1.0, delta = 0.0, name = "Normal case") ) for (test_case in test_cases) { cat("\n", test_case$name, ":\n") cat("True: α =", test_case$alpha, ", β =", test_case$beta, ", γ =", test_case$gamma, ", δ =", test_case$delta, "\n") # Generate data test_x <- rstable(2000, test_case$alpha, test_case$beta, test_case$gamma, test_case$delta) # Estimate using quantile method est_params <- mcculloch_quantile_init(test_x) cat("Estimated: α =", round(est_params$alpha, 3), ", β =", round(est_params$beta, 3), ", γ =", round(est_params$gamma, 3), ", δ =", round(est_params$delta, 3), "\n") } cat("\n=== McCulloch function tests completed ===\n") cat("Note: Replace the mock rstable() function with a proper stable distribution\n")