Code for the distribution "A New Topp--Leone Juchez Distribution with Applications to Accelerated Life-Test and Aircraft Window-Glass Strength Data" are provided:

Code for the distribution "A New Topp--Leone Juchez Distribution with Applications to Accelerated Life-Test and Aircraft Window-Glass Strength Data" are provided:rm(list = ls()) ############################################################# TOPP-LEONE JUCHEZ (TLJ) DISTRIBUTION - SECOND PARAMETER SCENARIO# Complete reproducible analysis for Scenario 2# True values: rho = 0.8 and alpha = 1.5# Monte Carlo replications: 500# Sample sizes: 25, 50, 100, 200, and 500# Parameters: rho > 0, alpha > 0; support: x > 0############################################################ if (!requireNamespace("coda", quietly = TRUE)) install.packages("coda")library(coda) round_df <- function(df, digits = 6) { df[] <- lapply(df, function(z) if (is.numeric(z)) round(z, digits) else z) df} ############################################################# USER SETTINGS############################################################ BAYES_NITER <- 50000BAYES_BURN <- 10000BAYES_THIN <- 10PROP_SD_MAIN <- c(rho = 0.055, alpha = 0.070) # Reproducible multi-chain diagnostic settingsMCMC_NCHAINS <- 4MCMC_DIAG_NITER <- 60000MCMC_DIAG_BURN <- 10000MCMC_DIAG_THIN <- 5MCMC_STARTS <- list( c(rho = 0.60, alpha = 0.65), c(rho = 0.90, alpha = 1.00), c(rho = 1.50, alpha = 1.70), c(rho = 2.00, alpha = 2.20)) MLE_LOWER <- c(rho = 1e-4, alpha = 1e-4)MLE_UPPER <- c(rho = 100, alpha = 100)MLE_NSTART <- 12MLE_HESSIAN_CONDITION_MAX <- 1e10 # Independent Gamma(shape=2, rate=1) priorsPRIOR_SHAPE <- c(rho = 2, alpha = 2)PRIOR_RATE <- c(rho = 1, alpha = 1) # Reviewer-requested final computations# TRUE uses the manuscript-quality settings reported in the paper.FINAL_REVIEWER_RUN <- TRUEMC_REPLICATIONS <- if (FINAL_REVIEWER_RUN) 500 else 100BOOTSTRAP_GOF_B <- 5000FISHER_MC_SIZE <- if (FINAL_REVIEWER_RUN) 50000 else 10000STRESS_STRENGTH_MC_R <- if (FINAL_REVIEWER_RUN) 500 else 30STRESS_STRENGTH_BOOT_B <- if (FINAL_REVIEWER_RUN) 1000 else 200RUN_FISHER_STABILITY <- TRUERUN_BAYES_ROBUSTNESS <- TRUERUN_STRESS_STRENGTH <- TRUE OUT_DIR <- "TLJ_Scenario2_rho0p8_alpha1p5_Outputs"if (!dir.exists(OUT_DIR)) dir.create(OUT_DIR, recursive = TRUE)OLD_WD <- getwd()setwd(OUT_DIR) ############################################################# TLJ DISTRIBUTION: CDF, PDF, SURVIVAL, HAZARD AND QUANTILE############################################################ # Baseline Juchez polynomial and log-survival:# S_J(x;rho) = exp(-rho*x) P_rho(x)/(rho^3+rho^2+6).tlj_logS_J <- function(x, rho) { D <- rho^3 + rho^2 + 6 P <- rho^3*x^3 + 3*rho^2*x^2 + (rho^3 + 6*rho)*x + D -rho*x + log(P) - log(D)} tlj_logg_J <- function(x, rho) { D <- rho^3 + rho^2 + 6 4*log(rho) + log(x^3 + x + 1) - rho*x - log(D)} pTLJ <- function(q, rho, alpha, lower.tail = TRUE, log.p = FALSE) { q <- as.numeric(q) if (length(rho) != 1L || length(alpha) != 1L || !is.finite(rho) || !is.finite(alpha) || rho <= 0 || alpha <= 0) return(rep(NA_real_, length(q))) ans <- numeric(length(q)) ans[q <= 0] <- 0 ok <- is.finite(q) & q > 0 if (any(ok)) { logS <- tlj_logS_J(q[ok], rho) S2 <- pmin(exp(2*logS), 1 - 1e-15) logB <- log1p(-S2) ans[ok] <- exp(alpha*logB) } ans[q == Inf] <- 1 ans <- pmin(pmax(ans, 0), 1) if (!lower.tail) ans <- 1 - ans if (log.p) ans <- log(ans) ans} dTLJ <- function(x, rho, alpha, log = FALSE) { x <- as.numeric(x) out <- if (log) rep(-Inf, length(x)) else numeric(length(x)) if (length(rho) != 1L || length(alpha) != 1L || !is.finite(rho) || !is.finite(alpha) || rho <= 0 || alpha <= 0) return(out) ok <- is.finite(x) & x > 0 if (any(ok)) { xx <- x[ok] logS <- tlj_logS_J(xx, rho) logg <- tlj_logg_J(xx, rho) S2 <- pmin(exp(2*logS), 1 - 1e-15) logB <- log1p(-S2) logf <- log(2) + log(alpha) + logg + logS + (alpha - 1)*logB out[ok] <- if (log) logf else exp(logf) } if (any(x == 0)) { if (alpha < 1) out[x == 0] <- Inf if (alpha == 1) { f0 <- 2*rho^4/(rho^3 + rho^2 + 6) out[x == 0] <- if (log) log(f0) else f0 } } out} sTLJ <- function(x, rho, alpha, log.p = FALSE) { pTLJ(x, rho, alpha, lower.tail = FALSE, log.p = log.p)} hTLJ <- function(x, rho, alpha) { logh <- dTLJ(x, rho, alpha, log = TRUE) - sTLJ(x, rho, alpha, log.p = TRUE) ans <- exp(logh) ans[!is.finite(ans)] <- NA_real_ ans} # F_TLJ(x)=p implies S_J(x)=sqrt(1-p^(1/alpha)).qTLJ <- function(p, rho, alpha, lower.tail = TRUE, log.p = FALSE) { p <- as.numeric(p) if (!is.finite(rho) || !is.finite(alpha) || rho <= 0 || alpha <= 0) stop("rho and alpha must be positive.") if (log.p) p <- exp(p) if (!lower.tail) p <- 1 - p if (any(!is.finite(p) | p < 0 | p > 1)) stop("p must be in [0,1].") out <- numeric(length(p)) out[p <= 0] <- 0 out[p >= 1] <- Inf ok <- p > 0 & p < 1 if (any(ok)) { out[ok] <- vapply(p[ok], function(pp) { target <- sqrt(pmax(1 - pp^(1/alpha), .Machine$double.xmin)) rootfun <- function(z) exp(tlj_logS_J(z, rho)) - target upper <- max(1, 2/rho) while (rootfun(upper) > 0 && upper < 1e8) upper <- 2*upper if (upper >= 1e8 && rootfun(upper) > 0) stop("The upper quantile bracket could not be found.") uniroot(rootfun, lower = 0, upper = upper, tol = 1e-11)$root }, numeric(1)) } out} rTLJ <- function(n, rho, alpha) { qTLJ(runif(n, 1e-10, 1 - 1e-10), rho, alpha)} ############################################################# MODEL VALIDITY CHECK############################################################ check_TLJ <- function(rho, alpha) { x_split <- qTLJ(0.001, rho, alpha) integral <- pTLJ(x_split, rho, alpha) + integrate(function(z) dTLJ(z, rho, alpha), lower = x_split, upper = Inf, rel.tol = 1e-8, subdivisions = 3000L)$value probs <- c(0.01, 0.10, 0.50, 0.90, 0.99) inv_error <- max(abs(pTLJ(qTLJ(probs, rho, alpha), rho, alpha) - probs)) data.frame(rho = rho, alpha = alpha, F0 = pTLJ(0, rho, alpha), F_Infinity = pTLJ(Inf, rho, alpha), PDF_Integral = integral, Inverse_Max_Error = inv_error, Valid = abs(integral - 1) < 1e-6 && inv_error < 1e-7)} ############################################################# MAXIMUM LIKELIHOOD ESTIMATION############################################################ nll_TLJ <- function(eta, data) { par <- exp(eta) lf <- dTLJ(data, par[1], par[2], log = TRUE) if (any(!is.finite(lf))) return(1e100) -sum(lf)} mle_once <- function(x, init = c(rho = 1, alpha = 1), nstart = MLE_NSTART) { stopifnot(all(is.finite(x)), all(x > 0), all(init > 0)) starts <- matrix(NA_real_, nstart, 2) starts[1, ] <- pmin(pmax(log(init), log(MLE_LOWER)), log(MLE_UPPER)) if (nstart > 1) { for (j in 2:nstart) { starts[j, ] <- pmin(pmax(log(init) + rnorm(2, 0, 0.70), log(MLE_LOWER)), log(MLE_UPPER)) } } fits <- vector("list", nstart) objective <- rep(Inf, nstart) for (j in seq_len(nstart)) { fits[[j]] <- tryCatch( optim(starts[j, ], nll_TLJ, data = x, method = "L-BFGS-B", lower = log(MLE_LOWER), upper = log(MLE_UPPER), control = list(maxit = 4000, factr = 1e7, pgtol = 1e-8)), error = function(e) NULL) if (!is.null(fits[[j]]) && is.finite(fits[[j]]$value)) objective[j] <- fits[[j]]$value } if (all(!is.finite(objective))) stop("No finite MLE solution was obtained.") fit <- fits[[which.min(objective)]] H <- tryCatch(optimHess(fit$par, nll_TLJ, data = x), error = function(e) matrix(NA_real_, 2, 2)) est <- setNames(exp(fit$par), c("rho", "alpha")) eig <- tryCatch(eigen(H, symmetric = TRUE, only.values = TRUE)$values, error = function(e) rep(NA_real_, 2)) eig_min <- if (all(is.finite(eig))) min(eig) else NA_real_ condition <- if (is.finite(eig_min) && eig_min > 0) max(eig) / eig_min else Inf cov_eta <- if (is.finite(eig_min) && eig_min > 1e-8) tryCatch(solve(H), error = function(e) matrix(NA_real_, 2, 2)) else matrix(NA_real_, 2, 2) se_eta <- if (all(is.finite(diag(cov_eta))) && all(diag(cov_eta) >= 0)) sqrt(diag(cov_eta)) else rep(NA_real_, 2) se <- est * se_eta lower <- pmax(est - 1.96 * se, 0) upper <- est + 1.96 * se ci_reliable <- condition <= MLE_HESSIAN_CONDITION_MAX && all(is.finite(se)) if (!ci_reliable) lower[] <- upper[] <- NA_real_ boundary <- est <= 1.001 * MLE_LOWER | est >= 0.999 * MLE_UPPER J <- diag(as.numeric(est), 2) cov_par <- J %*% cov_eta %*% J dimnames(cov_par) <- list(names(est), names(est)) list(est = est, se = se, lower = lower, upper = upper, cov = cov_par, cov_eta = cov_eta, logLik = -fit$value, convergence = fit$convergence, hessian_positive = is.finite(eig_min) && eig_min > 1e-8, hessian_condition = condition, ci_reliable = ci_reliable, at_boundary = boundary)} make_table_MLE <- function(rho, alpha, n_vals = c(25, 50, 100, 200, 500), R = 100) { set.seed(1123) truth <- c(rho = rho, alpha = alpha) ans <- list() for (n in n_vals) { cat("MLE running: n =", n, "\n") est <- lo <- up <- matrix(NA_real_, R, 2, dimnames = list(NULL, names(truth))) conv <- rep(NA_integer_, R) hessian_ok <- ci_ok <- rep(NA, R) hessian_condition <- rep(NA_real_, R) boundary <- matrix(NA, R, 2, dimnames = list(NULL, names(truth))) for (r in seq_len(R)) { x <- rTLJ(n, rho, alpha) fit <- tryCatch(mle_once(x, truth), error = function(e) NULL) if (!is.null(fit)) { est[r, ] <- fit$est; lo[r, ] <- fit$lower; up[r, ] <- fit$upper conv[r] <- fit$convergence; hessian_ok[r] <- fit$hessian_positive ci_ok[r] <- fit$ci_reliable; hessian_condition[r] <- fit$hessian_condition boundary[r, ] <- fit$at_boundary } } T <- matrix(truth, R, 2, byrow = TRUE) mean_est <- colMeans(est, na.rm = TRUE) mse <- colMeans((est - T)^2, na.rm = TRUE) ans[[length(ans) + 1]] <- data.frame( n = n, Parameter = names(truth), True_Value = truth, MLE = mean_est, Bias = mean_est - truth, Relative_Bias = (mean_est - truth) / truth, MSE = mse, RMSE = sqrt(mse), Coverage_95 = colMeans(lo <= T & up >= T, na.rm = TRUE), Lower = apply(lo, 2, median, na.rm = TRUE), Upper = apply(up, 2, median, na.rm = TRUE), Length = apply(up - lo, 2, median, na.rm = TRUE), Successful_Fits = colSums(is.finite(est)), Convergence_Rate = mean(conv == 0, na.rm = TRUE), Positive_Hessian_Rate = mean(hessian_ok, na.rm = TRUE), Reliable_CI_Rate = mean(ci_ok, na.rm = TRUE), Median_Hessian_Condition = median(hessian_condition, na.rm = TRUE), Boundary_Rate = colMeans(boundary, na.rm = TRUE)) } out <- do.call(rbind, ans); rownames(out) <- NULL; out} ######################################################################################################################## # All optimizations are performed on eta = log(rho, alpha), so positivity# is automatic. OLS and WLS use the expected uniform order-statistic# plotting positions p_i = i/(n+1). The WLS weights are the inverse of# Var[U_(i)] = i(n-i+1)/[(n+1)^2(n+2)]. safe_cdf_TLJ <- function(x, par) { pmin(pmax(pTLJ(x, par[1], par[2]), 1e-12), 1 - 1e-12)} objective_MPS <- function(eta, data) { x <- data par <- exp(eta) Fx <- safe_cdf_TLJ(sort(x), par) spacings <- diff(c(0, Fx, 1)) if (any(!is.finite(spacings)) || any(spacings <= 0)) return(1e100) -sum(log(pmax(spacings, 1e-14)))} objective_LS <- function(eta, x, weighted = FALSE) { par <- exp(eta) x <- sort(x); n <- length(x) p <- seq_len(n)/(n + 1) Fx <- safe_cdf_TLJ(x, par) w <- if (weighted) (n + 1)^2 * (n + 2) / (seq_len(n) * (n - seq_len(n) + 1)) else rep(1, n) sum(w * (Fx - p)^2)} # Population L-moments from the quantile representation:# lambda_1 = integral_0^1 Q(u)du,# lambda_2 = integral_0^1 Q(u)(2u-1)du.population_lmoments_TLJ <- function(rho, alpha, nodes = (seq_len(79) - 0.5)/79) { q <- qTLJ(nodes, rho, alpha) c(L1 = mean(q), L2 = mean(q * (2*nodes - 1)))} sample_lmoments <- function(x) { x <- sort(x); n <- length(x) b0 <- mean(x) b1 <- sum((seq_len(n) - 1) * x)/(n * (n - 1)) c(L1 = b0, L2 = 2*b1 - b0)} objective_LM <- function(eta, data) { x <- data par <- exp(eta) target <- sample_lmoments(x) model <- tryCatch(population_lmoments_TLJ(par[1], par[2]), error = function(e) c(NA_real_, NA_real_)) if (any(!is.finite(model))) return(1e100) # Scaling prevents the usually larger first L-moment from dominating L2. scale <- pmax(abs(target), 1e-6) sum(((model - target)/scale)^2)} estimate_TLJ_method <- function(x, method, init = c(rho = 1, alpha = 1), nstart = 5) { method <- match.arg(method, c("MLE", "MPS", "OLS", "WLS", "LM")) fn <- switch(method, MLE = nll_TLJ, MPS = objective_MPS, OLS = function(eta, data) objective_LS(eta, data, FALSE), WLS = function(eta, data) objective_LS(eta, data, TRUE), LM = objective_LM) starts <- matrix(NA_real_, nstart, 2) starts[1, ] <- log(init) if (nstart > 1) { for (j in 2:nstart) starts[j, ] <- log(init) + rnorm(2, 0, 0.45) } starts <- pmin(pmax(starts, log(MLE_LOWER)), log(MLE_UPPER)) fits <- lapply(seq_len(nstart), function(j) tryCatch( optim(starts[j, ], fn, data = x, method = "L-BFGS-B", lower = log(MLE_LOWER), upper = log(MLE_UPPER), control = list(maxit = 2500, factr = 1e8, pgtol = 1e-7)), error = function(e) NULL)) values <- vapply(fits, function(z) if (is.null(z) || !is.finite(z$value)) Inf else z$value, numeric(1)) if (all(!is.finite(values))) return(list(est = c(rho = NA_real_, alpha = NA_real_), convergence = 99L, objective = NA_real_, at_boundary = TRUE)) best <- fits[[which.min(values)]] est <- setNames(exp(best$par), c("rho", "alpha")) at_boundary <- any(est <= 1.001*MLE_LOWER | est >= 0.999*MLE_UPPER) list(est = est, convergence = best$convergence, objective = best$value, at_boundary = at_boundary) run_reviewer1_simulation <- function( rho_values = 0.8, alpha_values = c(0.5, 1.0, 1.5), n_values = c(25, 50, 100, 200, 500), R = 100, methods = c("MLE", "MPS", "OLS", "WLS", "LM"), seed = 20260917) { set.seed(seed) raw <- list(); summary_out <- list(); k <- 0L for (rho in rho_values) for (alpha in alpha_values) for (n in n_values) { truth <- c(rho = rho, alpha = alpha) cat(sprintf("Reviewer 1 simulation: rho=%.3f alpha=%.3f n=%d\n", rho, alpha, n)) estimates <- array(NA_real_, c(R, length(methods), 2), dimnames = list(NULL, methods, names(truth))) elapsed <- matrix(NA_real_, R, length(methods), dimnames = list(NULL, methods)) convergence <- matrix(NA_integer_, R, length(methods), dimnames = list(NULL, methods)) boundary <- matrix(NA, R, length(methods), dimnames = list(NULL, methods)) for (r in seq_len(R)) { x <- rTLJ(n, rho, alpha) for (m in methods) { tm <- system.time({ fit <- estimate_TLJ_method(x, m, init = truth) })[["elapsed"]] estimates[r, m, ] <- fit$est elapsed[r, m] <- tm convergence[r, m] <- fit$convergence boundary[r, m] <- fit$at_boundary } } for (m in methods) { for (p in names(truth)) { z <- estimates[, m, p] ok <- is.finite(z) k <- k + 1L raw[[k]] <- data.frame( rho_true = rho, alpha_true = alpha, n = n, Replication = seq_len(R), Method = m, Parameter = p, Estimate = z, CPU_seconds = elapsed[, m], Convergence = convergence[, m], Boundary = boundary[, m]) summary_out[[k]] <- data.frame( rho_true = rho, alpha_true = alpha, Alpha_Regime = if (alpha < 1) "alpha < 1" else if (alpha == 1) "alpha = 1" else "alpha > 1", n = n, Method = m, Parameter = p, True_Value = truth[p], Mean_Estimate = if (any(ok)) mean(z[ok]) else NA_real_, MRB = if (any(ok)) mean((z[ok] - truth[p])/truth[p]) else NA_real_, RMSE = if (any(ok)) sqrt(mean((z[ok] - truth[p])^2)) else NA_real_, Mean_CPU_seconds = mean(elapsed[, m], na.rm = TRUE), Median_CPU_seconds = median(elapsed[, m], na.rm = TRUE), Successful_Fits = sum(ok), Success_Rate = mean(ok), Convergence_Rate = mean(convergence[, m] == 0, na.rm = TRUE), Boundary_Rate = mean(boundary[, m], na.rm = TRUE)) } } } list(summary = do.call(rbind, summary_out), raw = do.call(rbind, raw))} ######################################################################################################################## fd_gradient <- function(fn, theta, rel_step = 1e-5) { h <- rel_step * pmax(abs(theta), 1) out <- numeric(length(theta)) for (j in seq_along(theta)) { tp <- tm <- theta; tp[j] <- tp[j] + h[j]; tm[j] <- tm[j] - h[j] if (tm[j] <= 0) tm[j] <- theta[j]/2 out[j] <- (fn(tp) - fn(tm))/(tp[j] - tm[j]) } out} fd_hessian <- function(fn, theta, rel_step = 2e-4) { p <- length(theta); H <- matrix(NA_real_, p, p) h <- rel_step * pmax(abs(theta), 1); f0 <- fn(theta) for (i in seq_len(p)) { tp <- tm <- theta; tp[i] <- tp[i] + h[i]; tm[i] <- tm[i] - h[i] if (tm[i] <= 0) tm[i] <- theta[i]/2 hi <- (tp[i] - tm[i])/2 H[i,i] <- (fn(tp) - 2*f0 + fn(tm))/hi^2 if (i < p) for (j in (i+1):p) { hj <- h[j]; pp <- pm <- mp <- mm <- theta pp[c(i,j)] <- theta[c(i,j)] + c(hi,hj) pm[c(i,j)] <- theta[c(i,j)] + c(hi,-hj) mp[c(i,j)] <- theta[c(i,j)] + c(-hi,hj) mm[c(i,j)] <- theta[c(i,j)] - c(hi,hj) if (any(c(pm,mp,mm) <= 0)) next H[i,j] <- H[j,i] <- (fn(pp)-fn(pm)-fn(mp)+fn(mm))/(4*hi*hj) } } H} fisher_comparison <- function(rho_grid = c(.05,.20,.80,3), alpha_grid = c(.05,.20,1,3), n_observed = 500, mc_expected = FISHER_MC_SIZE, seed = 8841) { set.seed(seed); out <- list(); k <- 0L for (rho in rho_grid) for (alpha in alpha_grid) { th <- c(rho, alpha); x <- rTLJ(n_observed, rho, alpha) nll_nat <- function(z) -sum(dTLJ(x,z[1],z[2],log=TRUE)) Iobs <- fd_hessian(nll_nat, th)/n_observed # Deterministic midpoint quadrature on U(0,1), X=Q(U), for E[ss^T]. xm <- qTLJ((seq_len(mc_expected)-.5)/mc_expected, rho, alpha) scores <- t(vapply(xm, function(xx) fd_gradient( function(z) dTLJ(xx,z[1],z[2],log=TRUE), th), numeric(2))) Iexp <- crossprod(scores)/nrow(scores) eo <- eigen(Iobs,symmetric=TRUE,only.values=TRUE)$values ee <- eigen(Iexp,symmetric=TRUE,only.values=TRUE)$values k <- k+1L; out[[k]] <- data.frame( rho=rho,alpha=alpha, Obs_Eigen_Min=min(eo),Obs_Eigen_Max=max(eo), Obs_Condition=if(min(eo)>0) max(eo)/min(eo) else Inf, Exp_Eigen_Min=min(ee),Exp_Eigen_Max=max(ee), Exp_Condition=if(min(ee)>0) max(ee)/min(ee) else Inf, Frobenius_Difference=sqrt(sum((Iobs-Iexp)^2))) } do.call(rbind,out)} ######################################################################################################################## gof_statistics_TLJ <- function(x, par) { u <- sort(pmin(pmax(pTLJ(sort(x),par[1],par[2]),1e-12),1-1e-12)) n <- length(u); i <- seq_len(n) KS <- max(max(i/n-u), max(u-(i-1)/n)) CvM <- 1/(12*n)+sum((u-(2*i-1)/(2*n))^2) AD <- -n-mean((2*i-1)*(log(u)+log1p(-rev(u)))) c(KS=KS, AD=AD, CvM=CvM)} parametric_bootstrap_gof_TLJ <- function(x, B=BOOTSTRAP_GOF_B, seed=9017, prefix="TLJ") { set.seed(seed); fit <- mle_once(x); obs <- gof_statistics_TLJ(x,fit$est) boot <- matrix(NA_real_,B,3,dimnames=list(NULL,names(obs))) for (b in seq_len(B)) { xb <- rTLJ(length(x),fit$est[1],fit$est[2]) fb <- tryCatch(mle_once(xb,fit$est,nstart=4),error=function(e) NULL) if (!is.null(fb)) boot[b,] <- gof_statistics_TLJ(xb,fb$est) if (b %% 250 == 0) cat("GoF bootstrap",b,"of",B,"\n") } pval <- (1+colSums(sweep(boot,2,obs,">="),na.rm=TRUE))/(1+colSums(is.finite(boot))) png(paste0(prefix,"_Bootstrap_AD_CvM.png"),1200,500,res=150) par(mfrow=c(1,2)); hist(boot[,"AD"],breaks=40,col="grey85",border="grey40", main="Bootstrap distribution of AD",xlab="AD"); abline(v=obs["AD"],lwd=2) hist(boot[,"CvM"],breaks=40,col="grey85",border="grey40", main="Bootstrap distribution of CvM",xlab="CvM"); abline(v=obs["CvM"],lwd=2) dev.off() list(summary=data.frame(Statistic=names(obs),Observed=obs, Bootstrap_p_value=pval,Successful=colSums(is.finite(boot))), bootstrap=boot, fit=fit)} ######################################################################################################################## stress_strength_R <- function(theta, nodes=(seq_len(399)-.5)/399) { xq <- qTLJ(nodes,theta[1],theta[2]) mean(pTLJ(xq,theta[3],theta[4]))} estimate_stress_strength <- function(x,y) { fx <- mle_once(x); fy <- mle_once(y) theta <- c(fx$est,fy$est); names(theta) <- c("rhoX","alphaX","rhoY","alphaY") Rhat <- stress_strength_R(theta) g <- fd_gradient(stress_strength_R,theta) V <- matrix(0,4,4); V[1:2,1:2] <- fx$cov; V[3:4,3:4] <- fy$cov se <- sqrt(drop(t(g)%*%V%*%g)) c(R=Rhat,SE=se,Lower=max(0,Rhat-1.96*se),Upper=min(1,Rhat+1.96*se))} stress_strength_bootstrap_t <- function(x,y,B=STRESS_STRENGTH_BOOT_B,seed=7419) { set.seed(seed); base <- estimate_stress_strength(x,y); fx<-mle_once(x); fy<-mle_once(y) tv <- numeric(B) for (b in seq_len(B)) { xb<-rTLJ(length(x),fx$est[1],fx$est[2]); yb<-rTLJ(length(y),fy$est[1],fy$est[2]) eb<-tryCatch(estimate_stress_strength(xb,yb),error=function(e) rep(NA_real_,4)) if (is.finite(eb[2]) && eb[2]>0) tv[b]<-(eb[1]-base[1])/eb[2] else tv[b]<-NA_real_ } q<-quantile(tv,c(.975,.025),na.rm=TRUE) c(base,Studentized_Lower=max(0,base[1]-q[1]*base[2]), Studentized_Upper=min(1,base[1]-q[2]*base[2]))} stress_strength_coverage <- function(theta=c(0.8,.7,.8,1.5),n=100, R=STRESS_STRENGTH_MC_R,B=STRESS_STRENGTH_BOOT_B) { trueR<-stress_strength_R(theta); delta<-boot<-logical(R) for (r in seq_len(R)) { x<-rTLJ(n,theta[1],theta[2]); y<-rTLJ(n,theta[3],theta[4]) d<-tryCatch(estimate_stress_strength(x,y),error=function(e) rep(NA_real_,4)) s<-tryCatch(stress_strength_bootstrap_t(x,y,B,seed=8000+r),error=function(e) rep(NA_real_,6)) delta[r]<-is.finite(d[3]) && d[3]<=trueR && trueR<=d[4] boot[r]<-is.finite(s[5]) && s[5]<=trueR && trueR<=s[6] } data.frame(n=n,Replications=R,Bootstrap_B=B,True_R=trueR, Delta_Coverage=mean(delta),Studentized_Bootstrap_Coverage=mean(boot))} ############################################################# BAYES: RANDOM-WALK METROPOLIS-HASTINGS ON LOG SCALE############################################################ logpost_TLJ <- function(eta, x, prior_shape = PRIOR_SHAPE, prior_rate = PRIOR_RATE) { par <- exp(eta); names(par) <- c("rho", "alpha") lf <- dTLJ(x, par[1], par[2], log = TRUE) if (any(!is.finite(lf))) return(-Inf) lp <- sum(dgamma(par, shape = prior_shape, rate = prior_rate, log = TRUE)) sum(lf) + lp + sum(eta)} bayes_once <- function(x, niter = BAYES_NITER, burn = BAYES_BURN, thin = BAYES_THIN, start = c(rho = 1, alpha = 1), prop_sd = PROP_SD_MAIN, prior_shape = PRIOR_SHAPE, prior_rate = PRIOR_RATE, return_draws = FALSE) { if (burn >= niter) stop("burn must be smaller than niter.") current <- log(start) lp_current <- logpost_TLJ(current, x, prior_shape, prior_rate) keep <- seq(burn + thin, niter, by = thin) draws <- matrix(NA_real_, length(keep), 2, dimnames = list(NULL, c("rho", "alpha"))) accepted_post <- saved <- 0 for (it in seq_len(niter)) { proposed <- rnorm(2, current, prop_sd) lp_proposed <- logpost_TLJ(proposed, x, prior_shape, prior_rate) accepted_now <- is.finite(lp_proposed) && log(runif(1)) < lp_proposed - lp_current if (accepted_now) { current <- proposed; lp_current <- lp_proposed if (it > burn) accepted_post <- accepted_post + 1 } if (it > burn && (it - burn) %% thin == 0) { saved <- saved + 1; draws[saved, ] <- exp(current) } } draws <- draws[complete.cases(draws), , drop = FALSE] HPD <- HPDinterval(as.mcmc(draws), prob = 0.95) sm <- cbind(Mean = colMeans(draws), Median = apply(draws, 2, median), SD = apply(draws, 2, sd), Q025 = apply(draws, 2, quantile, 0.025), Q975 = apply(draws, 2, quantile, 0.975), HPD_Lower = HPD[, 1], HPD_Upper = HPD[, 2]) out <- list(summary = sm, acceptance_rate = accepted_post / (niter - burn)) if (return_draws) out$draws <- draws out} make_table_BAYES <- function(rho, alpha, n_vals = c(25, 50, 100, 200, 500), R = 100, niter = BAYES_NITER, burn = BAYES_BURN, thin = BAYES_THIN) { set.seed(1456) truth <- c(rho = rho, alpha = alpha) ans <- list() for (n in n_vals) { cat("Bayes running: n =", n, "\n") est <- etl <- etu <- hl <- hu <- matrix(NA_real_, R, 2, dimnames = list(NULL, names(truth))) acc <- rep(NA_real_, R) for (r in seq_len(R)) { if (r %% 10 == 0) cat(" Replication", r, "of", R, "\n") x <- rTLJ(n, rho, alpha) fit <- tryCatch(bayes_once(x, niter, burn, thin, truth), error = function(e) NULL) if (!is.null(fit)) { est[r, ] <- fit$summary[, "Mean"] etl[r, ] <- fit$summary[, "Q025"]; etu[r, ] <- fit$summary[, "Q975"] hl[r, ] <- fit$summary[, "HPD_Lower"]; hu[r, ] <- fit$summary[, "HPD_Upper"] acc[r] <- fit$acceptance_rate } } T <- matrix(truth, R, 2, byrow = TRUE) mean_est <- colMeans(est, na.rm = TRUE) mse <- colMeans((est - T)^2, na.rm = TRUE) ans[[length(ans) + 1]] <- data.frame( n = n, Parameter = names(truth), True_Value = truth, Bayes = mean_est, Bias = mean_est - truth, Relative_Bias = (mean_est - truth) / truth, MSE = mse, RMSE = sqrt(mse), Equal_Tail_Coverage_95 = colMeans(etl <= T & etu >= T, na.rm = TRUE), HPD_Coverage_95 = colMeans(hl <= T & hu >= T, na.rm = TRUE), Equal_Tail_Lower = colMeans(etl, na.rm = TRUE), Equal_Tail_Upper = colMeans(etu, na.rm = TRUE), Equal_Tail_Length = colMeans(etu - etl, na.rm = TRUE), HPD_Lower = colMeans(hl, na.rm = TRUE), HPD_Upper = colMeans(hu, na.rm = TRUE), HPD_Length = colMeans(hu - hl, na.rm = TRUE), Acceptance_Rate = mean(acc, na.rm = TRUE), Successful_Fits = colSums(is.finite(est))) } out <- do.call(rbind, ans); rownames(out) <- NULL; out} ############################################################# MCMC DIAGNOSTICS############################################################ MCSE_coda <- function(z) sqrt(spectrum0.ar(as.numeric(z))$spec / length(z)) compute_mcmc_diagnostics <- function(draws, acceptance_rate) { ess <- effectiveSize(as.mcmc(draws)) mcse <- apply(draws, 2, MCSE_coda) data.frame(Parameter = colnames(draws), Posterior_Mean = colMeans(draws), Posterior_SD = apply(draws, 2, sd), ESS = as.numeric(ess), MCSE = mcse, MCSE_SD_Ratio = mcse / apply(draws, 2, sd), Acceptance_Rate = acceptance_rate)} compute_rhat <- function(chain_list) { g <- gelman.diag(mcmc.list(lapply(chain_list, function(z) mcmc(z$draws))), autoburnin = FALSE, multivariate = FALSE)$psrf data.frame(Parameter = rownames(g), Rhat = g[, 1], Upper_CI = g[, 2])} prior_sensitivity <- function(x, start) { priors <- list( Main_Gamma_2_1 = list(shape = c(rho = 2, alpha = 2), rate = c(rho = 1, alpha = 1)), Weak_Gamma_1_0.5 = list(shape = c(rho = 1, alpha = 1), rate = c(rho = 0.5, alpha = 0.5)), Diffuse_Gamma_0.5_0.1 = list(shape = c(rho = 0.5, alpha = 0.5), rate = c(rho = 0.1, alpha = 0.1))) out <- list() for (nm in names(priors)) { pr <- priors[[nm]] z <- bayes_once(x, start = start, prior_shape = pr$shape, prior_rate = pr$rate, return_draws = TRUE) out[[nm]] <- data.frame( Prior = nm, Parameter = rownames(z$summary), Posterior_Mean = z$summary[, "Mean"], Posterior_SD = z$summary[, "SD"], HPD_Lower = z$summary[, "HPD_Lower"], HPD_Upper = z$summary[, "HPD_Upper"], Acceptance_Rate = z$acceptance_rate) } ans <- do.call(rbind, out); rownames(ans) <- NULL; ans} ############################################################# REVIEWER 4: PRIOR ROBUSTNESS, FUNCTIONAL HPD AND WAIC############################################################ log_prior_spec <- function(par, spec) { if (spec$type == "gamma") return(sum(dgamma(par,spec$shape,spec$rate,log=TRUE))) if (spec$type == "lognormal") return(sum(dlnorm(par,spec$meanlog,spec$sdlog,log=TRUE))) if (spec$type == "scale_reference") return(-sum(log(par))) stop("Unknown prior type")} bayes_custom_prior <- function(x,spec,start,niter=BAYES_NITER,burn=BAYES_BURN, thin=BAYES_THIN,prop_sd=PROP_SD_MAIN) { logpost <- function(eta) { par<-exp(eta); lf<-dTLJ(x,par[1],par[2],log=TRUE) if(any(!is.finite(lf))) return(-Inf) sum(lf)+log_prior_spec(par,spec)+sum(eta) } cur<-log(start); lp<-logpost(cur); keep<-seq(burn+thin,niter,by=thin) dr<-matrix(NA_real_,length(keep),2,dimnames=list(NULL,c("rho","alpha"))) acc<-saved<-0L for(it in seq_len(niter)) { pr<-rnorm(2,cur,prop_sd); lpp<-logpost(pr) if(is.finite(lpp)&&log(runif(1))<lpp-lp){cur<-pr;lp<-lpp;if(it>burn)acc<-acc+1L} if(it>burn&&(it-burn)%%thin==0){saved<-saved+1L;dr[saved,]<-exp(cur)} } list(draws=dr,acceptance=acc/(niter-burn))} posterior_WAIC <- function(x,draws) { ll<-vapply(seq_len(nrow(draws)),function(i) dTLJ(x,draws[i,1],draws[i,2],log=TRUE),numeric(length(x))) lppd<-sum(log(rowMeans(exp(sweep(ll,1,apply(ll,1,max),"-"))))+apply(ll,1,max)) pwaic<-sum(apply(ll,1,var)); c(WAIC=-2*(lppd-pwaic),lppd=lppd,p_WAIC=pwaic)} functional_HPD <- function(draws,mission_times=c(.5,1,2)) { out<-list(); k<-0L for(t0 in mission_times) for(fun in c("Reliability","Hazard")) { z<-if(fun=="Reliability") mapply(sTLJ,t0,draws[,1],draws[,2]) else mapply(hTLJ,t0,draws[,1],draws[,2]) hp<-HPDinterval(as.mcmc(z),prob=.95) k<-k+1L; out[[k]]<-data.frame(Function=fun,t=t0,Mean=mean(z), SD=sd(z),HPD_Lower=hp[1],HPD_Upper=hp[2]) } do.call(rbind,out)} prior_robustness_extended <- function(x,start,mission_times=c(.5,1,2)) { specs<-list( Scale_reference=list(type="scale_reference"), Gamma_informative=list(type="gamma",shape=c(8,8),rate=c(8,8)), Gamma_moderate=list(type="gamma",shape=c(2,2),rate=c(1,1)), Gamma_ultradiffuse=list(type="gamma",shape=c(.5,.5),rate=c(.1,.1)), Lognormal_informative=list(type="lognormal",meanlog=log(c(1,1)),sdlog=c(.25,.25)), Lognormal_moderate=list(type="lognormal",meanlog=log(c(1,1)),sdlog=c(1,1)), Lognormal_ultradiffuse=list(type="lognormal",meanlog=log(c(1,1)),sdlog=c(2.5,2.5)) ) pars<-funcs<-list() for(nm in names(specs)) { fit<-bayes_custom_prior(x,specs[[nm]],start) hp<-HPDinterval(as.mcmc(fit$draws),prob=.95); w<-posterior_WAIC(x,fit$draws) pars[[nm]]<-data.frame(Prior=nm,Parameter=colnames(fit$draws), Mean=colMeans(fit$draws),SD=apply(fit$draws,2,sd),HPD_Lower=hp[,1], HPD_Upper=hp[,2],WAIC=w["WAIC"],Acceptance=fit$acceptance) ff<-functional_HPD(fit$draws,mission_times); ff$Prior<-nm; funcs[[nm]]<-ff } list(parameters=do.call(rbind,pars),functionals=do.call(rbind,funcs))} ############################################################# COLORED GRAPHS############################################################ COL_MLE <- "#0072B2"; COL_BAYES <- "#D55E00"; COL_TRUE <- "#009E73"COL_HPD <- "#CC79A7"; COL_FILL <- "#56B4E9" plot_trace_hist <- function(draws, true_values, file_name = "TLJ_Scenario2_MCMC_trace_hist.png") { png(file_name, 1500, 800, res = 150) old <- par(no.readonly = TRUE) on.exit({par(old); dev.off()}) par(mfrow = c(2, 2), mar = c(4, 4, 3, 1)) for (p in colnames(draws)) { z <- draws[, p]; mn <- mean(z); ci <- quantile(z, c(.025, .975)) plot(z, type = "l", col = COL_MLE, lwd = .8, xlab = "Iteration", ylab = p, main = paste("Trace plot:", p)) abline(h = mn, col = COL_BAYES, lwd = 2) abline(h = ci, col = COL_HPD, lwd = 2, lty = 2) abline(h = true_values[p], col = COL_TRUE, lwd = 2, lty = 3) hist(z, breaks = 40, probability = TRUE, col = adjustcolor(COL_FILL, .65), border = "white", xlab = p, main = paste("Posterior histogram:", p)) lines(density(z), col = COL_MLE, lwd = 2) abline(v = mn, col = COL_BAYES, lwd = 2) abline(v = ci, col = COL_HPD, lwd = 2, lty = 2) abline(v = true_values[p], col = COL_TRUE, lwd = 2, lty = 3) }} plot_compare <- function(mle, bayes, p, metric = "MSE", true_values = NULL) { m <- subset(mle, Parameter == p); b <- subset(bayes, Parameter == p) if (metric == "MSE") { y1 <- m$MSE; y2 <- b$MSE; ylab <- paste("MSE of", p); ttl <- paste("MSE:", p) } else { y1 <- m$MLE; y2 <- b$Bayes; ylab <- paste("Estimate of", p); ttl <- paste("Estimates:", p) } ylim <- range(c(y1, y2, if (metric == "Estimate") true_values[p]), na.rm = TRUE) plot(m$n, y1, type = "b", pch = 16, lwd = 2, col = COL_MLE, ylim = ylim, xlab = "Sample size", ylab = ylab, main = ttl) lines(b$n, y2, type = "b", pch = 17, lwd = 2, col = COL_BAYES) if (metric == "Estimate") abline(h = true_values[p], col = COL_TRUE, lwd = 2, lty = 2) legend("topright", c("MLE", "Bayes", if (metric == "Estimate") "True value"), col = c(COL_MLE, COL_BAYES, if (metric == "Estimate") COL_TRUE), pch = c(16, 17, if (metric == "Estimate") NA), lty = c(1, 1, if (metric == "Estimate") 2), bty = "n") grid(col = "grey85")} plot_interval <- function(mle, bayes, p) { m <- subset(mle, Parameter == p); b <- subset(bayes, Parameter == p) ylim <- range(c(m$Length, b$Equal_Tail_Length, b$HPD_Length), na.rm = TRUE) plot(m$n, m$Length, type = "b", pch = 16, lwd = 2, col = COL_MLE, ylim = ylim, xlab = "Sample size", ylab = paste("Interval length of", p), main = paste("Interval lengths:", p)) lines(b$n, b$Equal_Tail_Length, type = "b", pch = 17, lwd = 2, col = COL_BAYES) lines(b$n, b$HPD_Length, type = "b", pch = 18, lwd = 2, lty = 2, col = COL_HPD) legend("topright", c("MLE CI", "Bayes equal-tail", "Bayes HPD"), col = c(COL_MLE, COL_BAYES, COL_HPD), pch = c(16, 17, 18), lty = c(1, 1, 2), bty = "n") grid(col = "grey85")} plot_pdf_cdf_hazard <- function(rho, alpha, file_name = "TLJ_Scenario2_PDF_CDF_Hazard.png") { xmax <- qTLJ(.995, rho, alpha) x <- seq(max(xmax * 1e-5, 1e-7), xmax, length.out = 600) png(file_name, 1500, 500, res = 150) par(mfrow = c(1, 3), mar = c(4, 4, 3, 1)) plot(x, dTLJ(x, rho, alpha), type = "l", lwd = 3, col = "#E69F00", xlab = "x", ylab = "f(x)", main = "TLJ density"); grid(col = "grey88") plot(x, pTLJ(x, rho, alpha), type = "l", lwd = 3, col = "#0072B2", xlab = "x", ylab = "F(x)", main = "TLJ distribution"); grid(col = "grey88") plot(x, hTLJ(x, rho, alpha), type = "l", lwd = 3, col = "#D55E00", xlab = "x", ylab = "h(x)", main = "TLJ hazard"); grid(col = "grey88") dev.off()} ############################################################# RUN THE COMPLETE SIMULATION############################################################ rho_true <- 0.8alpha_true <- 1.5true_values <- c(rho = rho_true, alpha = alpha_true)n_vals_demo <- c(25, 50, 100, 200, 500)R_sim <- MC_REPLICATIONSR_REVIEWER1 <- MC_REPLICATIONS # Save the complete computational design so that every reported result is# reproducible and the reviewer can verify all simulation and MCMC settings.simulation_settings <- data.frame( Setting = c( "R_version", "Monte_Carlo_replications", "Sample_sizes", "Bayes_chains_for_diagnostics", "Bayes_iterations_per_chain", "Bayes_burn_in", "Bayes_thinning", "Saved_draws_per_chain", "Proposal_distribution", "Proposal_SD_rho", "Proposal_SD_alpha", "Prior_rho", "Prior_alpha", "Stress_strength_outer_replications", "Stress_strength_bootstrap_B", "Parametric_GoF_bootstrap_B", "MCMC_package" ), Value = c( paste(R.version$major, R.version$minor, sep = "."), R_sim, paste(n_vals_demo, collapse = ", "), MCMC_NCHAINS, MCMC_DIAG_NITER, MCMC_DIAG_BURN, MCMC_DIAG_THIN, (MCMC_DIAG_NITER - MCMC_DIAG_BURN) / MCMC_DIAG_THIN, "Bivariate normal random walk on the log-parameter scale", PROP_SD_MAIN["rho"], PROP_SD_MAIN["alpha"], "Gamma(shape=2, rate=1)", "Gamma(shape=2, rate=1)", STRESS_STRENGTH_MC_R, STRESS_STRENGTH_BOOT_B, BOOTSTRAP_GOF_B, paste0("coda ", as.character(packageVersion("coda"))) ), stringsAsFactors = FALSE)write.csv(simulation_settings, "TLJ_Computational_Settings.csv", row.names = FALSE)writeLines(capture.output(sessionInfo()), "TLJ_Session_Info.txt")cat("\n=== COMPUTATIONAL SETTINGS ===\n")print(simulation_settings, row.names = FALSE) cat("\n=== TLJ SCENARIO 2 MODEL VALIDITY CHECK ===\n")validity <- check_TLJ(rho_true, alpha_true)print(round_df(validity, 8))write.csv(validity, "TLJ_Scenario2_Distribution_Validity.csv", row.names = FALSE)plot_pdf_cdf_hazard(rho_true, alpha_true) # IMPORTANT SCENARIO DISTINCTION:# The primary true-parameter configuration of Scenario 2 is rho = 0.8 and# alpha = 1.5. All results labelled "Scenario 2" below are generated only# from this configuration for n = {25, 50, 100, 200, 500}.## The values alpha = {0.5, 1.0, 1.5} are examined separately as an additional# robustness and estimator-comparison experiment requested by Reviewer 1.# This supplementary experiment represents the regimes alpha < 1, alpha = 1,# and alpha > 1. MLE, MPS, OLS, WLS and L-moment estimators are compared.reviewer1_results <- run_reviewer1_simulation( rho_values = rho_true, alpha_values = c(0.5, 1.0, 1.5), n_values = c(25, 50, 100, 200, 500), R = R_REVIEWER1, methods = c("MLE", "MPS", "OLS", "WLS", "LM"))write.csv(reviewer1_results$summary, "TLJ_Additional_Robustness_Estimator_Comparison_Summary.csv", row.names = FALSE)write.csv(reviewer1_results$raw, "TLJ_Additional_Robustness_Estimator_Comparison_Raw.csv", row.names = FALSE)cat("\n=== ADDITIONAL ROBUSTNESS AND ESTIMATOR COMPARISON ===\n")print(round_df(reviewer1_results$summary)) if (RUN_FISHER_STABILITY) { fisher_table <- fisher_comparison(mc_expected=FISHER_MC_SIZE) write.csv(fisher_table,"TLJ_Observed_Expected_Fisher_Stability.csv",row.names=FALSE) cat("\n=== OBSERVED VERSUS EXPECTED FISHER INFORMATION ===\n") print(round_df(fisher_table))} table_mle <- make_table_MLE(rho_true, alpha_true, n_vals_demo, R_sim)write.csv(table_mle, "TLJ_Scenario2_MLE_results.csv", row.names = FALSE)cat("\n=== TLJ SCENARIO 2 MLE TABLE ===\n"); print(round_df(table_mle)) table_bayes <- make_table_BAYES(rho_true, alpha_true, n_vals_demo, R_sim)write.csv(table_bayes, "TLJ_Scenario2_Bayes_results.csv", row.names = FALSE)cat("\n=== TLJ SCENARIO 2 BAYES TABLE ===\n"); print(round_df(table_bayes)) combined_table <- merge(table_mle, table_bayes, by = c("n", "Parameter", "True_Value"), suffixes = c("_MLE", "_Bayes"))combined_table <- combined_table[order(combined_table$n, combined_table$Parameter), ]write.csv(combined_table, "TLJ_Scenario2_MLE_Bayes_combined.csv", row.names = FALSE) for (metric in c("MSE", "Estimate")) { png(paste0("TLJ_Scenario2_", metric, "_MLE_Bayes.png"), 1200, 500, res = 150) par(mfrow = c(1, 2)) for (p in names(true_values)) plot_compare(table_mle, table_bayes, p, metric, true_values) dev.off()} png("TLJ_Scenario2_Interval_Lengths.png", 1200, 500, res = 150)par(mfrow = c(1, 2))for (p in names(true_values)) plot_interval(table_mle, table_bayes, p)dev.off() ############################################################# DIAGNOSTICS, R-HAT AND PRIOR SENSITIVITY############################################################ set.seed(1789)x_demo <- rTLJ(300, rho_true, alpha_true)diagnostic_fit <- bayes_once(x_demo, start = true_values, return_draws = TRUE)plot_trace_hist(diagnostic_fit$draws, true_values) diagnostic_table <- compute_mcmc_diagnostics( diagnostic_fit$draws, diagnostic_fit$acceptance_rate)write.csv(diagnostic_table, "TLJ_Scenario2_MCMC_Diagnostics.csv", row.names = FALSE)cat("\n=== TLJ SCENARIO 2 MCMC DIAGNOSTICS ===\n"); print(round_df(diagnostic_table)) stopifnot(length(MCMC_STARTS) == MCMC_NCHAINS)chains <- lapply(MCMC_STARTS, function(s) bayes_once(x_demo, niter = MCMC_DIAG_NITER, burn = MCMC_DIAG_BURN, thin = MCMC_DIAG_THIN, start = s, return_draws = TRUE))rhat_table <- compute_rhat(chains)write.csv(rhat_table, "TLJ_Scenario2_Rhat.csv", row.names = FALSE)cat("\n=== TLJ SCENARIO 2 GELMAN-RUBIN R-HAT ===\n"); print(round_df(rhat_table)) prior_table <- prior_sensitivity(x_demo, true_values)write.csv(prior_table, "TLJ_Scenario2_Prior_Sensitivity.csv", row.names = FALSE)cat("\n=== TLJ SCENARIO 2 PRIOR SENSITIVITY ===\n"); print(round_df(prior_table)) if (RUN_BAYES_ROBUSTNESS) { robust_bayes <- prior_robustness_extended(x_demo,true_values,c(.5,1,2)) write.csv(robust_bayes$parameters,"TLJ_Extended_Prior_Robustness_WAIC.csv",row.names=FALSE) write.csv(robust_bayes$functionals,"TLJ_Functional_HPD_Reliability_Hazard.csv",row.names=FALSE)} if (RUN_STRESS_STRENGTH) { ss_table <- stress_strength_coverage( theta=c(0.8,.7,.8,1.5),n=100, R=STRESS_STRENGTH_MC_R,B=STRESS_STRENGTH_BOOT_B) write.csv(ss_table,"TLJ_Stress_Strength_Delta_Bootstrap_Coverage.csv",row.names=FALSE) print(round_df(ss_table))}