#مقاله ی دوم طرح دوم محاسبه MLE وMPS
install.packages("nleqslv")
library(nleqslv)

# پارامترها
n <- 30
m <- 21
T <- 200
alpha <- 0.3
beta <- 0.6
lambda <- 0.75
# تابع برای تولید داده‌ها و حل معادلات
generate_estimates <- function() {
  # تولید n داده تصادفی از توزیع یکنواخت (0, 1)
  u <- runif(n)
  
  # تبدیل داده‌های یونیفرم به داده‌های توزیع پاور لومکس سه پارامتری
  x <- ((lambda^alpha / (1 - u))^(1/alpha) - lambda)^(1/beta)
  
  # نمایش داده‌ها
  print("تمام داده‌های تولید شده:")
  print(x)
  # مرتب سازی داده‌ها
  x_sorted <- sort(x)
  # بررسی شرایط سانسور
  if (x_sorted[m] < T) {
    # سانسور نوع اول
    x_censored <- x_sorted[1:m]
    R <- rep(0, m)
    R[2] <- n - m
    print("حالت نوع اول رخ داده است:")
    # تعریف معادلات
    equationsML <- function(par) {
      alpha <- par[1]
      beta <- par[2]
      lambda <- par[3]
      
      # تعریف معادلات بر اساس مشتقات جزئی
      eq1ML1 <- sum(sapply(1:m, function(i) {
        1/alpha + log(lambda) - 2*log(lambda + x_censored[i]^beta)+log(lambda) - log(lambda + x_censored[i]^beta)
      }))
      
      eq2ML1 <- sum(sapply(1:m, function(i) {
        1/beta + log(x_censored[i]) - (alpha + 1) * (x_censored[i]^beta * log(x_censored[i])) / (lambda + x_censored[i]^beta) -
          alpha * x_censored[i]^beta * log(x_censored[i]) / (lambda + x_censored[i]^beta)
      }))
      
      eq3ML1 <- sum(sapply(1:m, function(i) {
        alpha / lambda - (alpha + 1) / (lambda + x_censored[i]^beta) -
          alpha / (lambda + x_censored[i]^beta)
      }))
      return(c(eq1ML1, eq2ML1, eq3ML1))
    }
    # تعریف سه معادله 
    equationsPS <- function(params) {
      alpha <- params[1]
      beta <- params[2]
      lambda <- params[3]
      # تعریف توابع مشتق جزئی برای alpha, beta, lambda
      # معادله 1
      eq1PS1 <- sum(vapply(2:m, function(i) {
        log(lambda) + 
          (lambda^alpha * log(lambda) * ((lambda + x_censored[i-1]^beta)^(-alpha) - (lambda + x_censored[i]^beta)^(-alpha))) / 
          (lambda^alpha * ((lambda + x_censored[i-1]^beta)^(-alpha) - (lambda + x_censored[i]^beta)^(-alpha))) - 
          (lambda^alpha * ((lambda + x_censored[i-1]^beta)^(-alpha) * log(lambda + x_censored[i-1]^beta) - (lambda + x_censored[i]^beta)^(-alpha) * log(lambda + x_censored[i]^beta))) / 
          (lambda^alpha * ((lambda + x_censored[i-1]^beta)^(-alpha) - (lambda + x_censored[i]^beta)^(-alpha)))
      }, numeric(1)))
      
      # معادله 2
      eq2PS1 <- sum(vapply(2:m, function(i) {
        (lambda^alpha * ((lambda + x_censored[i-1]^beta)^(-alpha) * log(lambda + x_censored[i-1]) - (lambda + x_censored[i]^beta)^(-alpha) * log(lambda + x_censored[i]))) / 
          (lambda^alpha * ((lambda + x_censored[i-1]^beta)^(-alpha) - (lambda + x_censored[i]^beta)^(-alpha))) -
          alpha * x_censored[i]^beta * log(x_censored[i]) / (lambda + x_censored[i]^beta)
      }, numeric(1)))
      
      # معادله 3
      eq3PS1 <- sum(vapply(2:m, function(i) {
        alpha / lambda +
          (lambda^(alpha-1) * ((lambda + x_censored[i-1]^beta)^(-alpha) - (lambda + x_censored[i]^beta)^(-alpha))) / 
          (lambda^alpha * ((lambda + x_censored[i-1]^beta)^(-alpha) - (lambda + x_censored[i]^beta)^(-alpha))) - 
          (lambda^alpha * (-alpha * (lambda + x_censored[i-1]^beta)^(-alpha-1) + alpha * (lambda + x_censored[i]^beta)^(-alpha-1))) / 
          (lambda^alpha * ((lambda + x_censored[i-1]^beta)^(-alpha) - (lambda + x_censored[i]^beta)^(-alpha))) +
          alpha / lambda - alpha / (lambda + x_censored[i]^beta)
      }, numeric(1)))
      
      
      # بازگشت مقادیر معادلات
      return(c(eq1PS1,eq2PS1,eq3PS1))
    }
    
    # مقدارهای اولیه برای پارامترها
    initial_guess <- c(0.3, 0.6, 0.75)
    # حل معادلات
    solution1 <- nleqslv(initial_guess, equationsML)
    solution2 <- nleqslv(initial_params, equationsPS)
    #solution2 <- nleqslv(x = initial_params, fn = equationsPS)
    # نمایش نتایج
    alphaML <- solution1$x[1]
    betaML <- solution1$x[2]
    lambdaML <- solution1$x[3]
    # نتایج
    alphamps <- solution2$x[1]
    betamps <- solution2$x[2]
    lambdamps <- solution2$x[3]
  }
  else {
    # بررسی شرایط حالت نوع دوم
    h <- sum(x_sorted < T)
    if (h < m && x_sorted[h + 1] > T) {
      x_censored <- x_sorted[1:h]
      R <- rep(0, h)
      if (h >= 2) {
        R[2] <- n - m
        print("حالت دوم رخ داده است:")
        # تعریف معادلات
        equationsML <- function(par) {
          alpha <- par[1]
          beta <- par[2]
          lambda <- par[3]
          
          h <- length(x_censored)
          R_h_star <- R[h]
          
          eq1ML2 <- h / alpha + h * log(lambda) - sum(log(lambda + x_censored^beta) * (1 + R)) +
            log(lambda) * sum(R) + R_h_star * log(lambda) - R_h_star * log(lambda + T^beta)
          
          eq2ML2 <- h / beta + sum(log(x_censored)) - sum((log(x_censored) * x_censored^beta * ((alpha + 1) + alpha * R)) / (lambda + x_censored^beta)) -
            alpha * R_h_star * (T^beta * log(T)) / (lambda + T^beta)
          
          eq3ML2 <- h * alpha / lambda - sum(((alpha + 1) + alpha * R) / (lambda + x_censored^beta)) -
            alpha * R_h_star / (lambda + T^beta) + alpha / lambda * sum(R) + alpha / lambda * R_h_star
          
          return(c(eq1ML2, eq2ML2, eq3ML2))
        }
        
        equationsPS <- function(params) {
          alpha <- params[1]
          beta <- params[2]
          lambda <- params[3]
          
          h = length(x_censored)
          R_h_star = R[h]
          # معادله اول 
          eq1PS2 <- (lambda^alpha * (log(lambda + x_censored[1]) - log(lambda))) / ((lambda + x_censored[1])^alpha - lambda^alpha) -
            log(lambda + x_censored[1]) + 2 * log(lambda) - R_h_star * log(lambda + T) +
            log(lambda) * sum(R_h_star) - sum(R_h_star * log(lambda + x_censored[1])) +
            sum(sapply(2:length(x_censored), function(i) {
              ((lambda + x_censored[i])^(-alpha) * log(lambda + x_censored[i]) -
                 (lambda + x_censored[i-1])^(-alpha) * log(lambda + x_censored[i-1])) /
                ((lambda + x_censored[i])^(-alpha) - (lambda + x_censored[i-1])^(-alpha))
            }))
          # معادله دوم 
          eq2PS2 <- (alpha * lambda^alpha * x_censored[1]^beta * log(x_censored[1])) /
            (1 - lambda^alpha * (lambda + x_censored[1])^(-alpha)) -
            (alpha * x_censored[1]^beta * log(x_censored[1])) / (lambda + x_censored[1]) -
            alpha * sum(R_h_star * log(x_censored) * x_censored^beta / (lambda + x_censored)) -
            alpha * R_h_star * T^beta * log(T) / (lambda + T) -
            sum(sapply(2:length(x_censored), function(i) {
              (x_censored[i]^beta * log(x_censored[i]) * (lambda + x_censored[i])^(-alpha-1) -
                 x_censored[i-1]^beta * log(x_censored[i-1]) * (lambda + x_censored[i-1])^(-alpha-1)) /
                ((lambda + x_censored[i])^(-alpha) - (lambda + x_censored[i-1])^(-alpha))
            }))
          # معادله سوم 
          eq3PS2 <- (-alpha * lambda^(alpha-1) * (lambda + x_censored[1])^(-alpha-1) * (lambda + x_censored[1] - lambda)) /
            (1 - lambda^alpha * (lambda + x_censored[1])^(-alpha)) +
            alpha * sum(sapply(2:length(x_censored), function(i) {
              ((lambda + x_censored[i])^(-alpha-1) - (lambda + x_censored[i-1])^(-alpha-1)) /
                ((lambda + x_censored[i])^(-alpha) - (lambda + x_censored[i-1])^(-alpha))
            })) -
            alpha / (lambda + x_censored[length(x_censored)]) + 2 * alpha / lambda -
            alpha * sum(R_h_star / (lambda + x_censored)) -
            alpha * R_h_star / (lambda + T) + alpha * sum(R_h_star) / lambda + alpha * R_h_star / lambda
          
          return(c(eq1PS2,eq2PS2,eq3PS2))
        }
        
        # مقدارهای اولیه برای پارامترها
        initial_guess <- c(0.3, 0.6, 0.75)
        
        # حل معادلات
        solution1 <- nleqslv(initial_guess, equationsML)
        solution2 <- nleqslv(initial_params , equationsPS)
        #solution2 <- suppressWarnings(nleqslv(initial_params , equationsPS))
        # نمایش نتایج
        alphaML <- solution1$x[1]
        betaML <- solution1$x[2]
        lambdaML <- solution1$x[3]
        # نتایج
        alphamps <- solution2$x[1]
        betamps <- solution2$x[2]
        lambdamps <- solution2$x[3] 
      } 
    }
   }
    
    # نمایش داده‌های سانسور شده و مقادیر R
    print("داده‌های سانسور شده:")
    print(x_censored)
    print("مقادیر R:")
    print(R)
    print(c(alphaML,
            betaML ,
            lambdaML,alphamps,
            betamps ,
            lambdamps))
    
    # فیلتر کردن نتایج غیرمنطقی
    if (alphaML > 0 && betaML > 0 && lambdaML > 0 && alphaML < 10 && betaML < 10 && lambdaML < 10 && alphamps > 0
        && betamps > 0 && lambdamps > 0 && alphamps < 10 && betamps < 10 && lambdamps < 10) {
      return(c(alphaML, betaML, lambdaML,alphamps, betamps, lambdamps))
    } else {
      return(c(NA, NA, NA,NA, NA, NA))
    }
  }
  
  # تکرار فرآیند 10000 بار و سرکوب اخطارها
  results <- suppressWarnings(replicate(10000, generate_estimates()))
  
  # تبدیل نتایج به یک data.frame برای تحلیل
  b <- as.data.frame(t(results))
  names(b) <- c("alpha_hat", "beta_hat", "lambda_hat","alpha_hatPS", "beta_hatPS", "lambda_hatPS")
  
  # حذف مقادیر NA
  b <- na.omit(b)
  b
  alpha_hatML = b[,1]
  alpha_hatML 
  beta_hatML = b[,2]
  beta_hatML
  lambda_hatML = b[,3]
  lambda_hatML
  
  alpha_hatPS = b[,4]
  alpha_hatPS 
  beta_hatPS = b[,5]
  beta_hatPS
  lambda_hatPS = b[,6]
  lambda_hatPS
  
  mse_alphahatML = mean((alpha_hatML - alpha)^2)
  mse_alphahatML
  mse_betahatML = mean((beta_hatML - beta)^2)
  mse_betahatML
  mse_lambdahatML = mean((lambda_hatML - lambda)^2)
  mse_lambdahatML
  
  EB_alphahatML = mean(alpha_hatML) - alpha
  EB_alphahatML
  EB_betahatML = mean(beta_hatML) - beta
  EB_betahatML 
  EB_lambdahatML = mean(lambda_hatML) - lambda
  EB_lambdahatML
  
  
  mse_alphahatPS = mean((alpha_hatPS - alpha)^2)
  mse_alphahatPS
  mse_betahatPS = mean((beta_hatPS - beta)^2)
  mse_betahatPS
  mse_lambdahatPS = mean((lambda_hatPS - lambda)^2)
  mse_lambdahatPS
  EB_alphahatPS = mean(alpha_hatPS) - alpha
  EB_alphahatPS
  EB_betahatPS = mean(beta_hatPS) - beta
  EB_betahatPS 
  EB_lambdahatPS = mean(lambda_hatPS) - lambda
  EB_lambdahatPS
  
  <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/frontend/web/index.php(18): yii\base\Application-&gt;run()
#11 {main}</pre>