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 > > # Generate synthetic data > data <- generate_synthetic_data() > > # Estimate bandwidth (if needed) > h_n <- tryCatch({ + stats::bw.SJ(data) + }, error = function(e) { + message("Bandwidth estimation failed: ", e$message) + 1.0 + }) > > # Run EM methods with error handling > safe_run_method <- function(method_func, data, method_name) { + cat("Running", method_name, "...\n") + tryCatch({ + result <- method_func(data) + cat("✓", method_name, "completed successfully\n") + return(result) + }, error = function(e) { + cat("✗", method_name, "failed:", e$message, "\n") + return(NULL) + }) + } > > methods <- list() > methods[["EM_Estimated"]] <- safe_run_method(em_fit_alpha_stable_mixture, data, "EM_Estimated") Running EM_Estimated ... Iteration 1: Log-Likelihood = -2716.034445 Iteration 2: Log-Likelihood = -2713.563292 Iteration 3: Log-Likelihood = -2713.350010 Converged. ✓ EM_Estimated completed successfully There were 50 or more warnings (use warnings() to see the first 50) > methods[["EM_ECF_recursive_ecf"]] <- safe_run_method(em_estimate_stable_recursive_ecf, data, "EM_ECF_recursive_ecf") Running EM_ECF_recursive_ecf ... Converged after 8 iterations. ✓ EM_ECF_recursive_ecf completed successfully > methods[["EM_ECF_Kernel"]] <- safe_run_method(em_estimate_stable_kernel_ecf, data, "EM_ECF_Kernel") Running EM_ECF_Kernel ... Converged after 10 iterations. ✓ EM_ECF_Kernel completed successfully > methods[["EM_ECF_weighted_ols"]] <- safe_run_method(em_estimate_stable_weighted_ols, data, "EM_ECF_weighted_ols") Running EM_ECF_weighted_ols ... Converged after 8 iterations. ✓ EM_ECF_weighted_ols completed successfully > methods[["EM_ECF_from_cdf"]] <- safe_run_method(em_estimate_stable_from_cdf, data, "EM_ECF_from_cdf") Running EM_ECF_from_cdf ... Converged after 7 iterations. ✓ EM_ECF_from_cdf completed successfully > > # Remove failed methods > methods <- methods[!sapply(methods, is.null)] > > if (length(methods) == 0) { + stop("All methods failed. Check your data and functions.") + } > > to_list <- function(p) { + if (is.list(p)) return(p) + names(p) <- c("alpha", "beta", "gamma", "delta") + as.list(p) + } > > # Plotting > x_vals <- seq(min(data), max(data), length.out = 1000) > > for (name in names(methods)) { + result <- methods[[name]] + + # Skip if result is NULL or doesn't have required structure + if (is.null(result) || is.null(result$params1) || is.null(result$params2)) { + cat("Skipping", name, "- invalid result structure\n") + next + } + + tryCatch({ + param1 <- to_list(result$params1) + param2 <- to_list(result$params2) + weight <- result$weight + + # Validate weight + if (is.null(weight) || !is.numeric(weight) || weight < 0 || weight > 1) { + cat("Skipping", name, "- invalid weight\n") + next + } + + p1_vals <- unpack_params(param1) + p2_vals <- unpack_params(param2) + + # Use explicit parameter names to avoid confusion + pdf_comp1 <- r_stable_pdf(x_vals, + alpha = param1$alpha, + beta = param1$beta, + scale = param1$gamma, + location = param1$delta) + + pdf_comp2 <- r_stable_pdf(x_vals, + alpha = param2$alpha, + beta = param2$beta, + scale = param2$gamma, + location = param2$delta) + + mixture_pdf <- weight * pdf_comp1 + (1 - weight) * pdf_comp2 + + cat(sprintf("\n🔍 Plotting %s\n", name)) + cat(sprintf("Component 1: α=%.3f, β=%.3f, γ=%.3f, δ=%.3f\n", + param1$alpha, param1$beta, param1$gamma, param1$delta)) + cat(sprintf("Component 2: α=%.3f, β=%.3f, γ=%.3f, δ=%.3f\n", + param2$alpha, param2$beta, param2$gamma, param2$delta)) + cat(sprintf("Weight: %.3f\n", weight)) + + plot_title <- sprintf("%s vs Data", gsub("_", " ", name)) + png_filename <- sprintf("mixture_alpha_stable_%s.png", tolower(name)) + + png(png_filename, width = 900, height = 500) + hist(data, breaks = 40, freq = FALSE, col = rgb(0.7, 0.7, 0.7, 0.4), + main = plot_title, xlab = "x", ylab = "Density") + lines(x_vals, mixture_pdf, col = "blue", lwd = 2) + legend("topright", + legend = c("Data", "Estimated Mixture"), + col = c("gray", "blue"), + lty = c(NA, 1), + pch = c(15, NA), + bty = "n") + dev.off() + + cat("✓ Plot saved as", png_filename, "\n") + + }, error = function(e) { + cat("Error plotting", name, ":", e$message, "\n") + }) + } Skipping EM_Estimated - invalid weight 🔍 Plotting EM_ECF_recursive_ecf Component 1: α=1.434, β=-0.139, γ=1.620, δ=4.987 Component 2: α=1.679, β=1.000, γ=0.760, δ=-4.226 Weight: 0.558 ✓ Plot saved as mixture_alpha_stable_em_ecf_recursive_ecf.png 🔍 Plotting EM_ECF_Kernel Component 1: α=1.561, β=-0.672, γ=1.818, δ=5.005 Component 2: α=2.000, β=1.000, γ=0.711, δ=-4.189 Weight: 0.584 ✓ Plot saved as mixture_alpha_stable_em_ecf_kernel.png 🔍 Plotting EM_ECF_weighted_ols Component 1: α=1.679, β=1.000, γ=0.760, δ=-4.227 Component 2: α=1.434, β=-0.113, γ=1.620, δ=4.973 Weight: 0.442 ✓ Plot saved as mixture_alpha_stable_em_ecf_weighted_ols.png 🔍 Plotting EM_ECF_from_cdf Component 1: α=1.881, β=-1.000, γ=2.212, δ=4.885 Component 2: α=1.504, β=-1.000, γ=0.502, δ=-4.271 Weight: 0.622 ✓ Plot saved as mixture_alpha_stable_em_ecf_from_cdf.png > > proc.time() user system elapsed 909.60 12.87 913.68