X <- c(53,55,85,87,22,23,25,93,51,62,53,32,43,47,30,88,59,49,42,71,41,82,75,93,37)
Y <- c(89,90,59,50,25,29,54,62,39,25,89,32,33,63,38,77,55,41,31,66,57,32,43,88,34)
nX <- length(X); nY <- length(Y)
set.seed(12345)

aic_val <- function(ll,k) -2*ll + 2*k
bic_val <- function(ll,k,n) -2*ll + log(n)*k
ks_fit  <- function(x, Fcdf) suppressWarnings(ks.test(x, Fcdf, exact=FALSE))

we_pdf <- function(x, a, l) ((a+1)/a)*l*exp(-l*x)*(1-exp(-a*l*x))
we_cdf <- function(x, a, l) 1-((a+1)/a)*exp(-l*x)+(1/a)*exp(-(a+1)*l*x)

nll_we <- function(p, x){ a <- exp(p[1]); l <- exp(p[2])
fx <- we_pdf(x,a,l); if(any(!is.finite(fx)|fx<=0)) return(1e12); -sum(log(fx)) }
fit_we_ms <- function(x){
  lam0 <- c(1/mean(x), 2/mean(x), 0.5/mean(x))
  a0   <- c(0.7,1.0,1.5,2.0)
  best <- Inf; par_best <- NULL
  for(l0 in lam0) for(a in a0){
    o <- optim(log(c(a,l0)), nll_we, x=x, method="L-BFGS-B",
               control=list(reltol=1e-10, maxit=2000))
    if(o$value < best){ best <- o$value; par_best <- o$par }
  }
  a <- exp(par_best[1]); l <- exp(par_best[2])
  list(alpha=a, lambda=l, loglik=-best)
}

nll_bwe_marg <- function(p, x, y){
  l1 <- exp(p[1]); l2 <- exp(p[2]); l3 <- exp(p[3])
  aX <- (l2+l3)/l1; aY <- (l1+l3)/l2
  fx <- we_pdf(x,aX,l1); fy <- we_pdf(y,aY,l2)
  if(any(!is.finite(fx)|fx<=0) || any(!is.finite(fy)|fy<=0)) return(1e12)
  -(sum(log(fx))+sum(log(fy)))
}
fit_bwe_joint <- function(x,y){
  l1_0 <- 1/mean(x); l2_0 <- 1/mean(y); l3_0 <- 0.5*min(l1_0,l2_0)
  starts <- rbind(c(l1_0,l2_0,l3_0),
                  c(2*l1_0,l2_0,0.3*l3_0+1e-3),
                  c(l1_0,2*l2_0,0.8*l3_0+1e-3),
                  c(0.5*l1_0,l2_0,1.2*l3_0+1e-3))
  best <- Inf; par_best <- NULL
  for(i in 1:nrow(starts)){
    o <- optim(log(starts[i,]), nll_bwe_marg, x=x, y=y, method="L-BFGS-B",
               control=list(reltol=1e-10, maxit=4000))
    if(o$value < best){ best <- o$value; par_best <- o$par }
  }
  l1 <- exp(par_best[1]); l2 <- exp(par_best[2]); l3 <- exp(par_best[3])
  aX <- (l2+l3)/l1; aY <- (l1+l3)/l2
  llX <- sum(log(we_pdf(x,aX,l1))); llY <- sum(log(we_pdf(y,aY,l2)))
  list(lambda1=l1, lambda2=l2, lambda3=l3,
       alphaX=aX, lambdaX=l1, loglik_X=llX,
       alphaY=aY, lambdaY=l2, loglik_Y=llY)
}

bwrr2_pdf <- function(x, xi, b, k){ t <- (xi*x^2)/(2*b^2); (x*xi*k/b^2)*t^(k-1)*exp(-t^k) }
bwrr2_cdf <- function(x, xi, b, k){ t <- (xi*x^2)/(2*b^2); 1-exp(-t^k) }
nll_wr <- function(p,x){ xi<-exp(p[1]); b<-exp(p[2]); k<-exp(p[3])
fx <- bwrr2_pdf(x,xi,b,k); if(any(!is.finite(fx)|fx<=0)) return(1e12); -sum(log(fx)) }
fit_wr <- function(x){
  starts <- rbind(c(1/mean(x), sd(x)/sqrt(2), 1.0),
                  c(2/mean(x), median(x)/sqrt(2), 0.8),
                  c(0.5/mean(x), IQR(x)/sqrt(2), 1.2))
  best <- Inf; par_best <- NULL
  for(i in 1:nrow(starts)){
    o <- optim(log(starts[i,]), nll_wr, x=x, method="L-BFGS-B",
               control=list(reltol=1e-10, maxit=2000))
    if(o$value < best){ best <- o$value; par_best <- o$par }
  }
  xi<-exp(par_best[1]); b<-exp(par_best[2]); k<-exp(par_best[3])
  list(xi=xi, b=b, k=k, loglik=-best)
}

bwe <- fit_bwe_joint(X,Y)

wrX <- fit_wr(X); wrY <- fit_wr(Y)

ks_we_X <- ks_fit(X, function(q) we_cdf(q, bwe$alphaX, bwe$lambdaX))
ks_we_Y <- ks_fit(Y, function(q) we_cdf(q, bwe$alphaY, bwe$lambdaY))
AIC_we_X <- aic_val(bwe$loglik_X, 2); BIC_we_X <- bic_val(bwe$loglik_X, 2, nX)
AIC_we_Y <- aic_val(bwe$loglik_Y, 2); BIC_we_Y <- bic_val(bwe$loglik_Y, 2, nY)

ks_wr_X <- ks_fit(X, function(q) bwrr2_cdf(q, wrX$xi, wrX$b, wrX$k))
ks_wr_Y <- ks_fit(Y, function(q) bwrr2_cdf(q, wrY$xi, wrY$b, wrY$k))
AIC_wr_X <- aic_val(wrX$loglik, 3); BIC_wr_X <- bic_val(wrX$loglik, 3, nX)
AIC_wr_Y <- aic_val(wrY$loglik, 3); BIC_wr_Y <- bic_val(wrY$loglik, 3, nY)

cat("\n=== BWE (WE marginals via joint λ1,λ2,λ3) ===\n")
cat(sprintf("λ1=%.6f, λ2=%.6f, λ3=%.6f\n", bwe$lambda1, bwe$lambda2, bwe$lambda3))
cat(sprintf("X ~ WE(alpha=%.6f, lambda=%.6f): KS=%.4f (p=%.4f) | AIC=%.4f, BIC=%.4f\n",
            bwe$alphaX, bwe$lambdaX, ks_we_X$stat, ks_we_X$p.value, AIC_we_X, BIC_we_X))
cat(sprintf("Y ~ WE(alpha=%.6f, lambda=%.6f): KS=%.4f (p=%.4f) | AIC=%.4f, BIC=%.4f\n",
            bwe$alphaY, bwe$lambdaY, ks_we_Y$stat, ks_we_Y$p.value, AIC_we_Y, BIC_we_Y))
cat("\n=== BWRR2 (Weibull–Rayleigh marginals) ===\n")
cat(sprintf("X: xi=%.6f, b=%.6f, k=%.6f | KS=%.4f (p=%.4f) | AIC=%.4f, BIC=%.4f\n",
            wrX$xi, wrX$b, wrX$k, ks_wr_X$stat, ks_wr_X$p.value, AIC_wr_X, BIC_wr_X))
cat(sprintf("Y: xi=%.6f, b=%.6f, k=%.6f | KS=%.4f (p=%.4f) | AIC=%.4f, BIC=%.4f\n",
            wrY$xi, wrY$b, wrY$k, ks_wr_Y$stat, ks_wr_Y$p.value, AIC_wr_Y, BIC_wr_Y))

out <- rbind(
  data.frame(Model="BWE (WE via joint)", Margin="X",
             KS=as.numeric(ks_we_X$stat), p=ks_we_X$p.value,
             AIC=AIC_we_X, BIC=BIC_we_X),
  data.frame(Model="BWE (WE via joint)", Margin="Y",
             KS=as.numeric(ks_we_Y$stat), p=ks_we_Y$p.value,
             AIC=AIC_we_Y, BIC=BIC_we_Y),
  data.frame(Model="BWRR2", Margin="X",
             KS=as.numeric(ks_wr_X$stat), p=ks_wr_X$p.value,
             AIC=AIC_wr_X, BIC=BIC_wr_X),
  data.frame(Model="BWRR2", Margin="Y",
             KS=as.numeric(ks_wr_Y$stat), p=ks_wr_Y$p.value,
             AIC=AIC_wr_Y, BIC=BIC_wr_Y)
)

cat("\n--- Comparison Table ---\n"); print(out, row.names = FALSE)

<pre>An Error occurred while handling another error:
yii\web\HeadersAlreadySentException: Headers already sent in  on line 0. in /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Response.php:366
Stack trace:
#0 /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Response.php(339): yii\web\Response-&gt;sendHeaders()
#1 /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/ErrorHandler.php(136): yii\web\Response-&gt;send()
#2 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/ErrorHandler.php(135): yii\web\ErrorHandler-&gt;renderException()
#3 [internal function]: yii\base\ErrorHandler-&gt;handleException()
#4 {main}
Previous exception:
yii\web\HeadersAlreadySentException: Headers already sent in  on line 0. in /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Response.php:366
Stack trace:
#0 /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Response.php(339): yii\web\Response-&gt;sendHeaders()
#1 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/Application.php(656): yii\web\Response-&gt;send()
#2 /var/www/html/prof-homepages/vendor/faravaghi/yii2-filemanager/models/Files.php(696): yii\base\Application-&gt;end()
#3 /var/www/html/prof-homepages/vendor/faravaghi/yii2-filemanager/controllers/FilesController.php(484): faravaghi\filemanager\models\Files-&gt;getFile()
#4 [internal function]: faravaghi\filemanager\controllers\FilesController-&gt;actionGetFile()
#5 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/InlineAction.php(57): call_user_func_array()
#6 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/Controller.php(180): yii\base\InlineAction-&gt;runWithParams()
#7 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/Module.php(528): yii\base\Controller-&gt;runAction()
#8 /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Application.php(103): yii\base\Module-&gt;runAction()
#9 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/Application.php(386): yii\web\Application-&gt;handleRequest()
#10 /var/www/html/prof-homepages/backend/web/index.php(16): yii\base\Application-&gt;run()
#11 {main}</pre>