# Clean environment and load libraries
rm(list = ls())  # Clear environment to avoid conflicts


# Source files
library(MixStable)

# Safe mean extractor with NA-check
safe_mean <- function(x) {
  vals <- unlist(x)
  if (!is.numeric(vals) || all(is.na(vals)) || length(vals) == 0) return(NA_real_)
  return(mean(vals, na.rm = TRUE))
}
# === Run Metropolis-Hastings and plot ===
set.seed(123)
xx <- seq(-8, 11, length.out = 200)
xx_true <- seq(-8, 11, length.out = 200)
yy_true <- r_stable_pdf(xx_true, 1.5, -0.5, 4, 5)

# Generate synthetic data
lok <- rnorm(1000, mean = 0, sd = 1)

# ⚠️ Increase iterations to ensure sufficient samples
mh_result <- metropolis_hastings(CDF, iterations = 1000, lok = lok, aa = c(1, 1), proposal_std = 0.1)


# Plot results using final MCMC samples
plot_results(mh_result$M2_w1, mh_result$M2_alpha1, mh_result$M2_beta1,
             mh_result$M2_delta1, mh_result$M2_omega1,
             mh_result$M2_w2, mh_result$M2_alpha2, mh_result$M2_beta2,
             mh_result$M2_delta2, mh_result$M2_omega2,
             xx, xx_true, yy_true)