# تعریف مقادیر پارامترها
theta <- 0.0571
c2 <- 1000

# تعریف مقادیر c1
c1_values <- c(10000, 20000, 30000, 40000, 50000, 60000, 70000, 80000, 90000, 100000, 110000, 120000, 130000, 140000, 150000)

# ایجاد یک data.frame برای ذخیره نتایج نهایی
results_df <- data.frame()

# محاسبه برای p = 1 (توزیع دوجمله ای منفی)
cat("محاسبه دقیق برای p = 1 (توزیع دوجمله ای منفی)\n")

# پارامترهای توزیع دوجمله ای منفی
r <- 2  # تعداد موفقیت های مورد نیاز
p_success <- theta  # احتمال موفقیت

# تابع برای محاسبه دقیق تابع بقا (Survival function)
survival_negbin <- function(n, r, p) {
  # P(X > n) = 1 - P(X <= n)
  # برای توزیع دوجمله ای منفی با پارامترهای r و p
  1 - pnbinom(n - r, size = r, prob = p)
}

# تابع برای محاسبه دقیق امیدریاضی min(T, n)
expected_min_negbin <- function(n, r, p) {
  sum <- 0
  for (i in 0:(n-1)) {
    sum <- sum + survival_negbin(i, r, p)
  }
  return(sum)
}

# تابع برای محاسبه تابع چگالی احتمال (PDF)
pdf_negbin <- function(n, r, p) {
  # P(X = n) برای توزیع دوجمله ای منفی
  dnbinom(n - r, size = r, prob = p)
}

# حلقه روی مقادیر c1 برای p = 1
for (c1 in c1_values) {
  best_n <- NULL
  best_cost <- Inf
  
  # جستجو در محدوده n گسترده تر
  for (n in 1:100) {
    # محاسبه P(X > n) برای توزیع دوجمله ای منفی
    survival_n <- survival_negbin(n, r, p_success)
    
    # محاسبه E(min(X, n)) برای توزیع دوجمله ای منفی
    expected_min <- expected_min_negbin(n, r, p_success)
    
    # محاسبه هزینه C(n)
    if (expected_min > 0) {
      cost_n <- (c1 + (c2 - c1) * survival_n) / expected_min
      
      # بررسی آیا این هزینه بهترین است
      if (cost_n < best_cost) {
        best_cost <- cost_n
        best_n <- n
      }
    }
  }
  
  # ذخیره نتایج برای p = 1
  results_df <- rbind(results_df, data.frame(
    p = 1,
    c1 = c1,
    n_star = best_n,
    C_star = round(best_cost, 4)
  ))
  
  cat("برای c1 =", c1, ": n* =", best_n, ", C(n*) =", round(best_cost, 4), "\n")
}

# حالا محاسبه برای p = 0.75 (توزیع فاز-نوع)
cat("\nمحاسبه برای p = 0.75 (توزیع فاز-نوع)\n")
p <- 0.75

# محاسبه ضرایب بر اساس فرمول های مقاله
a2 <- (p * theta)^2
b1 <- theta * (1 + p) - 2
b2 <- 1 + p^2 * theta^2 - theta * (1 + p)

# محاسبه ریشه های معادله مشخصه
discriminant <- b1^2 - 4 * b2
lambda0 <- (-b1 - sqrt(discriminant)) / 2
lambda1 <- (-b1 + sqrt(discriminant)) / 2

# محاسبه omega0
p1 <- 0
p2 <- a2 - b1 * p1
omega0 <- (p2 - p1 * lambda1) / ((1 - lambda0) * (lambda0 - lambda1))
omega1 <- 1 - omega0

# حلقه روی مقادیر c1 برای p = 0.75
for (c1 in c1_values) {
  best_n <- NULL
  best_cost <- Inf
  
  # جستجو در محدوده n
  for (n in 1:50) {
    # محاسبه P(X > n) = ω0 * λ0^n + ω1 * λ1^n
    survival_n <- omega0 * (lambda0^n) + omega1 * (lambda1^n)
    
    # محاسبه E(min(X, n)) = Σ_{i=0}^{n-1} P(X > i)
    expected_min <- 0
    for (i in 0:(n-1)) {
      survival_i <- omega0 * (lambda0^i) + omega1 * (lambda1^i)
      expected_min <- expected_min + survival_i
    }
    
    # محاسبه هزینه C(n)
    if (expected_min > 0) {
      cost_n <- (c1 + (c2 - c1) * survival_n) / expected_min
      
      # بررسی آیا این هزینه بهترین است
      if (cost_n < best_cost) {
        best_cost <- cost_n
        best_n <- n
      }
    }
  }
  
  # ذخیره نتایج برای p = 0.75
  results_df <- rbind(results_df, data.frame(
    p = 0.75,
    c1 = c1,
    n_star = best_n,
    C_star = round(best_cost, 4)
  ))
  
  cat("برای c1 =", c1, ": n* =", best_n, ", C(n*) =", round(best_cost, 4), "\n")
}

# نمایش نتایج
print("نتایج محاسباتی:")
print(results_df)

#نمودارها
# ایجاد داده های جدول ۱
data <- data.frame(
  c1 = c(10000, 20000, 30000, 40000, 50000, 60000, 70000, 80000, 90000, 
         100000, 110000, 120000, 130000, 140000, 150000),
  n_star_p075 = c(17, 10, 7, 6, 5, 5, 4, 4, 4, 4, 3, 3, 3, 3, 3),
  C_star_p075 = c(148.22, 224.62, 280.24, 324.91, 363.76, 397.04, 428.20,
                  453.96, 479.72, 505.48, 526.90, 544.64, 562.38, 580.12, 597.86),
  n_star_p1 = c(11, 7, 5, 4, 4, 4, 3, 3, 3, 3, 3, 2, 2, 2, 2),
  C_star_p1 = c(203.18, 300.81, 370.94, 428.14, 473.61, 519.08, 550.66,
                582.10, 613.54, 644.99, 676.43, 694.29, 710.61, 726.94, 743.27)
)

# ذخیره نمودار در فایل PNG
png("optimal_replacement_plots.png", width = 1200, height = 600)

# نمودار ۱: تاثیر c1 بر N*
png("plot_Nstar.png", width = 800, height = 600)
plot(data$c1, data$n_star_p075, type = "o", col = "red", lwd = 2, pch = 16,
     xlab = "هزینه تعویض پیشگیرانه (c₁)", ylab = "تعداد بهینه تعویض (N*)",
     main = "تاثیر هزینه تعویض پیشگیرانه بر N*")
lines(data$c1, data$n_star_p1, type = "o", col = "blue", lwd = 2, pch = 17)
legend("topright", legend = c("p = 0.75", "p = 1"), 
       col = c("red", "blue"), pch = c(16, 17), lwd = 2)
dev.off()

# نمودار ۲: تاثیر c1 بر C(N*)
png("plot_Cstar.png", width = 800, height = 600)
plot(data$c1, data$C_star_p075, type = "o", col = "red", lwd = 2, pch = 16,
     xlab = "هزینه تعویض پیشگیرانه (c₁)", ylab = "هزینه بهینه (C(N*))",
     main = "تاثیر هزینه تعویض پیشگیرانه بر C(N*)")
lines(data$c1, data$C_star_p1, type = "o", col = "blue", lwd = 2, pch = 17)
legend("topleft", legend = c("p = 0.75", "p = 1"), 
       col = c("red", "blue"), pch = c(16, 17), lwd = 2)
dev.off()

cat("دو نمودار جداگانه ذخیره شدند: 'plot_Nstar.png' و 'plot_Cstar.png'\n")

getwd()

# نمایش فایل های موجود در پوشه کار
list.files()
<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>