R Under development (unstable) (2025-09-21 r88861 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2025 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> # Clean environment and load libraries
> rm(list = ls())  # Clear environment to avoid conflicts
> 
> # Source files
> library(MixStable)

Attaching package: 'MixStable'

The following object is masked from 'package:graphics':

    clip

> 
> set.seed(123)
> 
> cat("=== Testing McCulloch Functions ===\n\n")
=== Testing McCulloch Functions ===

> 
> # Test 1: Generate test data with known parameters
> cat("Test 1: Generating test data\n")
Test 1: Generating test data
> 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")
Generated 5000 samples
> cat("True parameters: α =", true_alpha, ", β =", true_beta, ", γ =", true_gamma, ", δ =", true_delta, "\n")
True parameters: α = 1.5 , β = 0.3 , γ = 1.2 , δ = 0.5 
> cat("Data summary:", summary(test_data), "\n\n")
Data summary: -233.3332 -0.8735327 0.2190402 0.5226617 1.446061 498.4689 

> 
> # Test 2: Test compute_quantile_ratios
> cat("Test 2: Computing quantile ratios\n")
Test 2: Computing quantile ratios
> ratios <- compute_quantile_ratios(test_data)
> cat("v_alpha =", round(ratios$v_alpha, 4), "\n")
v_alpha = 3.1529 
> cat("v_beta =", round(ratios$v_beta, 4), "\n")
v_beta = 0.122 
> cat("v_gamma =", round(ratios$v_gamma, 4), "\n")
v_gamma = 2.3196 
> cat("v_delta =", round(ratios$v_delta, 4), "\n\n")
v_delta = -0.219 

> 
> # Test 3: Test McCulloch quantile initialization
> cat("Test 3: McCulloch quantile initialization\n")
Test 3: McCulloch quantile initialization
> init_params <- mcculloch_quantile_init(test_data)
> cat("Initialized parameters:\n")
Initialized parameters:
> cat("α =", round(init_params$alpha, 4), "(true:", true_alpha, ")\n")
α = 1.7707 (true: 1.5 )
> cat("β =", round(init_params$beta, 4), "(true:", true_beta, ")\n")
β = 0.122 (true: 0.3 )
> cat("γ =", round(init_params$gamma, 4), "(true:", true_gamma, ")\n")
γ = 2.3196 (true: 1.2 )
> cat("δ =", round(init_params$delta, 4), "(true:", true_delta, ")\n\n")
δ = 0.219 (true: 0.5 )

> 
> # Test 4: Generate McCulloch lookup table (small version for testing)
> cat("Test 4: Generating McCulloch lookup table\n")
Test 4: Generating McCulloch lookup table
> 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")
Grid sizes: α = 7 points, β = 5 points
> cat("Total combinations:", length(alpha_grid) * length(beta_grid), "\n")
Total combinations: 35 
> 
> # 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")
Generated lookup table with 35 entries
> 
> # Show a few entries
> cat("Sample entries:\n")
Sample entries:
> 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))
+ }
Key 0.8_-0.5: v_α=7.8943, v_β=-0.7015
Key 0.8_-0.25: v_α=9.8742, v_β=-0.3753
Key 0.8_0: v_α=10.6762, v_β=-0.0092
> cat("\n")

> 
> # Test 5: Build interpolators
> cat("Test 5: Building McCulloch interpolators\n")
Test 5: Building McCulloch interpolators
> interpolators <- build_mcculloch_interpolators(lookup_table)
> cat("Built interpolation functions\n")
Built interpolation functions
> 
> # 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))
Test interpolation: v_α=2.50, v_β=0.10 → α=2.0000, β=0.5000
> cat("\n")

> 
> # Test 6: Full McCulloch lookup estimation
> cat("Test 6: McCulloch lookup estimation\n")
Test 6: McCulloch lookup estimation
> lookup_params <- mcculloch_lookup_estimate(test_data,interp_alpha = interpolators$interp_alpha,interp_beta = interpolators$interp_beta)
> 
> cat("McCulloch lookup results:\n")
McCulloch lookup results:
> cat("α =", round(lookup_params$alpha, 4), "(true:", true_alpha, ", error:",
+     round(abs(lookup_params$alpha - true_alpha), 4), ")\n")
α = 1.4 (true: 1.5 , error: 0.1 )
> cat("β =", round(lookup_params$beta, 4), "(true:", true_beta, ", error:",
+     round(abs(lookup_params$beta - true_beta), 4), ")\n")
β = 0.5 (true: 0.3 , error: 0.2 )
> cat("γ =", round(lookup_params$gamma, 4), "(true:", true_gamma, ", error:",
+     round(abs(lookup_params$gamma - true_gamma), 4), ")\n")
γ = 2.3196 (true: 1.2 , error: 1.1196 )
> cat("δ =", round(lookup_params$delta, 4), "(true:", true_delta, ", error:",
+     round(abs(lookup_params$delta - true_delta), 4), ")\n\n")
δ = 0.219 (true: 0.5 , error: 0.281 )

> 
> # Test 7: Test with different parameter sets
> cat("Test 7: Testing with different parameter combinations\n")
Test 7: Testing with different parameter combinations
> 
> 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")
+ }

 Heavy-tailed, left-skewed :
True: α = 0.9 , β = -0.8 , γ = 0.5 , δ = 1 
Estimated: α = 0.6 , β = -0.786 , γ = 1.517 , δ = -1.757 

 Light-tailed, right-skewed :
True: α = 1.8 , β = 0.6 , γ = 2 , δ = -1 
Estimated: α = 2 , β = 0.07 , γ = 3.788 , δ = -1.237 

 Normal case :
True: α = 2 , β = 0 , γ = 1 , δ = 0 
Estimated: α = 2 , β = 0.003 , γ = 1.897 , δ = -0.046 
> 
> cat("\n=== McCulloch function tests completed ===\n")

=== McCulloch function tests completed ===
> cat("Note: Replace the mock rstable() function with a proper stable distribution\n")
Note: Replace the mock rstable() function with a proper stable distribution
> 
> proc.time()
   user  system elapsed 
   2.20    0.34    2.53