# نصب بسته‌ها (اگر نیاز باشد)
if(!require(evd)){install.packages("evd")} # بررسی می‌کند که آیا بسته 'evd' نصب شده است یا خیر. اگر نصب نشده باشد، آن را نصب می‌کند.
# بسته 'evd' برای کار با توزیع مقدار حدی تعمیم‌یافته (GEV) استفاده می‌شود.
if(!require(openxlsx)){install.packages("openxlsx")} # بررسی می‌کند که آیا بسته 'openxlsx' نصب شده است یا خیر. اگر نصب نشده باشد، آن را نصب می‌کند.
# بسته 'openxlsx' برای نوشتن داده‌ها در فایل اکسل استفاده می‌شود.

# فراخوانی کتابخانه‌ها
library(evd) # فراخوانی بسته 'evd' برای استفاده از توابع آن
library(openxlsx) # فراخوانی بسته 'openxlsx' برای استفاده از توابع آن

# تنظیم seed برای تکرارپذیری نتایج
set.seed(123) # تنظیم seed به عدد 123. این کار باعث می‌شود که هر بار کد اجرا شود، اعداد تصادفی یکسانی تولید شوند و نتایج قابل تکرار باشند.

# پارامترهای توزیع GEV (Generalized Extreme Value)
mu <- -0.35706 # پارامتر مکان (location) توزیع GEV
sigma <- 0.99171 # پارامتر مقیاس (scale) توزیع GEV
xi <- -0.26798 # پارامتر شکل (shape) توزیع GEV

# ماتریس انتقال (Transition Matrix)
# این ماتریس نشان می‌دهد که احتمال انتقال از هر حالت به حالت دیگر چقدر است.
# هر سطر نشان‌دهنده حالت فعلی و هر ستون نشان‌دهنده حالت بعدی است.
# به عنوان مثال، سطر اول نشان می‌دهد که اگر در حالت 1 باشیم، احتمال رفتن به حالت‌های 1 تا 7 چقدر است.
transition_matrix <- rbind(
  c(0.0512820512820513, 0.205128205128205, 0.025641026, 0.717948718, 0, 0, 0), # احتمالات انتقال از حالت 1 به حالات 1 تا 7
  c(0.183333333, 0.083333333, 0.05, 0.666666667, 0.016666667, 0, 0), # احتمالات انتقال از حالت 2 به حالات 1 تا 7
  c(0, 0.031496063, 0.165354331, 0.732283465, 0.05511811, 0.007874016, 0.007874016), # احتمالات انتقال از حالت 3 به حالات 1 تا 7
  c(0.026183283, 0.043303122, 0.085599194, 0.666666667, 0.123867069, 0.051359517, 0.003021148), # احتمالات انتقال از حالت 4 به حالات 1 تا 7
  c(0, 0, 0.023952096, 0.652694611, 0.137724551, 0.179640719, 0.005988024), # احتمالات انتقال از حالت 5 به حالات 1 تا 7
  c(0, 0, 0, 0.135416667, 0.510416667, 0.135416667, 0.21875), # احتمالات انتقال از حالت 6 به حالات 1 تا 7
  c(0.368421053, 0, 0, 0.631578947, 0, 0, 0) # احتمالات انتقال از حالت 7 به حالات 1 تا 7
)

num_states <- nrow(transition_matrix) # تعداد حالات (که برابر با تعداد سطرها یا ستون‌های ماتریس انتقال است)

# تابع چگالی احتمال GEV (Probability Density Function)
gev_pdf <- function(x) {
  evd::dgev(x, loc = mu, scale = sigma, shape = xi) # محاسبه چگالی احتمال توزیع GEV با استفاده از تابع 'dgev' از بسته 'evd'
}

# بازه‌های SPEI (Standardized Precipitation-Evapotranspiration Index)
# این لیست شامل بازه‌های SPEI برای هر یک از 7 حالت است.
spei_ranges <- list(
  c(-3, -2), # حالت 1: SPEI بین -3 و -2 (خشکسالی شدید)
  c(-2, -1.5), # حالت 2: SPEI بین -2 و -1.5 (خشکسالی متوسط)
  c(-1.5, -1), # حالت 3: SPEI بین -1.5 و -1 (خشکسالی خفیف)
  c(-1, 1), # حالت 4: SPEI بین -1 و 1 (نرمال)
  c(1, 1.5), # حالت 5: SPEI بین 1 و 1.5 (ترسالی خفیف)
  c(1.5, 2), # حالت 6: SPEI بین 1.5 و 2 (ترسالی متوسط)
  c(2, 3) # حالت 7: SPEI بین 2 و 3 (ترسالی شدید)
)

# تابع تولید SPEI تصادفی در بازه مشخص شده برای هر حالت
generate_spei_in_range <- function(state) {
  if (state < 1 || state > length(spei_ranges)) { # بررسی اینکه حالت ورودی معتبر است یا خیر
    stop(paste("Error: state must be between 1 and", length(spei_ranges))) # اگر حالت نامعتبر باشد، یک پیام خطا نمایش می‌دهد و اجرای کد متوقف می‌شود.
  }
  range <- spei_ranges[[state]] # دریافت بازه SPEI مربوط به حالت ورودی
  runif(1, range[1], range[2]) # تولید یک عدد تصادفی از توزیع یکنواخت بین مقادیر min و max بازه
}
# الگوریتم متروپلیس-هستینگز (با تغییرات)
metropolis_hastings_extended <- function(n_samples, initial_state) {
  samples <- numeric(n_samples + 1) # برداری برای ذخیره زنجیره مارکوف (دنباله‌ای از حالات) +1 برای ذخیره حالت اولیه
  all_proposed_states <- numeric(n_samples) # برداری برای ذخیره تمام حالات پیشنهادی در هر تکرار
  all_proposed_spei <- numeric(n_samples) # برداری برای ذخیره تمام SPEI های پیشنهادی در هر تکرار
  accepted_states <- numeric() # برداری برای ذخیره حالات پذیرفته شده توسط الگوریتم
  accepted_spei <- numeric() # برداری برای ذخیره SPEI های مربوط به حالات پذیرفته شده
  acceptance_indicator <- numeric(n_samples) # برداری برای ذخیره نشانگر پذیرش (1 برای پذیرش، 0 برای رد) در هر تکرار
  random_numbers <- numeric(n_samples) # برداری برای ذخیره اعداد تصادفی استفاده شده در تصمیم‌گیری پذیرش در هر تکرار
  acceptance_ratios <- numeric(n_samples) # برداری برای ذخیره نسبت‌های پذیرش محاسبه شده در هر تکرار
  
  current_state <- initial_state # تعیین حالت فعلی به عنوان حالت اولیه
  samples[1] <- current_state # ذخیره حالت اولیه در زنجیره مارکوف (اولین عنصر بردار samples)
  current_spei <- generate_spei_in_range(current_state) # محاسبه SPEI متناظر با حالت اولیه
  
  for (i in 1:n_samples) { # حلقه اصلی الگوریتم متروپلیس-هستینگز. این حلقه n_samples بار تکرار می‌شود.
    proposal_probabilities <- transition_matrix[current_state, ] # گرفتن احتمالات انتقال از حالت فعلی (current_state) از ماتریس انتقال
    proposed_state <- sample(1:num_states, 1, prob = proposal_probabilities) # انتخاب یک حالت پیشنهادی (proposed_state) با توجه به احتمالات انتقال. تابع sample یک نمونه تصادفی از اعداد 1 تا num_states با احتمالات مشخص شده در proposal_probabilities انتخاب می‌کند.
    proposed_spei <- generate_spei_in_range(proposed_state) # تولید یک مقدار SPEI تصادفی برای حالت پیشنهادی با استفاده از تابع generate_spei_in_range
    
    all_proposed_states[i] <- proposed_state # ذخیره حالت پیشنهادی در بردار all_proposed_states
    all_proposed_spei[i] <- proposed_spei # ذخیره SPEI پیشنهادی در بردار all_proposed_spei
    
    # محاسبه نسبت پذیرش (نسبت احتمال حالت پیشنهادی به حالت فعلی)
    # فرمول متروپلیس-هستینگز: min(1, (P(proposed) * Q(current|proposed)) / (P(current) * Q(proposed|current)))
    # در اینجا P همان gev_pdf و Q همان transition_matrix است.
    acceptance_ratio <- min(1, (gev_pdf(proposed_spei) * transition_matrix[proposed_state, current_state]) / (gev_pdf(current_spei) * transition_matrix[current_state, proposed_state]))
    acceptance_ratios[i] <- acceptance_ratio # ذخیره نسبت پذیرش در بردار acceptance_ratios
    
    random_number <- runif(1) # تولید یک عدد تصادفی از توزیع یکنواخت بین 0 و 1. این عدد برای تصمیم‌گیری در مورد پذیرش یا رد حالت پیشنهادی استفاده می‌شود.
    random_numbers[i] <- random_number # ذخیره عدد تصادفی تولید شده
    
    if (random_number < acceptance_ratio) { # شرط پذیرش: اگر عدد تصادفی تولید شده از نسبت پذیرش کمتر باشد، حالت پیشنهادی پذیرفته می‌شود.
      current_state <- proposed_state # به‌روزرسانی حالت فعلی به حالت پیشنهادی
      current_spei <- proposed_spei #به روز رسانی spei فعلی
      accepted_states <- c(accepted_states, current_state) # اضافه کردن حالت پذیرفته شده به بردار accepted_states
      accepted_spei <- c(accepted_spei, current_spei) # اضافه کردن SPEI پذیرفته شده به بردار accepted_spei
      acceptance_indicator[i] <- 1 # ثبت پذیرش (1) در بردار acceptance_indicator
    } else { # اگر شرط پذیرش برقرار نباشد، حالت پیشنهادی رد می‌شود و حالت فعلی بدون تغییر باقی می‌ماند.
      acceptance_indicator[i] <- 0 # ثبت رد (0) در بردار acceptance_indicator
    }
    
    samples[i + 1] <- current_state # ذخیره حالت فعلی (چه پذیرفته شده باشد چه رد شده باشد) در زنجیره مارکوف. این خط بسیار مهم است زیرا زنجیره مارکوف باید تمام حالات را شامل شود.
  }
  
  # خروجی تابع به صورت یک لیست شامل تمام داده‌های جمع‌آوری شده در طول اجرای الگوریتم
  list(samples = samples, # زنجیره مارکوف (دنباله‌ای از حالات نمونه‌برداری شده)
       all_proposed_states = all_proposed_states, # تمام حالات پیشنهادی
       all_proposed_spei = all_proposed_spei, # تمام SPEI های پیشنهادی
       accepted_states = accepted_states, # حالات پذیرفته شده
       accepted_spei = accepted_spei, # SPEI های پذیرفته شده
       acceptance_indicator = acceptance_indicator, # نشانگر پذیرش
       random_numbers = random_numbers, # اعداد تصادفی استفاده شده
       acceptance_ratios = acceptance_ratios) # نسبت‌های پذیرش
}
# اجرای الگوریتم
n_samples <- 10000 # تعداد نمونه‌هایی که توسط الگوریتم متروپلیس-هستینگز تولید می‌شود.
initial_state <- 1 # حالت اولیه زنجیره مارکوف. باید بین 1 و num_states باشد.

# بررسی اعتبار initial_state
if (initial_state < 1 || initial_state > length(spei_ranges)) { # بررسی می‌کند که آیا initial_state در محدوده معتبر قرار دارد یا خیر.
  stop(paste("Error: initial_state must be between 1 and", length(spei_ranges))) # اگر initial_state نامعتبر باشد، یک پیام خطا نمایش می‌دهد و اجرای کد متوقف می‌شود.
}

mh_results <- metropolis_hastings_extended(n_samples, initial_state) # اجرای الگوریتم متروپلیس-هستینگز با تعداد نمونه و حالت اولیه مشخص. نتیجه در متغیر mh_results ذخیره می‌شود.

# استخراج نتایج از لیست mh_results
all_proposed_states <- mh_results$all_proposed_states # استخراج تمام حالات پیشنهادی از خروجی تابع
all_proposed_spei <- mh_results$all_proposed_spei # استخراج تمام SPEI های پیشنهادی از خروجی تابع
accepted_states <- mh_results$accepted_states # استخراج حالات پذیرفته شده از خروجی تابع
accepted_spei <- mh_results$accepted_spei # استخراج SPEI های پذیرفته شده از خروجی تابع
acceptance_indicator <- mh_results$acceptance_indicator # استخراج نشانگرهای پذیرش از خروجی تابع
random_numbers <- mh_results$random_numbers # استخراج اعداد تصادفی استفاده شده در تصمیم‌گیری از خروجی تابع
acceptance_ratios <- mh_results$acceptance_ratios # استخراج نسبت‌های پذیرش محاسبه شده از خروجی تابع
sampled_states <- mh_results$samples[-1] #حذف مقدار اولیه از زنجیره نمونه برداری شده
# ایجاد دیتافریم‌ها برای ذخیره و سازماندهی نتایج
output_data_all <- data.frame(
  Iteration = 1:length(sampled_states), # شماره تکرار
  Proposed_State = all_proposed_states, # حالت پیشنهادی در هر تکرار
  Proposed_SPEI = all_proposed_spei, # SPEI پیشنهادی در هر تکرار
  Acceptance_Indicator = acceptance_indicator, # نشانگر پذیرش (1 برای پذیرش، 0 برای رد)
  Random_Number = random_numbers, # عدد تصادفی استفاده شده برای تصمیم‌گیری پذیرش
  Acceptance_Ratio = acceptance_ratios # نسبت پذیرش محاسبه شده
)

output_data_accepted <- data.frame(
  Accepted_State = accepted_states, # حالات پذیرفته شده
  Accepted_SPEI = accepted_spei # SPEI های پذیرفته شده
)

# ایجاد فایل اکسل برای ذخیره نتایج
wb <- createWorkbook() # ایجاد یک workbook جدید در اکسل
addWorksheet(wb, "All_Proposals") # اضافه کردن یک worksheet با نام "All_Proposals"
addWorksheet(wb, "Accepted") # اضافه کردن یک worksheet با نام "Accepted"
writeData(wb, sheet = "All_Proposals", x = output_data_all) # نوشتن داده‌های output_data_all در worksheet "All_Proposals"
writeData(wb, sheet = "Accepted", x = output_data_accepted) # نوشتن داده‌های output_data_accepted در worksheet "Accepted"
saveWorkbook(wb, "MH_results.xlsx", overwrite = TRUE) # ذخیره workbook در یک فایل اکسل با نام "MH_results.xlsx". اگر فایلی با این نام وجود داشته باشد، overwrite = TRUE آن را بازنویسی می‌کند.

# رسم نمودارها برای تجسم نتایج
hist(accepted_spei, main = "Histogram of Accepted SPEI Values", xlab = "SPEI") # رسم هیستوگرام SPEI های پذیرفته شده. main عنوان نمودار و xlab برچسب محور x است.
plot(1:length(sampled_states), sampled_states, type = "l", main = "Sampled States over Iterations", xlab = "Iteration", ylab = "State") # رسم نمودار خطی تغییرات حالات در طول تکرارها. type = "l" نمودار خطی رسم می‌کند.

# محاسبه و نمایش نرخ پذیرش
acceptance_rate <- sum(acceptance_indicator) / length(acceptance_indicator) # محاسبه نرخ پذیرش با تقسیم تعداد پذیرش‌ها بر کل تعداد تکرارها
cat("Acceptance Rate:", acceptance_rate, "\n") # نمایش نرخ پذیرش در کنسول